THE UNIVERSITY OF MICHIGAN COLLEGE OF ENGINEERING Department of Meteorology and Oceanography Technical Report No. 1 BAROCLINIC INSTABILITY OF ULTRA-LONG WAVES Perry W. Fisher Aksel C. Wiin-Nielsen Project Director.ORA. Project 00263 supported by: NATIONAL SCIENCE FOUNDATION GRANT NO. GA-16166 WASHINGTON, D.C. administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR November 1970

/CA' 4f e r4

TABLE OF CONTENTS Page LIST OF TABLES v LIST OF FIGURES vi LIST OF SYMBOLS x ABSTRACT xvi Chapter I. INTRODUCTION I 1,1. Background and Purpose of the Study 1 1.2. Brief Review of Previous Work 4 1-3. Outline of the Study.11 II. FORMULATION AND METHOD OF SOLUTION 13 2.1. Formulation 13 2.1.1 Assumptions and basic equations 13 2.1.2. Linearization and reduction of the basic equations 15 2.1,3. Determination of an exact solution 23 2.1.4. Derivation of a semi-circle theorem 27 2.1.5. Derivation of a sufficient condition for stability 33 2.1.6. The Newtonian heating model 36 2.2. Method of Solution 41 III. PRESENTATION AND DISCUSSION OF RESULTS 45 3.1. The Stability Analysis 45 3.2. The Phase Speed 66 3.3. The Wave Structure 71 IV. ENERGETICS 85 4.1. General Considerations 85 4.2. Conversion from Zonal to Eddy Available Potential Energy, C(AzAE) 87 4.3. Conversion from Eddy Available Potential to Eddy Kinetic Energy, C(AE,KE) 100 4.4. Comparison of Energy Strip and Full Energy Diagrams 111 iii

TABLE OF CONTENTS (Concluded) Page V. CONCLUSION 118 5.1. Summary of Results 118 5.2. Concluding Remarks and Suggestions for Future Work 123 APPENDIX 125 A.1. The Barotropic Limiting Case, A = A0o(p), az = ao = Constant 125 A.2. The Barotropic Limiting Case, A = A0(cp), z = o/P* 127 A.3. Comparison of the Quasi-Geostrophic and Geostrophic Models for a Basic State of No Motion and Constant Static Stability 129 A.4. Discussion of the Uniqueness of the Unstable Efgenvalue Solutions 138 A.5. Relationship Between Westward Wave Tilt and Northward Sensible Heat Transport 145 BIBLIOGRAPHY 150 iv

LIST OF TABLES Table Page Al. Difference Between Basic Current and Wave Speed for Barotropic Limiting Case, Cz = ao = Constant 126 A2. Difference Between Basic Current and Wave Speed for Barotropic Limiting Case, cz = fo/P* 128 A3. Wave Speed in Quasi-Geostrophic Model for Basic State of No Motion and Constant Static Stability 156 A4. r,,o Calculated From Numerical Unstable Solutions 142 v

LIST OF FIGURES Figure Page 1. Wave length as a fucr:iion of latitude for wave numbers 1, 2, and 3. 46 2. Zonal wind at p = 0 as a function of latitude as calculated for a linear profile by doubling 500 mb data of January, July, and the annual average, 1963. 47 3. Imaginary component of phase velocity as a function of latitude for adiabatic model-January, July, and annual average, 1963. 50 4. e-Folding time as a function of latitude for adiabatic model, n = 1-January, July, and annual average, 1963. 51 5. Imaginary component of phase velocity as a function of zonal wind at p = 0 for adiabatic model-January, cp = 45~; July, cp = 30~ and 60~; and annual average,c = 75~, 1963. 53 6. e-Folding time as a function of latitude for January, 1963adiabatic model at n = 1, 2, and 3. 55 7. e-Folding time as a function of latitude for July, 1963adiabatic model at n = 1, 2, and 3. 56 8. e-Folding time as a function of latitude for annual, 1963adiabatic model at n = 1, 2, and 3. 57 9. Imaginary component of phase velocity as a function of latitude for January, 1963-adiabatic model; and Newtonian model at n = 1, 2, and 3. 60 10. Imaginary component of phase velocity as a function of latitude for July, 1963 —adiabatic model; and Newtonian model at n = 1, 2, and 3. 61 11. Imaginary component of phase velocity as a function of latitude for annual average, 1963-adiabatic model; and Newtonian model at n = 1, 2, and 3. 62 12. Imaginary component of phase velocity as a function of latitude for Newtonian model at n = 1. 63 vi

LIST OF FIGURES (Continued) Figure Page 13. Imaginary component of phase velocity as a function of latitude for Newtonian model at n = 2. 64 14. Imaginary component of phase velocity as a function of latitude for Newtonian model at n = 3. 15. Phase speed as a function of latitude for adiabatic and Newtonian models-January, July, and annual average, 1963. 68 16. Phase angle of perturbation geopotential of most unstable wave as a function of pressure for adiabatic and Newtonian models. 76 17. Normalized amplitude of perturbation geopotential of most unstable wave as a function of pressure for adiabatic and Newtonian models. 78 18. Phase angle of perturbation vertical (pressure) velocity of most unstable wave as a function of pressure for adiabatic and Newtonian models. 80 19. Normalized amplitude of perturbation vertical (pressure) velocity of most unstable wave as a function of pressure for adiabatic and Newtonian models. 82 20. Phase angle of perturbation temperature of most unstable wave as a function of pressure for adiabatic and Newtonian models. 83 21. Normalized amplitude of perturbation temperature of most unstable wave as a function of pressure for adiabatic and Newtonian models. 84 22. Normalized conversion from zonal to eddy available potential energy as a function of pressure for January, 1963-adiabatic model at n = 1, 2, and 3. 90 23. Normalized conversion from zonal to eddy available potential energy as a function of pressure for January, 1963 —Newtonian model at n = 1, 2, and 3. 91 24. Normalized conversion from zonal to eddy available potential energy as a function of pressure for July, 1963 —adiabatic model at n = 1, 2, and 3. 93 vii

LIST OF FIGURES (Continued) Figure Page 25. Normalized conversion fromzonal to eddy available potential energy as a function of pressure for July, 1963-Newtonian model at n = 1, 2, and 3. 94 26. Normalized conversion from zonal to eddy available potential energy as a function of pressure for annual, 1963-adiabatic model at n = 1y 2, and 3. 96 27. Normalized conversion from zonal to eddy available potential energy as a function of pressure for annual, 1963 —Newtonian model.at n = 1) 2, and 3. 97 28. Normalized conversion from zonal to eddy available potential energy as a function of pressure for adiabatic model at n =. 98 29. Normalized conversion from zonal to eddy available potential energy as a function of pressure for Newtonian model at n = 1. 99 30. Normalized conversion from eddy available potential energy to eddy kinetic energy as a function of pressure for January, 1963-adiabatic model at n = 1, 2, and 3. 103 31. Normalized conversion from eddy available potential energy to eddy kinetic energy as a function of pressure for January, 1963-Newtonian model at n = 1, 2, and 3. 104 32. Normalized conversion from eddy available potential energy to eddy kinetic energy as a function of pressure for July, 1963adiabatic model at n = 1, 2, and 3. 106 33. Normalized conversion from eddy available potential energy to eddy kinetic energy as a function of pressure for July, 1963Newtonian model at n = 1, 2, and 3. 107 34. Normalized conversion from eddy available potential energy to eddy kinetic energy as a function of pressure for annual, 1963adiabatic model at n = 1, 2, and 3. 108 35. Normalized conversion from eddy available potential energy to eddy kinetic energy as a function of pressure for annual, 1963 — Newtonian model at n = 1, 2, and 3. 109 viii

LIST OF FIGURES (Concluded) Figure Page 36. Normalized conversion from eddy available potential energy to eddy kinetic energy as a function of pressure for adiabatic model at n = 1. 112 37. Normalized conversion from eddy available potential energy to eddy kinetic energy as a function of pressure for Newtonian model at n = 1. 113 ix

LIST OF SYMBOLS a = mean radius of the earth except in Section 2.1.4, 6371 km a = minimum value of A at given cp in Section 2.1.4 b = maximum value of A at given cp in Section 2.1.4 c = phase velocity of the perturbations c = c for adiabatic model *a * c = c for Newtonian heating model N = s h c = specific heat of air at constant pressure p c = specific heat of air at constant volume v d = upper boundary of y in Appendix A.3 f = Coriolis parameter except in Appendix A.4, 2Qsincp 2 rf = o o + ccosc fJ U sinp'dcp' in Appendix A.4 2 2 2 A d 2QA a. sin cp U sin c L o o -4 f = Coriolis parameter f evaluated at central latitude, o, 10 0o1 o sec -2 g = acceleration of gravity, 9-8 msec h = equivalent depth k = 2T/(wave length) m = arbitrary positive or negative integer n = zonal wave number except in Sections 2.1.3 and 2.2 n = arbitrary positive integer greater than 1 in Sections 2.1.3 and 2.2 p = pressure p = standard pressure, usually 1000 mb x

LIST OF SYMBOLS (Continued) 2 2 2 q = a p /2.A a sin (p in Section 2.1.3 o00 oe e q = heating or cooling coefficient in Section 2.1.6 r = 41 in Section 2.1.3 r = iT in Appendix A.4 r = value of r around which Markushevich's implicit function theorem is applied in Appendix A.4 s = +1 for A > 0, -1 for A < 0 o o t = time u = component of speed along the x-axis v = component of speed along the y-axis w = vertical component of speed x = distance coordinate increasing to the east y = distance coordinate increasing to the north z = height above mean sea level A = available potential energy C(AE K) = conversion from eddy available potential to eddy kinetic E E energy C(A,AE) = conversion from zonal to eddy available potential energy F = p-dependent function in $ in Appendix A.3 2 1(' = Gauss' function G = normalized W at p = 1 in Chapter II G = y-dependent function in t in Appendix A.3 H = vertical scale of motion in Chapter II xi

LIST OF SYMBOLS (Continued) H = northward sensible heat transport averaged over the latitude circle at p = p in Appendix A.5 K = kinetic energy L = horizontal scale of motion M = A-c N =- 1 ^ ^ 2 Q -= | a /IMl in Section 2.1.4 Q = diabatic heating per unit time and unit mass in Section 2.1.6 -1 -l R = gas constant for air except in Section 2.1.4, 287 kj t deg R = sI |"12I /IMi in Section 2.1.4 Ri = Richardson number Ro = Rossby number 2 2 2 S a== p /2Ra sin p in Chapter II S = surface area in Chapter IV T = temperature T = e-folding time e TE = convective-radiative equilibrium temperature except in Chapter IV T = eddy or perturbation temperature in Chapter IV E U = velocity component (in the basic state) along the x-axis V = characteristic horizontal wind speed W = characteristic vertical wind speed on first page of Chapter II W = M 1/2G ^ in the rest of Chapter II xii

LIST OF SYMBOLS (Continued) a = specific volume in Chapter II a = phase angle of temperature on p = pi in Appendix A.5 -12 D3 ~= df/dy, the Rossby parameter, except in Appendix A.5, 16 x 10 m-1 sec" 3B fi= phase angle of geopotential on p = p in Appendix A.5 56 =b= maximum allowable G in eigenvalue computations in Chapter II 6 = phase angle in Chapters III and V C = cVc, 5/7 for air Q = potential temperature x\ -= wave length in Chapters I and III - = longitude in Chapters II and IV = p*/(1-sc.) e -= 5 at p = 1 i^5= -value of ~o around which Markushevich's implicit function theorem is applied in Appendix A.4 p = density except in Appendix A.4 p = radius of 0o-neighborhood of analyticity in Markushevich's implicit function theorem in Appendix A.4 p' = radius of 0o-neighborhood in which Markushevich's implicit function theorem is applicable in Appendix A.4 a = -aalnO/.p (where a is the specific volume), a measure of the static stability o = GzL at p = I CP = latitude xiii

LIST OF SYMBOLS (Continued) Jt ~ = radius of r-neighborhood of analyticity in Markushevich's implicit function theorem in Appendix A.4 t~' == radius of r-neighborhood in which Markushevich's implicit function theorem is applicable in Appendix A.4 X = vertical (pressure) velocity, dp/dt 2 22 A = f k c in Appendix A.3 1Ap - P2 in Appendix A.5 A = basic state zonal angular velocity b ^ = gz, geopotential -5 ]-1 0.D = angular velocity of the earth, 7.3 x 10 sec ~( ) = the steady, axially symmetric basic state component of ( ) in Chapter II ( ) =,the zonal average of ( ) in Chapter IV z ( )' = ( ) - ( ) = the perturbation component of ( ) in Chapter II z ()E - = ( ) - ( ) = the eddy component of ( ) in Chapter IV A ( ) = that part of ( )' which depends on the latitude and pressure only ( ) = area average of ( ) over a given isobaric surface ( " = ( ) - () = the deviation of ( ) from its area average = complex conjugate of ( ) ( ) = the real component of ( ) ( ). = the imaginary component of ( ) ( ) = nondimensional ( ) ( )L = ( ) at the lower latitude limit, PL xiv

LIST OF SYMBOLS (Concluded) ( )e = ( ) for exact solution, except T e e ( )o = ( )at p = 0, except p, a, f, and ( )1 = ( ) at p = 1 in Chapter II ( ) = ( ) at p = Pi in Appendix A.5 ( )2 = ( ) at p = P2 in Appendix A.5 ( ) = maximum value of ( )I at a given cp for all p max ( ) m = minimum value of ( )I at a given cp for all p min ( )N = normalized ( ) at a given cp for all p N ( )N = normalized ( ) for all cp and p''~~~NO ~~ - xv

ABSTRACT Baroclinic instability of ultra-long waves (wave numbers 1, 2, and 3) is studied with a 40-layer, geostrophic, hydrostatic, adiabatic model which is later modified to include Newtonian heating. The assumption of a meridionallyvarying zonal current linear in pressure upon which are superimposed wavelike perturbations traveling in the west-east direction and static stability dependent only upon the inverse of pressure allows the reduction of the basic equations to a single second-order homogeneous ordinary differential equation in vertical pressure velocity having latitude as an implicit parameter. The "shooting method, a numerical search procedure to determine the eigenvalue, perturbation phase velocity, is employed under the simplifying boundary conditions that vertical pressure velocity is zero at the top and bottom of the atmosphere. This method requires one exact solution, and a procedure is given to determine it. A discussion is also included of the uniqueness of the unstable eigenvalues. Three separate cases for both the adiabatic and Newtonian models are considered, based upon zonal wind means of January, July, and the annual average, 1963 from 20~N to 85~N. Investigation of instability as measured by the imaginary component of phase velocity reveals no dependence upon wave number for adiabatic flow, while instability increases as wave number increases for Newtonian heating, which itself exerts a stabilizing effect, becoming more pronounced with increasing instability. At all wave lengths in adiabatic and Newtonian flow greater instability is found over most latitudes in January, than in July or for the xvi

annual average, having a maximum at 25~N in January, at 40"N in July, and at 35~N for the annual average. In the adiabatic model this corresponds to efolding times for the third harmonic of 4, 7 1/2, and 6 days, respectively. The phase speed, defined as the real part of the phase velocity of the perturbation wave, is the same for both adiabatic and Newtonian flow and is independent of wave number. Values of the phase speed for the most unstable waves are 8.6, 5.2, and 6.3 m/sec in January, July, and for the annual average, respectively. Vertical wave structure of the most unstable waves is investigated and is found strikingly similar for the adiabatic and Newtonian models. Virtually no variation with wave number or time is observed. Phase angle and normalized amplitude variation with pressure of perturbation geopotential, vertical pressure velocity, and temperature reveal a westward tilt with decreasing pressure, Newtonian phase angles 5~ west of their adiabatic counterparts at every level, and negligible difference in normalized amplitude between the models. Instantaneous energetics of the adiabatic and Newtonian models is studied (1) by calculation of the vertica:L variation and total over the entire atmosphere of the normalized conversion from zonal to eddy available potential energy, C(A,AE), and normalized conversion from eddy available potential z E energy to eddy kinetic energy, C(AE,KE), and (2) by calculation of the vertical variation and total in a. latitudinal strip 100 wide centered on the latitude of maximum instability of the same energy conversions. In this strip it is found that, with the sole exception of C(AE,KE) < 0 in the January-Newtonian model, C(A,AE) > 0 and C(A,KE) > 0 which agrees with the time-averaged z xvii xvii

energetics as observed in the atmosphere. This is not so for the energy conversions over all latitudes, where the energetics of the stable perturbations apparently dominate that of the unstable perturbations and the models' validity may be questioned in the high-latitude regime. Here, however, addition of Newtonian heating does cause the instantaneous energetics to conform more closely to the time-averaged empirical energy conversions, providing some support for the inclusion of Newtonian heating in an ultra-long wave model. xviii

CHAPTER I INTRODUCTION 1.1. BACKGROUND AND PURPOSE OF THE STUDY On hemispherical weather charts atmospheric motion can be described as zonal (west-east) flow on which are superimposed disturbances. These disturbances represent the superposition of components having horizontal scales in a very wide range. On hemispheric 500 mb maps of mid-latidudinal westerly flow one clearly observes the existence of very long waves having wave lengths from 10,000-20,000 km. In this study ultra-long waves are synonymous with very long waves, specifically those waves having wave numbers 1, 2, or 3, where wave number indicates the number of waves of given wave length around a latitude circle. While some authors use the term "planetary" waves to describe waves with length comparable to the planet's radius, i.e., ultralong waves, others use it to describe any wave occurring in a planetary atmosphere. To avoid confusion we shall refrain from using it in this study. The observed ultra-long waves are either stationary or slow-moving. Eliasen (1958) describes these waves as slowly-moving systems whose amplitude doubles in a few days, oscillating around certain preferred geographical positions. For the purpose of mathematically modeling the ultra-long waves we consider a superposition of two distinct wave types-the stationary or quasi-stationary waves and the transient, traveling, or free waves. The former necessarily result from constant large-scale external forcing such as terrain or heat sources and sinks, while the latter may be subdivided into barotropic and 1

2 baroclinic modes, the barotropic describing the vertically-averaged component and the baroclinic the deviation from the vertical average. It is the purpose of this study to show how a baroclinic, adiabatic, frictionless, geostrophic model describes the transient ultra-long wave dynamics of the atmosphere and how the addition of linear diabatic heating modifies the solution. Analysis of our models with regard to the stability and structure of ultra-long waves reveals the existence of unstable solutions. Thus, the presence of baroclinic instability in the atmosphere may be a factor accounting for the fact that ultra-long waves often seem to be excited. Mathematically, the hydrodynamic instability problem is as follows: Suppose the system of hydrodynamic equations has a steady-state solution for the components of velocity, pressure, and temperature. We consider an initialvalue problem with these variables slightly different from those in this timeindependent solution. If the solution approaches this steady-state solutior as time t -+ o, the motion is stable. Otherwise, it is unstable. Baroclinic instability is a hydrodynamic instability arising from the existence of a meridional temperature gradient (and thus, a thermal wind) in the atmosphere in which potential energy of the basic flow is converted into kinetic energy of the unstable perturbation. We distinguish this from barotropic instability, a hydrodynamic instability arising from certain distributions of relative vorticity in two-dimensional flow in which kinetic energy is the only form of of energy transferred between the basic current and perturbation. Baroclinic instability is considered to be one of the factors responsible for the development of wave disturbances within strong westerly wind flow which frequently

occurs in middle and high latitudes. Growth of the disturbances is characterized by ascent of warmer, and descent of colder, air masses, representing a decrease in disturbance total potential energy and its associated release of kinetic energy. In order to formulate the mathematical problem of hydrodynamic instability we take, from the many possible steady-state solutions of the equations of motion, one which is mathematically expedient and superimpose upon it a disturbance of a suitable kind, generating a set of nonlinear disturbance equations governing the behavior of the disturbance. By assuming the disturbance or perturbation amplitude small, i.e., the square of the amplitude is negligible in comparison with the amplitude, and derivatives of the perturbation of the same order of magnitude as the perturbation, we employ the method of small perturbations to linearize the disturbance equations. This eliminates the nonlinear advective terms in the perturbation hydrodynamic equations. The resultant linear system of equations contains time t only through derivatives with respect CTt to t, and hence solutions containing an exponential time factor e may be expected. Boundary conditions on the disturbance equations vary considerably with the nature of the problem but usually require vanishing of at least one of the perturbation velocity components at the boundaries. In general the boundary conditions, like the disturbance equations, are homogeneous, i.e., each term is proportional to a dependent variable or one of its derivatives, and we therefore have an eigenvalue problem for the determination of a. If a has a positive real component, the flow is unstable according to linear theory, the greatest instability occurring at the maximum value of a; otherwise, the

flow is stable. When a is known, the structure of the perturbations can usually then be determined. An unstable nonlinear system may or may not approach another steady state. The method of small perturbations is incapable of making this prediction. Despite this and other limitations the method outlined is the most common way of treating an instability problem and has met with considerable success. Cowling (1957) has pointed out: "The perturbation method is not altogether adequate, particularly in problems of the stability of flow: flow is sometimes stable for small disturbances but unstable for large. But no method superior to the perturbation method has yet been devised" (P. 57). 1.2. BRIEF REVIEW OF PREVIOUS WORK In 1939 Rossby developed a theoretical treatment of the motion of long waves. He assumed uniform zonal flow and wave perturbations independent of latitude in a barotropic atmosphere. To simplify the problem he introduced the:-plane in which the only effect of the earth's curvature considered appears in the north-south variation of the vertical component of the earth's rotation. This model gives stable waves moving at speed 2 c - U - 2 X 4n2 where U is the basic current, X is the wave length, and v = 2Qcosp/a. Here, a is the earth's radius, Q the earth's rotation rate, and p the latitude. In comparing his results with observations we find reasonable agreement for moderately long waves. His predictions for the very long waves, however, of strong

retrograde motion and c highly sensitive to small changes in wave length are not confirmed. In delineating behavior of very long waves, Eliasen (1958) describes them from zonal harmonic analysis as slowly moving systems which oscillate around certain preferred geographical positions. Haurwitz (1940) modified Rossby's model by considering his problem on a sphere. The first analysis to include all the essential features of zonal-wind instability was undertaken by Charney (1947). It is probably safe to say that subsequent analyses are elaborations of this model. He incorporated Rossby's:-plane approximation, but assumed a baroclinic hydrostatic atmosphere in which a perfect gas was undergoing adiabatic motion having small Rossby number (Ro) and large Richardson number (Ri). Traveling wave disturbances were assumed to move along parallels, while vertical or meridional variation of certain parameters was eliminated by taking appropriate averages. Solutions were obtained in terms of confluent hypergeometric functions. Among his important results he showed the existence of large-scale unstable waves in a zonal current having vertical shear (or, equivalently, horizontal temperature gradient) exceeding a critical value and described the structure of these waves. He did, however, find his unstable regime to be bounded by a long wave limit due to the. -effect. By neglecting the 3-term Eady (1949) presented a short wave limit to instability due to static stability and obtained the most preferred scale of the unstable waves. A numerical calculation to extend the class of stability problem as posed by Charney and Eady was performed by Green (1960). He examined stability properties of a baroclinic zonal flow bounded above and below by two rigid boundaries, and concluded that there are neither

6 short nor long wave limits to instability. This was supported by Burger (1962) who showed that extrapolation of linear vertical wind shear from the troposphere into the stratosphere leads to unstable disturbances for all but a finite number of isolated wave lengths. Further, Miles (1964), in a general study of baroclinic instability in which he assumed an obliquely traveling small disturbance relative to a zonal wind, found disturbances of all wave lengths in typical flows always unstable for sufficiently small wind speeds. He suggested that small disturbances of typical zonal-wind configurations are unstable for almost all wave lengths. Other related studies were undertaken by Kuo (1952), investigating generalizations of Rayleigh's Theorem, which states that an inflection point in a velocity profile is a necessary condition for instability of an inviscid, homogeneous shear flow, and by Charney and Sterr (1962), who proved Rayleigh-like theorems for jets. Another class of theorems of unstable disturbances referred to as semicircle theorems relates to their rate-of-growth. Howard (1961) constructed a general quadratic integral for the instability of gravity waves in stratified shear flows from which he derived a semi-circle theorem stipulating that complex wave speeds for unstable gravity waves in stratified shear flows must lie within the semi-circle based on the range of the basic zonal wind speed as diameter. Shortly after, Eckart (1963) extended Howard's theorem to adiabatic jets. His extensions were less explicit than the several semi-circle theorems for unstable disturbances contained in a study by Miles (1964) devoted to the general baroclinic instability problem. A year later Miles (1965), considering limiting cases of his 1964 paper, proved a theorem for very long waves in a

7 zonal flow nearly identical to Howard's for gravity waves in a stratified shear flow: complex wave speeds for unstable waves must lie within the upper half of a circle having its center at the minimum wind speed and its radius equal to the range of wind speeds. Study of the effects of large-scale external forces on the existence of quasi-stationary long waves was pioneered by Charney and Eliassen (1949), who showed the importance of large-scale mountains and friction, and Smagorinsky (1953), who studied the ultra-long waves created by large-scale heat sources and sinks and modified by frictiorn Later, long wave studies by Doos (1962), who included diabatic heating; Sa-ltzman (1965), who incorporated mountains and heating; and Derome (1968), who investigated the effect of lower boundary topography, distribution of heat sources and sinks, and friction on the maintenance of the time-averaged standing eddies, augmented these works. Recent studies of transient ultra-long waves and their representation by a barotropic model include those of Deland (1964, 1965) and Eliasen and Machenhauer (1965, 1969). Deland's (1964) study is strictly observational. A Fourier zonal harmonic analysis of the 500 mb height field is made, with quadrature spectra calculated for the amplitudes of the cosine and sine harmonic components. This is amplified in his 1965 study in which he makes a comparison between wave speeds calculated from daily surface-spherical harmonic expansions of 500 mb height in the Northern Hemisphere and theoretical nondivergent Rossby-Haurwitz wave speeds. It is observed that there exists a systematic difference between these wave speeds, i.e., all of the observed waves show greater eastward, or less rapid westward, motion than the Rossby

8 Haurwitz predicted wave speeds. A similar study over the Northern Hemisphere by Eliasen and Machenhauer (1965) uses the spherical harmonic expansions of 500- and 1000-mb height to obtain a representation of the nondivergent part of the large-scale horizontal motion in terms of the spherical-harmonic components of the stream function. Mean values of the velocity of propagation obtained from the stream field for the different harmonic components are then compared with the velocities determined by the Rossby effect modified by weak divergence and found nearly in accordance. In Eliasen and Machenhauer (1969) the large-scale waves are represented by spherical harmonic components of the height field for the whole earth as well as for each of the two hemispheres and no attempt is made to compute stream function components. Velocities of propagation and amplitudes of fluctuations of the very long waves are then, as in their earlier work, compared with the Rossby effect combined with a divergence effect and essential agreement realized. These studies of Deland, Eliasen, and Machenhauer thus conclude that the transient part of the ultralong waves to some extent can be described by a divergent barotropic model. In attempting to explain characteristics of the very long transient waves in the baroclinic mode through the use of a mathematical model, both the analytical and numerical approaches have been applied. Burger (1958) derived through the vehicle of scale analysis a quasi-stationary vorticity equation in the very long wave regime, an expression derivable from the geostrophic relations. This effectively threw doubt on the simplified vorticity equation from which most earlier baroclinic theories began. Phillips (1963) later amplified Burger's conclusions in a review of the general problem of baroclinic instability within

9 the framework of geostrophic motion. Weland-er (1961) found analytically explicit solutions for only neutral traveling waves of very long wave length in a model in which both wind speed and density were linear in pressure. Assuming strong rotational constraint (Rossby number (Ro) < < 1), strong dynamic stability (Richardson number (Ri) > > 1), and weak stratification or static stability (fractional change with height of potential temperature < < 1), he obtained for waves having wave length > 10,000 km, a general wave equation for small, frictionless, adiabatic perturbations in zonal atmospheric flow. His solution, worked out for the beta-plane approximation, consisted of slowly-moving waves, independent of wave number. Another analytical attack on the internal dynamics of ultra-long baroclinic transient waves was launched by Wiin-Nielsen (1961) which might be considered an extension of Welander's (1961) study containing more simplifications but no restriction to neutral waves. Utilizing Burger's (1958) quasistationary vorticity equation he investigated the stability and structure of the waves. In the two-parameter model he showed that the stationarity of the vorticity equation (neglect of 5c/ f) acts as a filter to eliminate Rossby waves from the solution of the linearized equations, while the remaining waves move with a speed independent of wave number. Expanding to the threeparameter model, the smallest resolution allowing for the possibility of unstable solutions, he found only slowly-moving waves possible and the appearance of unstable waves for sufficiently large vertical wind shears. The studies by Miles (1964, 1965) referred to above are both analytical. In the earlier study, for adiabatic, nonviscous motion of a perfect

10 gas with small Ro and large Ri, the equations are reduced to a singular Sturm-Liouville problem and approximations deduced for small and large wave numbers. The latitude (cp) appears as an implicit parameter. Much attention is devoted to the singular Sturm-Liouville problem, but no physical interpretations of unstable waves such as their structure or energy conversion process are offered. Miles (1965) formulated the eigenvalue problem governing baroclinic disturbances relative to zonal flow of waves having length comparable to the circumference of the earth. Several theorems governing the existence of such disturbances were deduced, but calculations of the complex phase velocity were made in only two limiting cases: 3 -+ 0, oo where = Ro Ricot cp. Most recently Hirota (1968), using the quasi-geostrophic equations, investigated numerically the stability of a zonal current and characteristics of wave perturbations by a finite difference approximation of the linearized perturbation equations applied to a multi-layer model. He found three types of baroclinic transient wave solution: (1) a pair of amplifying and decaying waves, (2) neutral waves without steering levels (levels where the basic zonal current and wave speed are equivalent) in the basic current, and (3) neutral waves with a steering level. Thorough discussion is given to the nature of the unstable waves-their vertical structure and energy conversion process-and to the comparison between the characteristics from theory and those of the ultra-long waves observed in the atmosphere.

11 1.3. OUTLINE OF THE STUDY Of the research papers referred to in the previous section this study most nearly resembles the most recent, Hirota's (1968). Like him a numerical solution to thelinearized perturbation equations for a wave perturbation on a zonal current in a multi-layer model is obtained and much effort given to the vertical structure and energy conversion process of the unstable, baroclinic transient waves. Unlike him, however, we assume geostrophy rather than quasigeostrophy and, besides investigating adiabatic flow, obtain complete solutions for the unstable waves modified by the addition of Newtonian heating. Detailed comparison between the two models is made with the intent of determining the effect of Newtonian heating on the wave stability, structure, and energetics. Although the assumption of a basic zonal wind speed linear with pressure breaks no new ground, the consideration of its meridional variation as determined from observation for three specific cases —January, July, and the annual average, 1963 —between 20~N and 85~N forges into unfurrowed terrain. This, incorporated into both the adiabatic and Newtonian heating models in the ultra-long wave regime (wave number = 1, 2, 3) is the major contribution of this study to the state-of-the-art. In Chapter II the assumptions and basic equations to be used for adiabatic flow are presented, manipulated, linearized, and reduced to a single secondorder ordinary differential equation. Numerical solution of the equation requires one exact solution, and a procedure is given to determine it. Derivations of a semi-circle theorem and sufficient condition for stability follow. Modification of the equations to include Newtonian heating and a simple

12 technique to compute its eigenvalues from the adiabatic solutions are then given. The chapter concludes with a description of the numerical method used to compute the adiabatic eigenvalues. In Chapter III the results excluding energetics are presented and discussed. Stability of the waves as a function of latitude is first considered followed by the variation of phase speed with latitude. The chapter concludes with the vertical structure of the most unstable waves in January, July, and the annual average as indicated by the phase angle and amplitude variation with pressure of geopotential, vertical pressure velocity, and temperature. Chapter IV focuses on the energetics of the two models, specifically the pressure variation of two parameters: conversion from zonal to eddy available potential energy, C(A, A), and conversion from eddy available potentialto4 eddy kinetic energy, C(AE, K). In Chapter V, the results of the study are reviewed and some suggestions for future improvements are offered.

CHAPTER II FORMULATION AND METHOD OF SOLUTION 2.1. FORMULATION 2.1.1. Assumptions and Basic Equations We assume that the gravitational field, which includes centrifugal forces due to the earth's basic rotation, is uniform and directed vertically toward -2 the earth's center; thus, we let the acceleration of gravity, g = 9.8 m sec. Consistently, we imply a spherical earth. We neglect mountains and friction, not that these are of no importance, but partly for convenience in solving the equations and partly because we desire to isolate the system's internal dynamical factors. For ultra-long waves we consider wave number n = 1, 2, 3, giving a horizontal scale of motion, L 10 m, the order of magnitude of the radius of the earth, a = (2/Ar) x 10 m. The vertical scale of motion, H, 4 approximates 10 m, while we assume the characteristic horizontal and vertical -1 speeds, from observation, to be V 15 m sec and W < 0.1 V, respectively. For these scale values the equation for the vertical component of velocity can then be approximated by the hydrostatic relation, in which case it is convenient to use pressure instead of height as vertical coordinate. By using pressure as vertical coordinate we eliminate density from the equations of motion and continuity. For Rossby number, Ro = V/fL ~ 1, where the Coriolis parameter f = 2 0. -4 -1 - -1 sin cp - 10 sec for mid-latitudes and. = 7.3 x 10 sec is the earth's 13

14 angular velocity, the horizontal pressure and Coriolis forces in the horizontal equations of motion nearly balance each other so that the geostrophic approximation is well satisfied as shown by Burger's (1958) and Miles' (1965) scale analyses. Further, the type of motion we consider is an example of geostrophic motion of type 2 (see Phillips (1963)) and, as shown by Phillips, use of the P-plane is invalid. In spherical coordinates the basic equations take the form 1 a~ u - -- (2.1) 2 0 a sin cq (. 1 C v- ~~ (2.2) 2 Q a sin p cos cp (2 = - a (2.3) a 1 - u + vcosc (2.4) ap a cos cp a p a-a u a Bk +eC = 0 (2.5) at ra a. cos Tp 2aP,, a.t a C ap where we assume the equation of state is the ideal gas law, p = RT (2.6) Equations (2.1) and (2.2) express the geostrophic relations, (2.3) the hydrostatic equation, (2.4) the continuity equation, expressing conservation of mass, and (2.5) the thermodynamic energy equation, expressing conservation of energy as embodied in the First Law of Thermodynamics, where we first assume adiabatic flow. The independent variables are time (t), latitude (cp), longitude (X), and pressure (p), where latitude and longitude increase northward

15 and eastward, respectively. The dependent variables are the velocity components (u,v,co), where u increases eastward, v northward, and w = dp/dt is the "vertical velocity"; geopotential C = gz, where z is the height above mean sea level of isobaric surface p; and specific volume a. T represents the temperature and R = 287 kj t deg is the gas constant for air. In Eq. (2.5) a, a measure of the static stability, is expressed as o = - n 0 (2.7) where 0 is potential temperature, the temperature a. parcel of dry air would have if brought adiabatically from its initial state to standard pressure p, usually 1000 mb. From Poisson's equation e T( pR/ (2.8) where c is the specific heat capacity at constant pressure for air. Substituting (2.3), (2.6), and (2.8) into (2.7) we obtain 2o = -- (2.9 a -+ - - 9 o 2 c p (2.9) ~p p where c = c - R is the specific heat capacity at constant volume for air. v p 2.1.2. Linearization and Reduction of the Basic Equations Inserting (2.1) and (2.2) into (2.4) we find __ 1 _a 2 (2.10) 2 2Q a sin cp In solving (2.1), (2.2), (2.5), (2.9), and (2.10) we desire to linearize the

16 equations by perturbation theory in order to simplify the mathematics such that, with appropriate additional assumptions, we can reduce the set of equations to a. single second-order ordinary differential equation in one unknown. Consider a steady, axially symmetric basic state having only zonal wind with perturbations superimposed upon it, u = U(p,p) + u'(t,,cp,p) (2.11a) v = v'(t,,c,,p) (2.11b) X = o'(t,\,q,p) (2.11c) ~ = $ (yp) + D'(t,,cp,p) (2.l1d) o = a (Cp) + a'(t,,qp,p) (2.11e) where the primed terms are perturbations. Because we work with spherical coordinates, it is convenient to define the basic state angular velocity, A(cpp) -.= (2 (2.12) a cos cp In determining the basic state solution we assume U(q,p) is given and express 0 (gp) and a (q,p) in terms of it. From (2.11) into (2.1) and (2.9), we find z Z z 2 a_ = - 2 X a sin p U -2.0 a sin qp cos p A (2.13) 20 c ac za + (2.14) z 62 c p ap ap p Differentiating (2.13) with respect to pressure gives Za~z~) = -2 a sin p cos ( ) <^p \. 6p 6p

17 while by integrating (2.13) with respect to latitude we get Z (Pp) z (P) - 2 0 a f U(cp',p) sin cp' d p' (2.16) i Z L Cp L where cp is the lower latitude limit, either 20~N or 25~N depending on the data L used, cp' the integration variable, and zL = at cpL. Inserting (2.16) into (2.14) completes the basic state solution, zo = z 2 2 a a sincp' + -v 1 d c' (2.17) z zL 2L 2 cp 1 where GzL = at qPL To obtain perturbation equations we substitute (2.11) and (2.12) into (2.5) and (2.10), giving _ ( -) + A )+a c = 0 (2.18) at \ ap/ c^ \6ap a acp amp z ~c' 1'D ap — 2 2 — ax (2.19) 2 2 a. sin cp After re-expressing the third term of (2.18) using (2.2) and (2.15), (2.18) becomes, + + a = 0 (2.20) at a~ ~P ~p a z Using the method of the separation of variables we introduce perturbations of the form ( ) = ()e in(X-ct) (

18 where ( ) is assumed to be a function of cp and p, n is the zonal wave number, and c is the phase velocity of the perturbations and is in general complex, c = cr + i c. We thus assume that the perturbations are wavelike disturbances, periodic in the east-west direction and time, and that the flow can be represented by a single component of a Fourier series in X. For ultra-long waves we consider only n = 1, 2, and 3. For c. > 0 the rate of growth is n ci, the n c* t motion being unstable through the growth factor, e i, multiplying all the dependent variables. For c. = 0, we term the flow neutral and for c. < 0, stable. Equation (2.21) into (2.19) and (2.20) gives A _ in (2.22) ap 2 2 2 0 a sin cp inL(A-c) ^ +- C A = 0 (2.23) which, after substituting (2.22) into (2.23) reduces to the single equation, 2A A C (A- c) a- - W 0 (2.24) (A-c) 2 + 2 U 0 ap2 ~ p up 2 R a sin2 c In Eq. (2.24) cp enters only as an implicit parameter, not as an independent variable. In physical terms this means the behavior of the ultra-long waves is determined essentially by the vertical and not the horizontal structure of the atmosphere. Because we solve (2.24) or related forms for given cp and U(cp), we may consider l at the specified latitude to be solely a function of pressure, and thus could replace the partial derivatives by total derivatives. Equation (2.24) and related forms would thus become ordinary differential equations in

19 this localized sense. We shall, however, continue to express (2.24) and related forms as partial differential equations to emphasize the in-the-large A dependence of c on latitude and pressure. We specify the boundary conditions for cw to be a) = 0 at p = 0 and p = p = 1000 mb. The lower condition is obtained by expanding w in the z-coordinate system, removing the horizontal pressure advection through the geostrophic approximation, and assuming a negligible pressure tendency to arrive at ap a w W - 6z as shown in deriving Eq. (4.5). Assumption of level terrain at pressure p implies w(p=p ) = 0 from which we get the desired boundary condition. 0 In seeking a numerical solution we enhance the accuracy of our calculations by eliminating in (2.24) the first derivative, recasting the eigenvalue problem in canonical form. This is achieved by letting X = M W for M = A-c and substituting in (2.24) to get, after dividing by M' 2a I I 2M 3 1 aM 2 z 1 _ +21"~~2 + z~2~ M Iw = 0 (2.25) ap L 2 M 4M2 2 sn (p ~ 6p M "p 2 P a sin (p It is convenient to nondimensionalize pressure and phase velocity by p = c c (2.26) 0 0 where IA = sA, implying that s = +1 for A > 0 and s = -1 for A < 0. We assume a linear variation of U with pressure at each latitude, We assume a linear variation of U with pressure at each latitude, U(q),p) = U (q))(i-p) or A(q,p,) = A(q)(l-p,) (2.27) " 0 o o

20 0 I I / i/'lI /^^ I /I Given cp // I I 0 Ao UO U,A To illustrate the validity of this, we compare below (see Welander (1961)) normalized profiles of U*, where U. = 1 - p. is the assumed linear profile, and the observed is computed from mean winter conditions at 50~N. P* O y~bS Obs U* k From observations a is proportional to p where k varies between -1 and -2. z For mathematical expediency we choose the former, =a - 0 (2.28) zL p* We use throughout this study a = 1 in the MTS system, this being based upon results obtained by Peixoto (1960) in which the difference between computed values of a at 45~N from average data for the winter and summer of 1950 was very small compared with the average of the two seasons. It being always

21 desirable to have analytical solutions wherever possible to check our numerical results, we consider in Appendix A.1 the barotropic limiting case, A = A (qy), a =a constant, while in Appendix A.2 we investigate a more realistic baroz 0 tropic limiting case more similar to our model, A = A(cp), a = /p The o z 0/P The results from these, pertaining to perturbation phase speed, will be compared with our model in Section 3.2. It has earlier been stated that in the limiting case of Ro = V/fL < 1, i.e., for the horizontal scale of wave motion-here assumed in the west-east direction-comparable to the radius of the earth, scale analyses of Burger (1958) and Miles (1965) showed the geostrophic approximation to be well satisfied. That is the assumption used throughout this study. Assuming the earth's radius to be approximately 6.6 x 103 kmn, Figure 1 illustrates that, for 20~ < cp 85~, the wave length of wave numbers 1, 2, and 3 is greater than or equal to the earth's radius over much of the latitude range considered, but that, especially for wave number 3 in the higher latitudes, the radius of the earth does exceed the wave length and the geostrophic approximation thus becomes strained. This is one justification for not including wave number 4 in this study. Another approach is to obtain analytical solutions for simple quasi-geostrophic and geostrophic models and to compare the phase velocities as a function of wave length. This is demonstrated in Appendix A.3 for the basic state of no motion and it is shown that at 450N the percentage errors in assuming geostrophy for wave numbers 1, 2, and 3 are 1.6o, 4.6% and 10o, respectively. Thus, for the limiting case of no motion in the basic state, wave number 4 must be excluded in the geostrophic model at 450N if errors less than 10t are to be

22 obtained, further justification to restrict the wave numbers in this study to the first three. Now substituting (2.17), (2.26), (2.27), and (2.28) into (2.25) gives us, where M = Ao(l-p -scC), __ 1 o_______ _ + ~ cos 0 ~ U sin p' d 2 *-(1-sc 2 a sin cp U sin p L p. * 2 A Aasin 2 o o + -1 W = 0 (2.29) 4 [p (1-sc 2 where c = c /c P The boundary conditions are W = 0 at p = 0, 1. We avoid regularity difficulties at the equator by confining the range of latitudes to 20-85 N. Further, we examine only cases having Uo # 0, precluding a similar type of singularity. Finally, investigation of (2.29) reveals the existence of singularities at p* = 0 and p* = 1 - sc.. The former is avoided by virtue of the fact that it coincides with the upper boundary where (2.29) is not applied. The latter will only occur when c* is real. For all the cases considered only three yielded real c*; Period Latitude Uo c* Annual 25~ 14 m/sec (-.3521928,0) Winter 20~ 25.60 (-.1873492,0) Summer 300 2.16 (-.7394156,0) In each s = +1. The eigenvalue c* < 0 as above implies that singularity can only occur for p, > 1. Thus, all singularities are avoided in the numerical

23 solution of (2.29). For given qp and U (q'), qL < q' < p, we solve for c (qp) using the "shooting method" (see 2.2, Method of Solution). This is a numerical search procedure in the complex plane which requires an exact solution c *e for foe some U at latitude cp to start it and provide a test condition for the search oe e procedure. We now desire to find this exact solution. 2.1.3. Determination of an Exact Solution Inserting (2.17), (2.26), (2.27), and (2.28) into (2.24) we obtain 2 [p*-(1-sc*)] - - ap. A2 a s in2 ~ 6p - 2 ^ * 2 P A a sin2 + _" S 2_ J0L U sin p' d p' = 0 (2.30) U sin cp L Transforming to a new independent variable, 5 = p0 /l-sc*, (2.30) becomes 2 2A rA a p ^,^aw aw ^ _______ oo.... Gau'2 QA a sin cp A + "E- a2+ [(O c+ U sin+)p + y] + c; = 0 (2.321) Gauss's differential equation, known also as the hypergeometric equation, a special case of the generalized hypergeometric equation. The standard form of Gauss's equation is ^5-1) ^ + [(a-l +1) a + +3 0 (2.32) 2$^ Equating (2.31) and (2.32),

24 y = 0, + + 1 = -1 (2.33) 2 3 = -O_ + cos c u sin d p' (2.4) 2 Q A a sin cp U sin c L 0o For q = 0e = qL' (2.34) becomes a - = -q (2.35a) where 2 o P oq (2.35b) - 2 2 2 Q A a sin qP oe e Combining (2.33) and (2.35a), = -1 + r, = -1 -r (2.36) where r = JT+ (2.37) From Kamke (1961), page 467, for y = 0 and a and? as determined in (2.36), the solution to (2.32) is A A + C A (2.38) C 1 2 o2 where A w = 2F (r,-r;2;) (2.39a) X1 2 F1 o (1-r ) in I F (r,-r;2;S) + power series in -2 1 starting with "" (2.1 9b)

25 Gauss's function is O0 (a) (b) x 2F (a,b;d;x) (d) n (2.40) n=O (dn n where (a) = a(a+l) (a+2)... (a+n-l) and a, b, d, and x can be real or complex, n but d cannot be a negative integer. A From (2.39a) and (2.39b), for the upper boundary condition at g =, = 0 and = 1. Since a = 0 here, C must be 0. To satisfy the lower boundary condition we seek e = e such that 0 2F1(r,-r;2;) = 0 (2.41) where c. is determined from e = l/l-sc. Thus, -1 c = (2.42) o From Kamke (1961), for integer n, 1-c( -c-a n c+n-1 a+nF (n+a,-n;c;) = ( )- nl(1)a+ (2.43) 21 c(c+l)... (c+n-1) dn Combining (2.41) and (2.43), n = r, a = O, c = 2 (2.44) which, with (2.37), gives n = or q = n2 (2.45) Replacing q from (2.35b),

26 2 ap 00 2 = n -1 (2.46) 2 2 2 Q A a sin cp oe e which, using (2.12), enables us to calculate U oe 2 o p cos (p U - ~0 e (2.47) 2 fi a sin2 q (n -1) Given ce =,L' this varies only with arbitrary integer n > 1. Combining (2.39a), (2.43), and (2.44) we now get A C (1-) ) d Fn+l (1_ )n-2] W0 C1 (n+l)(:) j-n (2.48) from which we can easily determine c. Let us consider several special cases. (1) n= 2 Equation (2.48) simplifies to Ax <(1-). At the lower boundary condiA tion, o = 0 obtains from e = 0, 1 which substituted into (2.42) and discarding 0 the meaningless solution at e = 0 gives c = 0. 0 *e (2) n =3 Equation (2.48) simplifies to oC S(1-2) (1- 2 5). To satisfy the lower boundary condition, e = 0, 1, 2/5 which gives c = 0, -1.5. 0 *e (3) n=4 A 0~2 2 Equation (2.48) simplifies to wAc> (1-) (7 2 - C6 + 1). The lower boundary condition is satisfied when o = 01, 3/7 + 2/7 which gives c* 0 G0, -2 + 2.

27 2.1.4. Derivation of a Semi-Circle Theorem For M = A-c and p = p p*, (2.24) becomes 2A A w c)M ) A M a.+ Xs 6 X = o (2.49) where 2 P z o s = (2.50) 2 2 2 0 a sin cp and 62 + N 6 a = a zL 2 -2 a ( sin cp' ( 2 d (2.17) z zL CP 2 P L 6p Condition (2.27) gives au Uo(") a2U as Po ap2 while condition (2.28) is azL = a /p. As stated earlier a = 1 for all calculations while U (cp) > 0 in all cases except cp = 20~N, 25~N, and 85~N for July, -1 -1 where IU (qc)l < 10 m sec, while U (q) > 10 m sec for cp = 35~N - 70~N in July. In proving a > 0 for all qp and p, we thus need only consider in detail z the July case at 250N, as qL = 20~N for July. For U (20~N) = -8.06 m/sec, U (250N) = -6.08 m/sec, and a = 1 in MTS (meter-ton-second) units, we find, on evaluation of (2.17), the amplitude of the first term several orders of magnitude greater than the second, confirming the assertion that, for every case considered, a > O, implying S > 0 from (2.50). Combining the first two terms of (249) we get terms of (2.49) we get

28 2 a6 (-Ma) + S At = O (2.51) In assuming there are unstable solutions, i.e., c. > 0, we express c c + i c. (2.52a) r 1 A A A ( = 1+ i i (2.52b) r i A A A ) = ) - i L). (2.52c) r 1 M = (A-c ) - i c. (2.52d) r 1 M = (A-c ) + i c. (2.52e) r 1 where ( ) represents the complex conjugate of ( ) and ( ) ~ (~) = ( )I 2 A Dividing (2.51) by M, multiplying by w and integrating over p* from 0 to 1, A we obtain, after applying the boundary conditions for w, 1 C- 2 11 2 A - dp M S|cu dp + 0 (2.55) which can be alternately expressed as A^ 2 - M p o M SII4 dp* = 0 (2-54) MI 2 +f MI where -2 [(-_2 2 M = [(A-c ) - ci] + i 2 c (A-c ) (2.-5) r 1 1 r Let 2 ^Q = ~ >>0 (2.56a) IMI2

29 and,A2 R = S 2->0 (2.56b) |M|4 where we apply the fact that S > 0 in (2.56b). Dividing (2.54) into real and imaginary parts, after inserting (2.52e), (2.55), (2.56a), and (2.56b), we obtain the set, 1 1 2 2 - f (A-c )Q dp* + f [(A-c ) - c_]R dp* = 0 (2.57) 0 r 0 r * i [- / Q dp + f 2(A-c )R dp ] = 0 (2.58) For c. > 0 and Q > 0, (2.-8) now becomes JO (A-c)Rdp = f Q dp >0 (209) o r 2 * which implies, for R > 0, c < A (2.60) r max Further, A Rdp = c R dp + Qdp (2.61) o r o 2 o From (2.57) we get 1l A2 21 221 1 J A R dp, = -c R dp* + c R dp + 2c J A R dp o r 0 i r 0 + Q dp* -c Q Q dp (2.62) 0w r o which, applying (2.61), simplifies to

30 R dp = (cr+) f Rdp + 1 A Q dp* (2.63) d r i o o Denoting A = b and A = a and recalling that R > 0, we have for A $ max min max min' f (A-a)(A-b)R dp < 0 (2.64) o0 which, expanded, gives 1 A2 R dp - (a+b) f A R dp + ab J R dp* < 0 (2.65) Substituting (2.61) and (2.63) into (2.65), we find 2 2 1 1 [c + c - (a+b) c + ab] J R dp* + 1 A Q dp. r o 0 a+b 1 2 — - Q, d < 0 (2.66) 2 o By suitable manipulation the following identity can be established, 2 2 2 2 2 c + c. - (a+b) c + ab = (c -a) + c. - (b-a) r i r r i + (b-c )(b-a) (2.67) r which, inserted in (2.66), gives 2 2 2 1 1 [(c -a) + c. - (b-a) ] R dp + (b-c )(b-a) f R dp r 0 r 0o 1 (A- a+b)Q dp < 0 (2.68) o 2 * From the definition of a as A min

31 1 af b a+b b-a 1 o ( - 2)Q dp > f (a - )Q dp - - dp (2.69) o0 ^ d2..2 Now, substituting (2.61) into the right side of (2.69), recalling that b = A, we obtain max O (A - -)Q dp* > - (b-a) fo (A-cr)R dp* > (b-a)(b-c ) r Rdp which transforms to f1 (A - a)Q dp + (b-a)(b-c) If R dp > 0 (2.70) It now follows from (2.68) that 2 2 2 (cr-a) + c < (b-a) (2.71) Thus, for ci > 0 the point (c,ci) must fall inside a semi-circle with center (a,o) and radius, b-a, as illustrated below. Ci (a,b-a) (2a-b,o) (a,o) (b,o) Cr For our linear profile A = A (1-p.) (2.27) we investigate two cases.

(1) A > 0 o Here a = 0, b = A and we find 0 -A < c < +A, c. < A o r o 1 o which gives, for c = c/IA I -1 < <, cci < 1 We conclude that, in westerlies, dimensionless complex wave speeds for unstable waves must lie within the upper half of the unit circle centered at the origin. (2) A < 0 Here a = A b = O and we find 0) o 2 A < c < 0 c. < -A O r 1 o which gives, for c* = c/jA -2< < c*r < c <1 Thus, in easterlies, dimensionles complex wave speeds for unstable waves must lie within the upper half of the unit circle centered at -1 on the real axis. While this semi-circle theorem specifies that unstable dimensionless wave speeds must lie within a. particular region of the complex c*-plane, it makes no statement about the uniqueness of the solution of (2.29) for given cp and U subject to the specified boundary conditions. Discussion of the uniqueness of our unstable eigenvalue solutions is shown in Appendix A.4.

33 2.1.5. Derivation of a Sufficient Condition for Stability We define stability by the condition c. = 0. For M = A-c, (2.22) and (2.23) may be re-expressed as 1 a in A p aP 2 asi (2.22)' Po *2 2 a sin cp in 3 MA^ A in ( ~r M a~ —^ ]+ o = 0 (2.23) Po ~ap* p z A A Whereas we earlier eliminated $ from this set, we now choose to eliminate o, giving ~p. 2 = 2 I ~p 22 2 ap p J 2 Q a sin cp z o which is easily transformed to M M A aP a p2 aI —+~ O (2.72) zo 2 2 2 Q a sin cp subject to the boundary conditions, using (2.23)', ~M a = 0 for p =0 and p* =1 (273) Substituting (2.50) in (2.72), we now find _~ M1 1- ~M) + 0 P* S 6p* S p* which expands to

34 M a+~ 41 p aM 1 a 1 aM a a( a;+ 0 o 6p* s -p* s aap ap ap S- s p ) Cancelling and collecting terms, we arrive at M. ).\+ [ l ( e o which may be put in the form a (i a + = 0 (2.74) where aN 1 ap (2-75) After multiplying by $ and integrating over p* (2.74) becomes 1 a (D 1 1 2 d I f ap s dp S ap dp* A 2 + f N M L dp o (2.76) IMi2 Integrating the first term and applying (2.73), we get 1 1 aM jAj2 1 1 aM\ A 2 S1 1 p 1 So M a0 1 2 d + 1 N M 2 wdNM he (277) where

35 ( )1 = ( ) at p = 1 ( ) = ( ) at p = 0 Combining (2.50), (2.17), and (2.28), we arrive at the boundary condition S + oas p* + 0 (2.78) o Taking the imaginary part of (2.77) and applying (2.78), (2.52d) and (2.52e), we obtain i S.c )1 IP + l MI2 p 0 (2.79) 1~ S a 1iM 2 o 2 Thus, if N has the same sign as (Ma/ap ) in the interval 0 < p* < 1, (2.79) can only be satisfied if c. = 0. This gives us the sufficient condition for stability that N and (/ap*)l, have the same sign for 0 < p* < 1. For our linear profile, a/Vap = -A and is negative for westerlies, positive for 0o easterlies. From (2-75) we find N =1 +A - (2.80) o 6ap S which, after substituting (2.50), (2117), (2.28), and (2.27) gives us 2 2 (2 0 a sin cp)A, N = 1+ 2 (2.81) p2 + 2 a U (') sin d 1 caos 2 CP o where S > 0 for all cases implies that the denominator is positive for all cases. Thus, N is positive for westerlies and either positive or negative for easterlies, depending on the relative magnitudes of the two terms in (2.81).

36 To conclude, for the linear profile the sufficiency condition for stability can never be satisfied for westerlies but may be for easterlies. 2.1.6. The Newtonian Heating Model Here we retain all of the assumptions made previously except the adiabatic assumption. Thus, the thermodynamic energy equation, (2.5), now becomes a aD u a c au v aft ^ Rl, at,p, a cos p \app a, p\p/ c p where Q is the diabatic heating per unit time and unit mass. We consider the diabatic heating to be Newtonian in form so that Q = c q(T -T) (2.82) p E where q is the heating or cooling coefficient and T (cp,p) is the equilibrium E temperature toward which the diabatic processes are driving the atmosphere. This is the hypothetical temperature which would be established in the absence of large-scale motion but in the presence of radiation and small-scale convec-6 -1 tion. A value of 0.4 x 10 sec for q was obtained by Wiin-Nielsen, Vernekar, and Yang (1967) using the calculated value of the generation of eddy available potential energy from observations. This constant value for q will be used throughout this study. The Newtonian heating assumption, Eq. (2.82), says that the atmosphere is heated in regions where the temperature is below the equilibrium temperature, while it is cooled where the temperature exceeds the equilibrium temperature. It certainly represents an oversimplification of the heating in the atmosphere,

37 especially so because of the role of condensation in large-scale flow patterns as studied by Fjortoft (1959). However, Eq. (2.82) does agree with the average behavior of the atmosphere, and it is on this basis that we justify its use. Expressing temperature as a steady, axially symmetric basic state upon which is superimposed a perturbation, we have T = T (c,p) + T'(t,\,p,p) (2.11f) Inserting (2.11a)-(a.llf) and (2.82) into (2.5)' we obtain at8 Z + + a o o a ( + at' + V ( + at p p a cos Tp C ap apG al ap ap p + (az+a') = -(TE -T'-) (2.83) which gives the basic state expression (T-T) = 0 (2.84) P -E z In seeking a revised perturbation energy equation we first express T' as a function of 6a'/ap through the ideal gas law, Eq. (2.6), and hydrostatic equation, Eq. (2.3), - TTt =-=p (2.85) R R R p and then substitute (2.12) and (2.85) into (2.83) to get _'\+ (oA A', (21 _t _\p) apem I - A- (2. 18)

Introducing perturbations of the form (2.21) into (2.18)' we arrive at A^A in (A-c) ap- ap + a = -q P P~ z 6P Eliminating $ by (2.22), it follows that 0-++ Wu= 0 L iq ap p — Ap 2.2 i 6p2 6p 6p 2 X a sin p which is easily transformed, after applying (2.17), (2.2) (2.2(.7), and (2.28) to 2 2 P, - (l_-scN) + q i 0 0 L o p2 2 a A a sin q) O + e CfS J U sin p' d q = 0 (2.30)' U sin qp L 0 where cN is the nondimensional phase velocity for the Newtonian heating model. In seeking to determine c it is very useful to compare (2.30)' with its adiabatic equivalent, 2 2A pA 0P [p. f (-sc o 1 o0 + COS qp * 2 _Sa)p p2 2 ~2 pa p2 p* 2 A a sin qp U sin cp 7~ U sin p' d qcp = 0 (2.30) 0L where c is the nondimensional phase velocity for the adiabatic model. We note the equivalence of the two equations in all terms but the coefficient of 2 2 a /ap* which, for the Newtonian heating model, may be expressed as (p* - [I - s(cN +(q/sn A )i)]} as compared with [p - (1-sc )] for the adiabatic model. It is therefore clear tha.t, for given cp, n, and A, a solution c*a

39 will suffice to calculate c*N from the equality, c + q i = *N sn A *a 0 or, applying (2.12), c c* q a cos cp (2.86) *N - a- sn U o This precludes the necessity to apply the "shooting method" for the Newtonian heating model once the adiabatic phase velocities have been obtained. We note that, whereas c is independent of zonal wave number n, cN is dependent upon it in the sense that increasing n tends to cause perturbations in both westerlies and easterlies to become more unstable, as s U > 0 in either case. 0 In contrasting the two models we shall show that complex values of phase velocity necessarily occur in complex conjugate pairs for the adiabatic model but not for the Newtonian heating model. As derived earlier for adiabatic flow, 2A 2 2 c Ac a sin q U sin q L * r *i'r i so that (2.30) can be expressed as two separate equations, one for the real so that (2.50) can be expressed as two separate equations, one for the real

40 components, 2A 2A A 00 a6 W. aA (p*l+scr) 2 sc F = 0 (2.87) ap ap a * and the other for the imaginaries, 2A 2A A a l. a AO. (i+ ) 1 r 1 A (p -+sc-r) a2+ SC- - F W. 0 (2.88) (P*-+SC r) 2 *i 2 ap p 1 If we now consider the complex conjugate eigenvalue and eigenfunction, c = A A A c* - i c*i, = c - i i. and substitute into (2.30), replacing c* by c and C-#rr 1r A A X by c, and resolve into real and imaginary components, we arrive at (2.87) explicitly and (2.88) with each term multiplied by -1. Thus, in the adiabatic A model for given cp and U, each eigenfunction c with its eigenvalue c* has a corresponding complex conjugate eigenfunction with its eigenvalue c corresponding complex conjugate eigenfunction W with its eigenvalue c., illustrating that to each damped stable wave there corresponds an amplified unstable wave, and vice versa. For Newtonian heating Eq. (2.86) gives us, for each adiabatic pair of complex conjugate eigenvalues, the Newtonian counterpart in which only the imaginary components are affected. We thus obtain two eigenvalues which necessarily are not complex conjugates, one of which must represent damped stable wave motion. The other eigenvalue may represent either damped stable or amplified unstable wave motion, depending upon the degree of instability of its unstable adiabatic counterpart and the magnitude of the stabilizing effect of Newtonian heating.

41 2.2. METHOD OF SOLUTION For given pq and U (qp'), qpL < q < cp, we desire to solve o L _~2 CY/ p2 2W _____oo C_ cos K* [p.-(l-sc 2)] 2 2;3p.2 a a sin q2 U sin 0 4 4 U sin cp d w = o (2.29) L o / [p.-(i-sc.)] under the boundary conditions, W = O at p = 0, 1, for c (cp) by a numerical search procedure in the complex plane. This search procedure is referred to as both the "shooting method" and the "Biblical method," the latter deriving from the Scriptural text, Matthew 7:7, "Seek and ye shall find," for the process is based upon an iterative algorithm whereby a crude eigenvalue is continuously changed until the boundary conditions are "satisfied." After the finite difference equivalent of (2.29) is obtained and the atmosphere divided into N levels-N = 40 was chosen for all eigenvalue calculations-W = 0 is specified at the upper boundary and an arbitrary value at the first level below the upper boundary. Successive values of W are then calculated from the finite difference equation at lower levels until W at p* = 1 is determined. This absolute value is then normalized by dividing by the maximum absolute value of W. We notate this normalized W at the lower boundary by G. In the absence of computer round-off and truncation errors, G = 0 for each eigenvalue. This not being the situation, however, a small number, 6, must be calculated to serve as an upper limit for the lower boundary condition in eigenvalue calculations. It is in the determination of 6 that an exact solution, c* for some U at latitude cq, to (2.29) is required. To eliminate the integral in (2.29) oe e

42 Pe is chosen as cpL In calculating U from (2.47), oe 2 o p~ cos qe U e (2.47) oe 2 2 2 Q a sin 2e (n-1) for positive integer n > 1, there is a. choice to be made in our selection of n. As shown in Section 2.1.3., c*e is also a. function of n whose dependence is summarized below; c n *e 2 0 3 0, -1.5 4 0, -2 +2, -2 -2 In determining the optimum selection of n and c we tabulate 6 and, for = 25~N and U = 14 m/sec as the annual average of zonal velocity at cpL' c* for each combination of n and c *e e 6 c(cp = 25, U = 14 m/sec) 2 0.523 x 10-2 -13499997 3 0.157 x 10-1 -.3499997 3 -1.5.518 x 10-4 -.1524882 4 0.327 x 10-1 -.3499997 4 -2 + 2.560 x 10-l -.3521928 4 -2 - 42.882 x 10-4 -.3524940 It is clear that the smaller the value of 6, the more accurate the computed eigenvalue will be. At the same time errors in the value of U preclude more o precise calculations of c than to three significant figures. As computation time will increase for a more restrictive convergence criterion, it is desired

43 to select the largest 5 giving accuracy of c* to the third significant figure. From the table we see that this condition is satisfied for n= 4, c e =-2 +2, 6 =.560 x 10. All eigenvalue calculations were based upon this convergence criterion. The search procedure follows for latitude cL and zonal velocity U: (1) For grid increments Ac =.05 and Ac*i =.01 integrate the finite difference equation and evaluate G at each of the eight points surrounding c. (2) Find the point for which G is a minimum. Call this P1 and the value of G, G1. (3) Test G1 < (a) If true, P1 is the desired eigenvalue and we are through. (b) If false, repeat step #1, replacing c* by P1 (4) Repeat step #2, replacing ( ) by ( )2 (5) Test G2 < 5 (a) If true, P2 is the desired eigenvalue and we are through. (b) If false, test G2 < G1. [1] If true, repeat step #1, replacing c*e by P. [2] If false, halve Ac* and Ac. and repeat step #1, replacing c*e by P Continue in this manner retaining the most recent values of Ac* and Ac*i about P +if G+ < G, or halving the most recent values of Ac and Ac.i about P if G > G until G < 6. n+l n Once an eigenvalue has been determined, it is used as the starting point in computing the next, for new p and U, while refinements in Acr and Ac*i, based upon the relative positions of the two preceding eigenvalues, can be

44 applied to quicken convergence. In the event that convergence does not occur, implying that the real eigenvalue is in another portion of the c -plane, the values of Ac and Ac*i used at the beginning of the unsuccessful computation are increased and the search procedure repeated. If convergence still eludes us, we increase Ac and Ac i once more and repeat. In every situation this technique succeeded in zeroing in on the elusive eigenvalue.

CHAPTER III PRESENTATION AND DISCUSSION OF RESULTS 3. 1. THE STABILITY ANALYSIS It is useful, in order to more readily ascertain the wave length of the perturbations considered, to plot wave length ( ) vs. latitude (cp) as shown in Figure 1 for the range of latitude, 20-85~N, of the data and wave numbers, n = 1,2, and 3, for the ultra-long waves. Clearly, A = 2 J a cos gp / n where 2 i a cos cp is the length of the latitude circle at (p. Having assumed a linear variation of basic zonal velocity with pressure, U(, p) = U((p) (l-p*) (2.27) where U (cp) = basic zonal velocity at latitude cp and p* = O, we obtain the o0 meridional distribution of U by doubling values of the zonal average of U at o 500 mb from data taken during 1963. We consider three separate cases- January, July, and the annual average. Figure 2 displays U vs. p for these cases. It should be noted that, in a large range of latitude (45~ < 9 > 75~) U for the annual average is not sandwiched by corresponding U in January and 0 0 July. If the January and July profiles represented six-month averages for winter and summer, this would certainly be the case, but the fact that profiles of ten other months are averaged with those of January and July to yield the annual average evidently causes the anamoly in which, at certain latitudes, 45

46 40 38 36 34 32 30 28 26 - 24 22 E 0....... 18 - 10 -3.... 4 B....~~~''.. x 2- ""~"'. 0 20 25 30 35 40 45 50 55 60 65 7075 8085 LATITUDE Figure 1. Wave length as a function of latitude for wave numbers 1, 2, and 3.

47 6050 40: T / \ E \/ January E / \ 30 Annual 201/ \0 -10 /July LATITUDE Figure 2. Zonal wind at p = 0 as a function of latitude as calculated for a linear profile by doubling 500 mb data of January, July, and the annual average, 1963.

48 U for the annual average exceeds corresponding values in January and July and, 0 at other latitudes, is exceeded by both. Thus, we should not expect the degree of instability for the annual average to fall between values for January and July at every latitude. Future graphs confirm this. The greater resemblance between U (cp) for July and the annual average than either with January hints o that instability results should in general portray the same agreement, a conjecture which is also confirmed in future graphs when comparing the three cases. In expressing degree of instability we choose to use two parameters based upon the assumed exponential form of the perturbations, ( ), ( ") ein(-ct) (2.21) (A (2.21) They are c.i the imaginary component of the generally complex phase velocity of the perturbations, and T, the e-folding time or time it takes for the amplitude of the perturbations to increase by a factor of e (22.7). Some authors prefer the parameter T = (in 2) T, the time for the perturbation amplitude to double, but we shall utilize the more common T. From (2.21) it follows that T = l/nc.. Here, c. = 1A c*i is expressed in radians/ e i i o sec. To convert into m/sec we multiply by a cos c, noting that IU I = a cos p 0 IA. From (2.29), inclusion of / U sin y' d p' indicates that, for given ~L latitude, the computed eigenvalue will be dependent upon, not just U at that latitude, but on wind profiles to its south. Thus, the eigenvalue calculations are not independent of each other. One sees further from (2.29) that c and thus ci, in the adiabatic case is independent of wave number n, while 1

49 (2.86) relating c* in the Newtonian and adiabatic models, expresses dependence of the Newtonian c* upon n. If T is chosen as the measure of instability, it is obvious from its expression above that, whether the model be adiabatic or Newtonian, there will be dependence upon n. Figure 3 compares instability as measured by c. in m/sec as a function of latitude for adiabatic flow in the three cases-January, July, and annual average. Its equivalent in terms of T for n = 1 is shown in Figure 4, where the e vertical scale is inverted in order that instability increases as the upper part of the graph is approached, as in Figure 3. In specifying the latitude of maximum instability we are faced with a selection of basing our criterion on the maximum positive value of ci or on the minimum value of T. Comparing i e Figures 3 and 4 one readily sees that the criteria give, in general, different latitudes. In determining which is more valid, we first express T in the following way; ^1 _ 1 _ a cos cp (5 1) e n ic nA olci nJUJ Ic (30 *i o *i This indicates that, regardless of ci., T +e- as the North Pole is approached which is physically untenable. Thus, maximum instability will be determined by Ci. Examining Figure 3 we note greater instability occurring in January than in July or for the annual average over most latitudes, with a maximum of c. = 5.76 m/sec at 25~N. In comparison, the July curve peaks at 40~N with c. = 2.46 m/sec while the curve for the annual average peaks at 35~N with

50 6 I 5-1 I \ 5 x I \ I I (U I I!'' I -, - I/ / I' J ". annual ^ I ^^- \ *' I f ~I I: ". 20 25 30,35 40 45 50 55 60 65 70 75 80 85 LATITUDE Figure 3. Imaginary component of phase velocity as a function of latitude for adiabatic model —January, July, and annual average, 1963.

2 4 6 8 I0 12 X-X, it 14- January 16- 18 _ o Annual"' 20 22 24-. 26 36 _. 4 I \/:I \ 32 34:July 36' 38- a 40 20 2530 35 4045 50 55 60 65 70 75 8085 LATITUDE Figure 4. e-Folding time as a function of latitude for adiabatic model, n = l-January, July, and annual average, 1963.

52 c. = 3.33 m/sec. This is not surprising in view of the fact that, from Figure 2, U (January) surpasses U (July) and U (annual) for all cp except 450 < q < 60~ and that vertical wind shear tends to make waves unstable, a result found or stated in many works on baroclinic instability, in particular Derome and Wiin-Nielsen's (1966) recent study. In each case there is an approximate monotonic structure on both sides of the peak. It is of interest, in comparing Figures 2 and 3, to note that c. is a maximum, not at the latitude where U is a maximum, but 5~ to the south. This is undoubtedly due to the effect of all profiles to the south of cp, although the particular reason why is not clear. Figure 5, which plots c vs. U, shows more vividly the influi o ence of vertical shear on instability. Here, the distribution of U with cp' o as displayed in Figure 2 is applied for all cp'< cp, while U at cp is varied 0 arbitrarily between 0 and 50 m/sec. Eigenvalues are then computed for January at 445~N, July at 30~N and 60~N, and the annual average at 75~N. Clearly, increasing U at cp brings about greater instability. o Returning to Figure 3 we observe absence of instability at only low latitudes, in particular 20~N for January, 30~N for July, and 25~N for the annual average. This corresponds to infinite T and is not represented in Figure 4. e From Figure 2 it is apparent that the only easterly zonal wind profiles occur in July at cp = 20~N, 250N, and 85~N where, from Figure 3, instability is slight. This supports the assertion by Kuo (1952) that "a negative vertical shear is relatively more stable than a positive shear." These results also indicate that the amplitude of the very long waves should be a maximum in

53 6.0 5.7 5.4 5.1. 4.8 July, =4 300/ 4.5 / 4.2 / / 3.9 ~ 5 10152025303540 45 50 3.6 O 3.3 / x anuary, =450 E 3.0 - 2.7 g 2.14 ~/2.1-// July, ~.600...."" 1.8 / / * * /I...o..o 1.5 1.2 Annual, 4,75@ 0.6 0.3 I 0 5 10 15 20 25 30 35 40 45 50 Uo (msec') Figure 5. Imaginary component of phase velocity as a function of zonal wind at p = 0 for adiabatic model-January, p. 45 0; July, cp =. 30~ and 600; and annual average, cp = 750, 1963.

middle latitudes and be of very small magnitude in low latitudes. Perhaps the most startling feature of Figure 4 is the pronounced maximum in T at 55~N for January which is hardly expected from Figure 3. The explae nation lies in the relative magnitude of two effects working against each other — T Ocos qp and T Ccl/c (m/sec) as seen in Eq. (3.1). Thus, as the North Pole e e i is approached the first effect, one of pure geometry, tends to decrease T while the second is inverse to the magnitude of ci from Figure 3. The fact that there is a secondary minimum in ci at 55~N dominates the latitudinal effect and we have a maximum in T. For q) > 55, the latidudinal effect domine ates (decreases more rapidly than) ci, inducing a tendency toward smaller T 1i~~~~~~ e The reverse takes place for cp < 55~. For this reason there is a strong resemblance between Figures,3 and 4 for January, July, and the-annual average in low and middle latitudes. We note also in Figure 4 that T (July) exceeds at every latitude T (annual average) and at all but 50~N and 55~N, T (January), e e indicating greater stability during the summer, which is seen clearly from Figure 3. Figures 6, 7, and 8 illustrate the inverse relationship of T to n for January, July, and the annual average profiles. In each the shape of the curves is very similar for all n. Thus, as the perturbation wave length decreasesthe instability as measured by T increases, with a minimum T of e e approximately four days for January, 7 1/2 days for July, and six days for the annual average, all of course, for the third harmonic. It is interesting to compare this with the results from Wiin-Nielsen's (1961) three-parameter

55 0 26 ~4 - IQ{...Q o..o 8 _..X.20 _..16 C,18 ~~ 20 22 26 2130 32 34 36 38 40 20 25 30 35 40 45 50 55 60 65 70 75 80 85 LATITUDE Figure 6. e-Folding time as a function of latitude for January, 1963adiabatic model at n = 1, 2, and 3.

8 33 n3 10 12.. 14- / 18 - X 20 _.. 022e 24 28 30- n 32 34 x 36 38 2020 25 30 35 40 45 50 55 60 65 70 75 80 85 LATITUDE Figure 7. e-Folding time as a function of latitude for July, 1963adiabatic model at n =- 1 2, and 3.

57 O j, 2 4 6 _... --..... 3 10,##X~~ m..'1 ~.. n 3 8 **0000 —4c 14 " 16 14 18 20 22 24 26 28 30 32 34 36 38 40~I I I 20 25 30 35 40 45 50 55 60 65 70 75 80 85 LATITUDE Figure 8. e-Folding time as a function of latitude for annual, 1963 — adiabatic model at n = 1, 2, and 3.

58 model having linear vertical wind shear and increasing static stability with decreasing pressure, where he predicts a doubling-time of the order of two to three days for n = 1, 2, and 3. When advection of relative vorticity is included in his vorticity equation, he predicts a minimum e-folding time of about five days for the meteorological range of vertical wind shear, very much in agreement with results from this study. In many baroclinic instability studies in which meridional variation of instability is not considered, characteristics of the instability are portrayed in what is called a stability diagram in which contours of T are drawn on a graph having vertical wind shear e as the ordinate and wave length as the abscissa. On this a curve representing wave length of maximum instability is often constructed, generally a function of vertical shear. Due to the meridional variation of vertical shear in the real atmosphere which we consider, it is not possible to construct a stability diagram representing significant information. Suffice it to say that at any latitude, of the three wave numbers we consider, wave number 3 has maximum instability. In comparison with this representation of instability by T we recall that instability for adiabatic flow as measured by ci is independent of wave number. It follows that Figure 3 represents n = 1, 2, and 3. In determining eigenvalues for Newtonian heating, cN, from those computed for adiabatic flow, c*, we recall q a cos cp c c a c i (2. 86) *N *a - s n U6) o For q > 0 and s U > o in westerlies or easterlies it is clear that (c N)i < (c. ),' indicating that heating tends to stabilize the flow. The degree of

59 stabilization is directly proportional to q, the heating or cooling coefficient, and cos cp, while inversely proportional to wave number and amplitude of basic zonal velocity. Figures 9, 10, and 11 compare c. for the adiabatic and Newtonian models in the three cases-January, July, and the annual average, respectively. We immediately note the similarity in shape between the adiabatic and Newtonian waves for each case. As wave number increases, the Newtonian curves approach the limiting adiabatic curve, crowding closer together. Further; whereas the minimum ci in adiabatic flowwasnon-negative, for Newtonian heating ci can be negative, indicating damping or stable flow at particular latitudes and wave lengths. Maximum damping occurs of course for n = 1 and approximates c. = -2.4 m/sec at 20~N for January, -2.25 m/sec at 20~N for July, and -2.3 m/sec at 25~N for the annual average. In comparing the adiabatic and Newtonian curves we observe that the maximum stabilizing effect of heating occurs when the degree of instability is maximum, while the difference between adiabatic and Newtonian instability diminishes with increasing latitude to less than.5 m/sec at 850N. Clearly, all values of T for Newtonian heating exceed e corresponding values for adiabatic flow. For ci < 0, T cannot be defined. Its upper limit of infinity is approached as c. + 0 from above. Figures 12, 13, and 14 illustrate, for Newtonian heating, the variation of c. with latitude for the three cases —January, July, and the annual averi age-at each wave number. As with adiabatic flow c. (January)> ci (July), Ci (annual) for all n and all latitudes except 50~N and 55~N where ci (January) dips sharply. Despite the variation of ci with n the curve shapes for each case are strikingly similar, approximately independent of wave number.

60 6Adiabatic!i \v ^\\ Newtonian 3-f \ =2 b 11qn 3 E2 2 I- 0o 3 I 1 I January, 1963-adiabatic model; and Newtonian model at n = 1 2, and 3 -3___________________________ 202530354045 50556 6570 5808'1 x ~'~'......IT... Figure~ ".Iaiaycmoeto hs eoiya ucino a~u~ o Janary 193~aiaai moel'".Newonamoeatn=1,2ad 3...

61 2.5 2.0 2.0 oAdiabotic 1.5?/?..... Newtonian 15J^ -x-~ x.....ol...1.5 1 I / 1)%./'....Q I -1.5 X- + I I/.x -2.0 1 -2.5 Figure 10. Imaginary component of phase velocity as a function of latitude for July, 1963-adiabatic model; and Newtonian model at n = 1, 2, and 3.

62 5. 4 3 Adiabatic Newtonian I Ss X ^^ X *~~~^^1. "^^TL _? 2 - -4/ 20 25 30 35 40 45 50 55 60 65 70 75 80 85 LATITUDE anuual aversge, 1963 —adiabatic model; and Newtonian model at n = 1, 2, and 3.

63 4l- x I \ 21 \January - I Annual.D......D! I..a July -2 - -3 -4-5 -.! i; I I I I I i... 20 25 30 35 40 45 50 55 60 65 70 75 80 85 LATITUDE Figure 12. Imaginary component of phase velocity as a function of latitude for Newtonian model at n = 1.

64 5 41,,- I \ I I \January 21 Annual..El..... I E 0 -2 -3 -4 -.5, I I, I I.... I,,I,,I I,, I_ I, _. 20 25 30 35 40 45 50 55 60 65 70 75 80 85 LATITUDE Figure 13o Imaginary component of phase velocity as a function of latitude for Newtonian model at n = 2.

65 51~^ I \January 4 j I \ I I 2 - | /Annual f~/ /~Jul 33 -2/.' -X"x I -4 I z -5-5 I! I I I I I! 1. 1. 1._ 20 25 30 35 40 45 50 55 60 65 70 75 80 85 LAT ITUDE Figure 14. Imaginary component of phase velocity as a fU:ll1iti-: of latitude for Netonian model at n = 3.

66 In discussing the dependence of instability on wave length (or wave number) for assumed qp and U (cp), if c. is used as a measure of instability, we 0 1 recall from Figure 3 no dependence for adiabatic flow. This is in agreement with Wiin-Nielsen's (1961) three-parameter model having zonal wind increasing linearly and static stability increasing with decreasing pressure. From Figures 9-11 instability or ci increases as n increases or wave length decreases for Newtonian heating. If we select T as our instability parameter, (3.1) expresses the decrease of T (or increase of instability) as n increases or wave length decreases for both adiabatic flow and Newtonian heating for given cp and U (p). 3.2. THE PHASE SPEED We define the phase speed, c, as the real part of the phase velocity of the perturbation wave. It is the speed with which the wave travels eastward relative to the rotating earth. It is not the speed of the wave relative to the basic zonal velocity, U. A wave is termed retrogressive if c and U are of opposite sign, i. e., if the wave moves in a direction opposite to that of the flow in which it is embedded. Based upon statistical studies by Eliasen (1958) and Martin (1958) of motion based on zonal harmonic analysis of observed patterns of geopoter;ial in the troposphere, planetary or ultra-long waves appearing in the westerly flow of middle latitudes are slowly moving systems which oscillate around certain preferred geographical positions, moving to the east or west. Eliasen shows that at 500 mb and 50~N for wave numbers 1, 2, and 3, the deviation of

67 the wave from its mean position exceeds one-fourth of the wave length in only a few cases. Occasionally displacement during one day is relatively large, but this always coincides with small values of amplitude. As a consequence of the semi-permanence these waves are maintained on monthly mean charts as well as on normal charts for individual months, which indicates that the mean positions of the waves are nearly the same from year to year. Motion of waves at 40~N and 60~N parallels rather closely motion at 50~N. Hirota (1968) suggests from these studies that there are two kinds of ultra-long waves in the middle-latitude westerlies. One is a quasi-stationary or forced standing wave and the other is a traveling or free wave. Due to the absence of any geographical or terrain features in this study tending to induce standing waves, we conclude that the phase speeds computed refer to the traveling waves. Figure 15 displays the relationship of phase speed to latitude for January, July, and the annual average. The curves represent all wave lengths of the ultra-long waves, both for adiabatic flow and Newtonian heating. This is apparent when we recall that c* for adiabatic flow is independent of wave number and, from (2.86), that the real parts of c*N and c* are identical. Only their imaginary parts differ by a factor dependent upon wave number. Comparing Figures 15 and 2 we note a striking resemblance between U vs. cp and c vs. cp. Maximum values of phase speed are 13.9 m/sec for January at r 35~N, 6.7 m/sec for July at 45~N, and 8.2 m/sec for the annual average at 40~N. For the most unstable waves as determined by ci, phase speeds are 8.6 m/sec at 25~N in January, 5.2 m/sec at 40~N in July, and 6.3 m/sec at 35~N

68 14 I I I 13/ \ 12 / \Jonuary Io,I I \ 97 \ 6 I6S 6' ^ X12XIl X X mr:July /X I 2 I \ - - I' \' 3-4r\! 025 3035 40-45 505560-65~7075 8085 LAT ITDE -J__________________ 20 5 0 3 04 05 06 07 08 -2 / ~ ~~LAITD Figure 1^. Phas~~espeaa 1011oflfcLfucefoaca'afianNwtin moes~auay Juy'n"ana vrae 93

69 for the annual average. These values are in good agreement with those obtained in the empirical study by Bradley and Wiin-Nielsen (1968) where the phase speeds ranged from 5-15 m/sec. While most of the waves move eastward relative to the earth, some do not, i.e., the waves in January at 20~N, in July at 20~, 25~, 30~, and 850N, and for the annual average at 25~N. Of these only the waves in January at 20~N, in July at 30~N, and for the annual average at 25~N are retrogressive. The maximum speed to the west is 5 m/sec. In Appendixes A. 1 and A. 2 phase speed of stable perturbations in the barotropic limiting cases, A = A (p), (1)a = a = constant and (2) a = o /p is determined analytically and it is found in both models that the waves move at a slower eastward (or faster westward) speed than the basic current at any latitude, the tendency toward westward motion increasing with decreasing latitude. High latitudes evidence little discrepancy from the basic current. Comparing the models one finds that inclusion of inverse pressure dependence of static stability increases the tendency toward westward motion by a factor of about 3. In contrasting these limiting cases with our results as presented in Figure 15, where the basic zonal wind profiles are shown in Figure 2, it is clear that, while the shapes of the curves in Figures 2 and 15 are similar for the respective January, July, and annual average cases, the magnitudes of phase speeds are in general considerably less than the basic zonal winds, most of which are westerlies. The vertical scale shows this most graphically. This completely agrees qualitatively with the barotropic model results. To obtain a quantitative estimate of this tendency, we tabulate U (Figure 2) - c -1 (Figure 15) in m sec as follows;

70 c (deg) 20 30 45 60 75 January 30.4 31.2 14.6 8.6 10.7 July - 3.6 3.8 14.0 8.3 5.3 Annual average -- 17.0 16.0 9.3 5.1 This reveals, when compared with Tables Al and A2, baroclinic tendencies toward westward motion considerably greater than those in the limiting barotropic cases with the sole exception of 20~N in the second barotropic model at mode number 1. We note decreasing discrepancies at higher latitudes in consonance with the barotropic models. In the quasi-geostrophic 10-layer model of Hirota (1968), having TJU/6p = 2 m/sec/100 mb and constant static stability the amplifying waves with X > 5000 km, which he terms the second type of unstable waves, have phase speed almost independent of wave length, c - U = -5 m/sec. Here U is the basic zonal velocity at 500 mb and is thus, 10 m/sec, giving c = 5 m/sec. In our notation, the U corresponding to this would be 20 m/sec. The closest 0 vertical wind shears we have to this occur in July at 40~N, where U = 18.5 m/ sec and for the annual average at 30~N, where U = 20 m/sec. The phase speeds 0 for each are 5.2 and 3 m/sec, respectively, which is in good agreement with Hirota. While phase speed is very nearly constant with wave length for U = 0 20 m/sec, Hirota finds that, as the vertical wind shear increases in magnitude, c becomes less constant with wave length, increasing as k increases. This is r exemplified in his 20-layer model where, for U = 40 m/sec, c = 12 m/sec for = 12,000 km and 1 /sec for = 20,000 km. In comparison, we find, for k = 12,000 km and 13 m/sec for k = 20,000 km. In comparison, we find, for

71 U = 39 m/sec in January at 25~N, c 8.6 m/sec, not as good agreement as in O r the former comparison. 3.3 THE WAVE STRUCTURE If our model is at all an approximation to reality, the unstable waves which are solutions to the linearized equations should have a vertical structure which is in a gross sense comparable to the structure of ultra-long waves in the atmosphere. In particular we investigate the vertical structure of the most unstable waves as determined by a maximum in c.. We recall that these occur, for adiabatic flow and Newtonian heating, at 25~N in January, 40~N in July, and 35~N for the annual average at n = 1, 2, and 3. One of the striking results of this study is the pronounced similarity in the structure of the waves, especially in regard to the vertical variation of amplitude, and approximately for the vertical variation of phase angle. The assumed form of the perturbations, from Chapter II, is ( ), = ( A ) ein(x-ct) (2.21) A where ( ), a function of cp and p, describes the variation in perturbation amplitude and is in general complex. Given eigenvalue c, due to the fact that ( )', multiplied by any constant will satisfy the perturbation equations, we choose to normalize the amplitude function at given qp by dividing by the maximum amplitude of ( ) for all p. Let us denote this by I ( ) f and the max normalized amplitude function and perturbation, ( ) - (... N I( )m max

72 ( ), = ( ) N ein( -ct) ( ^N NN For a = Re[ ( )N] b Im[ ( )N], we obtain, where c = c + ici, Re[ ( )N] = exp(ncit). [a cos n(X-c t)-b sin n(X-c t)] which represents the -N i r r physical quantity. It is an easy matter to show that Re[ ( )N] = A exp(nc.t) ~ cos (X - 8) N i where X = n (X-c t) r A = a2+b2 = (^)N and = tan ( Thus, at time t = 0, Re[( )] = A cos(nX-3) For given latitude, amplitude A and phase angle 8 depend solely on pressure. It is clear that, on a given latitude circle and isobaric surface, 5 is the value of nA at which the real component of the normalized perturbation will be a maximum or, according to some authors a ridge. In describing the wave structure we plot A and 5 vs. pressure for the normalized perturbation geopotential, vertical velocity, and temperature. It should be noted that there are several ways we could present phase information. As stated above cos(n\-5) is maximum when n\ = 6or k = -. For n = 1 the position of the ridge is specified by 8, while - and - represent ridges for the higher wave numbers 2 and 3. In addition for n = 2 a ridge is located at + 180~ while for n = there are addition, for n = 2 a ridge is located at - + 180o, while for n = 3 there are 2

73 two additional ridges, at - + 120~ and - + 240~ around a latitude circle. ReaS 3 sonably enough for wave number n we have n ridges. Analogously, we could depict troughs, where cos (n-6b) is a minimum. Here, n\ = 6 + 180~ or \ = 6/n + 180~/n. Thus, for wave number 1 a trough is displaced 180~ from the ridge, while for wave numbers 2 and 3, troughs are located 90~ and 60~ from ridges, respectively. For wave number n, there are n troughs. Because knowledge of the vertical distribution of 6 is sufficient to determine ridge and trough positions for any wave number, we elect to use it to display phase information. Other writers choose to plot ridges and/or troughs. This latter selection does have the advantage of presenting the cellular structure of perturbations around a latitude circle which is most evident for n.> 1. In this study it happens that 5 is virtually independent of n, enabling us to represent all phase information for a particular parameter in one figure rather than three, one for each wave number. Eliasen (1958) explains the existence of the long, semi-permanent waves as resulting from the combined effort of the large-scale distribution of heat sources and large mountain barriers, i.e., the Rockies and Himalayas. These will effect both the phase angle and amplitude of the ultra-long waves. Considering first the phase angle, in the winter continents largely act as cold sources and oceans as warm sources, while the reverse holds during the summer. This should cause a phase shift for corresponding waves from January to July. From observations we find in the lowest layers of the troposphere phase

74 differences to be considerable between January and July. This qualitatively agrees with the phase shift due to heat sources. However, in the atmospheric layer from 700-500 mb phase differences between winter and summer become much smaller. This constitutes an essential argument for ascribing an important role to large mountain barriers which would not tend to induce any phase differences during the year. Effect of diabatic heating and cooling on the amplitude of ultra-long waves is included in Derome's (1968) study on the basis of quasi-geostrophic, linearized, steady state equations. Here, the most prominent feature, a trough, was found to be about three times deeper at 25 cb than at 75 cb. Perturbations caused by mountain barriers have been investigated for the equivalent-barotropic level by Charney and-Eliassen (1949) and Bolin (1950) and considerable agreement with the real 500 mb perturbations demonstrated. Eliasen (1958) reports that by extending the computations to a baroclinic atmosphere and utilizing the confluent hypergeometric functions tabulated by Charney (1947), it is found that amplitudes of orographically-induced waves increase with height, the amplitudes at the 500 mb level being between two and three times the 1000 mb level amplitude. To conclude, large-scale distribution of heat sources is important for the formation and structure of the quasi-stationary, ultra-long waves in the lowest part of the troposphere, while the effect of large mountain barriers is predominant in the middle and upper troposphere. As discussed, however, in an earlier section there is another kind of ultralong wave, the traveling or free wave. Our neglect of large-scale heat sources

75 and mountain barriers in the models studied indicates that it is the structure of this second kind of wave disturbance which we examine. Figures 16-21 portray the vertical variation of the phase angle and amplitude of normalized geopotential, vertical velocity, and temperature, respectively, for the most unstable waves in January, July, and for the annual average in the adiabatic and Newtonian models. As alluded to earlier there is striking similarity between the profiles, so much so that seasonal variations and changes in wave number are deemed negligible in each model, a single curve on each plot representing January, July, and the annual average for n = 1, 2, and 3. Further, amplitude deviations between the models are effectively infinitesinal, leaving us with just one curve on the amplitude figures and two on the phase angle graphs, where the differences between the adiabatic and Newtonian models are large enough to be significant on the plots. In Figure 16 depicting phase angle of geopotential, the most striking feature in the east-west tilt with decreasing pressure is the maximum rate of change for.3 < p* <.5 in both models, with less rapid change occurring in the lower and upper atmosphere. As shown in Appendix A. 5 the existence of a westward tilt in the most unstable waves indicates northward sensible heat transport. The range of 8 is approximately 220~. This compares very favorably with an empirical study on the transient part of the ultra-long waves undertaken by Bradley and Wiin-Nielsen (1968) in which the tilt of surface harmonics of both forced and free modes of geopotential was calculated to be westward from 850 to 500 mb and approximately vertical from 500 to 200 mb.

0.1 \ Adiabatic.2 -.3.4 0..6 Newtonian'..7.8.9 1.0 I I I I I 40 60 80 100 120 140 160 180 -160 -140 -120 -100 -80 -180 PHASE ANGLE OF 4 Figure 16. Phase angle of perturbation geopotential of most unstable wave as a function of pressure for adiabatic and Newtonian models.

77 Extrapolating their values in the 850 to 500 mb regime for the free mode to the rest of the atmosphere we find the range of 6 to be 233" for n = 1 and 2 and 156~ for n = 3, in good agreement with our computed value of 220~. Phase angle for Newtonian heating is about 5~ west of that for adiabatic flow. This agrees wellwi-th.Hirota's (1968) waves in his second (long waves and strong shear) and third (weak shear) unstable regions, where waves exhibit a striking westward tilt in the middle layers, and the upper and lower half of the wave are out of phase with each other. Figure 17 depicts the vertical variation of normalized geopotential, the single curve representing adiabatic and Newtonian heating conditions for either January, July or the annual average and is independent of wave number. Here we note a maximum in the wave's intensity at the top of the atmosphere, a secondary maximum at the ground, and a minimum around p* =.4. We thus see in comparing Figures 16 and 17, that the wave exhibits amplitude maxima where the tilt is least, and an amplitude minimum where the phase angle change in the vertical is greatest. These amplitude characteristics agree well with Hirota's (1968) unstable wave solution for weak shear conditions (.4 m sec / 100mb) for X = 5000 km, referred to as the third unstable region, in his 20layer model, while there is less qualitative resemblance with his 10-layer model unstable wave solutions under strong shear (3.0 m sec /100mb) where the waves with X = 10,000 and 15,000 km have monotonically increasing amplitudes from the ground to the top of the atmosphere. We recall that in both of these cases there was good agreement with our phase angle variation.

0.2.3.4 0a..5.6.7.8.9 1.0 I I I I I I I I I I I I 0.05.10.15.20.25.30.35 40 45.50.55.60.65.70.75.80.85.90.95 1.00 I$/t lmax Figure 17. Normalized amplitude of perturbation geopotential of most unstable wave as a function of pressure for adiabatic and Newtonian models.

79 Wiin-Nielsen's (1961) three-parameter model including the advection of relative vorticity displays vertical amplitude variation very similar to our own but having variation with wave number. Finally, in Derome and Wiin-Nielsen's (1966) quasi-geostrophic model having constant static stability and Coriolis parameter one finds normalized amplitude of unstable stream functions or geopotential waves (since they are assumed porportional) maximum at p* = 0 and 1 and minimum at p. = 1/2 with virtually no wave number dependence, excellent agreement with our own results. The phase angle variation in this model also agrees well with ours with the exception that the tilt increases with increasing wave length. In Figure 18 we plot phase angle of vertical velocity, o, vs. pressure. Here as with., both.the tilt is greatest in mid-troposphere and the phase angle for Newtonian heating is about 50 west of that for adiabatic flow, but the range of 8 is very much less, approximately 40~, indicating a slight eastA A west slope of c compared with D. Derome and Wiin-Nielsen (1966), in their constant static-stability and Coriolis parameter quasi-geostrophic model alluded to above, found a very similar vertical variation, the primary disparity being that the range of 5 increased with increasing wave length, approximating 60" at 18,000 km. As with our study, the slope was nearly constant as pressure decreased, having a slight maximum in mid-troposphere. They concluded that, although the effect of static stability is to decrease the slope of the waves, this effect is felt to a smaller degree as the wave length increases such that the slope of the ultra-long waves is not affected appreciably.

0, Adiabotic.2.3 -.4 c*.5 *' Newtonian'.6.8 1.01 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 PHASE ANGLE OF b Figure 18. Phase angle of perturbation vertical (pressure) velocity of most unstable wave as a function of pressure for adiabatic and Newtonian models.

81 A Turning to normalized amplitude of c in Figure 19, we observe a variation very nearly opposite to that of 0. Here, minima occur at the bottom and top of the atmosphere to satisfy the lower and upper boundary conditions, respectively, while the maximum appears at p* =.4. Derome and Wiin-Nielsen (1966) find the same structure, with a maximum, however, at p* =.5 and negligible variation with wave length. They conclude, as above, the static stability of the basic flow does not affect appreciably the variation with pressure of the unstable X wave amplitude. Like the phase angle variation of $ and c, we observe in Figure 20 a A maximum slope of T in mid-troposphere, with Newtonian heating phase angle 5~ west of that for adiabatic flow. The range of 5 is 85-90". In contrast with A A X, T has a maximum tilt for.5 < p. <.7, below the layer of maximum slope for ~~~~~~A A ~. The variation of the normalized amplitude of T is seen in Figure 21, where there is a maximum at the bottom boundary and minimum at the top of the atmosphere. This is reasonable due to of the fact that T p from the ideal gas law. Further, we observe a secondary maximum around p. =.4, which, from Figure 17, is the level of minimum in 1l1. Combining the perfect gas and hydrostatic equations, we obtain T "cr/~p or,alternatively,T o M/6(ln p) which supports this behavior.

0.1.2.3.4.6.85.. 7....9!.O - _ I I, I, I i I I i i i i. 0.05.10.15.20.25.30.35 40 45.50.55.60.65.70.75.80.85.90.95 1.00 I /I Cimax Figure 19. Normalized amplitude of perturbation vertical (pressure) velocity of most unstale wave as a function of pressure for adiabatic and Newtonian models.

0 2. kAdiabatic.2.3:.4- \,\. *.5-..6- Newtoan -..7-.8 1.0 1 1 1 I! 1 1 1 1 1 I 1 1!! I X i.9- 1. 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 PHASE ANGLE OF T Figure 20. Phase angle of perturbation temperature of most unstable wave as a function of pressure for adiabatic and Newtonian models.

0.1.2.4.5.6.7.8.9 I.0.9 i I I I I...I. I i I I I I I! I i I 0.05.10.15.20.25.30.35.40.45.50.55.60.65.70.75.80.85.70.95 1.0 ITI / ITimax Figure 21. Normalized amplitude of perturbation temperature of most unstable wave as a function of pressure for adiabatic and Newtonian models.

CHAPTER IV ENERGETICS 4.1. GENERAL CONSIDERATIONS In the study of unstable waves it is of paramount importance to consider the energetics or energy transfer processes involved, in particular how the amplifying perturbation obtains its energy. The foundation of most current energy studies was laid by Lorenz (1955), included in which he considered energy conversion processes of unstable long waves associated with maintenance of the general circulation. He divided the atmosphere's energy into four forms: zonal kinetic energy (K ); eddy kinetic energy (K); zonal available potential energy (A ); and eddy available potential energy (A ). Kinetic z E energy (KE), expressed as one-half of the square of the speed of a fluid 1 2 parcel of unit mass, v, is the energy the parcel possesses as a consequence of its motion. Potential energy is the energy which a body possesses as a consequence of its position in the field of gravity and is given by gz for a parcel of unit mass. Internal energy, written as c T for a parcel of unit mass of perfect v gas, is defined as the sum of the electromagnetic potential energy of electrons and nuclei and the kinetic energy due to the molecular motions. Under hydrostatic equilibrium Haurwitz (1941) has shown that the potential and internal energy of an atmospheric column are proportional. Therefore in an hydrostatic atmosphere it is convenient to combine these into one form, the sum being called total potential energy. Lorenz (1955) introduced the concept of tj

86 available potential energy (APE), defining it as that portion of the total potential energy which may be converted to kinetic energy in a. system across whose boundaries no heat or mass is transported. The minimum state of total potential energy is characterized by horizontal stratification of the potential temperature surfaces. Any parameter a(\,cp,p,t) can be expressed a.s the sum of the zonal average, a (c,p,t), and eddy component, aE(k,c,p,t), where a is defined as 1 2t a = - ad (4.1) z 2- o Lorenz (1955) showed that the energy transfer scheme for an atmosphere in baroclinic instability in its normal state is; Az K C(Az,AE) Here, C represents an energy conversion quantity and has the units of time rate of change of energy; thus, C(A,AE) is the conversion from zonal to eddy APE and C(AE,KE) is the conversion from eddy APE to eddy KE. Potential energy of the basic flow (A ) is therefore converted into kinetic energy of the unstable perturbation (Ks) through the intermediary, AE. This stands in

87 marked contrast with the barotropic instability mechanism, where kinetic energy of the basic current (K ) is the energy source and is converted directly into z KE. The mechanism for energy conversions in baroclinic instability is related to the three-dimensional structure of the unstable wave. This implies a, definite relationship between the results in Chapter III and in this chapter which will be examined in following pages. Saltzman (1959) suggested that ultralong waves receive part of their kinetic energy from small-scale motion in addition to the baroclinic mechanism referred to above. The absence of nonlinear interactions in this study precludes examination of his thesis. 4.2. CONVERSION FROM ZONAL TO EDDY AVAILABLE POTENTIAL ENERGY, C(A,AE) For an hydrostatic atmosphere it can be shown that p a" ((aO"v) cos m zc(AA) p1 i a~ dS dp (4.2) C(AAE) p S G a cos qp (p for a layer of atmosphere between pressures p1 and P2 bounded by the latitude limits pL and pU, where cp = 20~N, U = 85~N for the January and July cases, L U'LU while pL = 25 N, U = 80oN for the annual average. Here S represents surface area, a differential element of which satisfies 2 dS = a cos dp d cp d while ( ) = area average of ( ) over a given isobaric surface and ( )" = ( )( ), the deviation of ( ) from its area average. Because all perturbation parameters may be expressed in terms of eigensolution c, by normalizing with respect to the maximum of 131 over all cp and p*, and then calculating the other

88 A perturbation parameters in terms of the normalized a, we are enabled to compare, not just the signs of C(AZ,AE), but also the magnitudes for the adiabatic and Newtonian heating models and their respective seasonal cases. The term in (4.2), (a"v)z, can be re-expressed and simplified by recalling that v = (2.2) 2 Q a sin cq cos cp ax from which we find vz = 0. It follows that (za"v) = (0w Va + ivE) Eliminating the first term by [( ] = 0, we arrive at (o"v) = (c vE)z which, on an isobaric surface, is proportional to the meridional transport of sensible heat by the wave perturbation, as these eddy parameters are synonymous with the perturbations in Eqs. (2.11), denoted there by primes. We retain the ( )E notation to conform with current practice in energetics. Graphs of normalized C(Az,AE), denoted CNO(AzAE), vs. pressure are presented in Figures 22-29 illustrating the dependence upon the model, either

89 adiabatic or Newtonian heating, the season, and the wave length of perturbation. The expression in Eq. (4.2) is numerically integrated for a given isobaric surface over the entire range of latitude and then computed for ten layers in the vertical, each of 100 mb thickness. In order to eliminate the exponential dependence upon time in the perturbation parameters we have taken t = 0, thus giving instantaneous values of CNO(Az,AE). Figures 22 and 23 depict the adiabatic and Newtonian models in January for n = 1, 2, and 3. The difference is striking. CNO(Az,AE) under adiabatic flow is < 0 for p* <.88, all but the lower troposphere where there is conversion from Az to AE. The minimum at the top of the atmosphere results from the inverse dependence of cz on pressure, while the secondary minimum in ICNO(Az,AE) around p =.55 may be explained A by the secondary minimum in normalized ITI from Figure 21. The similar shapes of CNO(AZAE) for n = 1, 2, 3 follow from the equivalence of the zonal or basic state contributions for each model and wave length in (4.2), while 2 2 A A 2 Q a sin qp 0 0 = in~ a~-(2.22) in p

.0O ~~.3~~ ~ ~x: n=3:-9.4-10 n=2/ C(Ax 102) of pressure for January.5-adiabatic model at n = 1, 2 and.6 n =I.' l\ i i I J_"_ 0.7.8.9 1.0. 10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 CNO( Az, AE) (x10'2) Figure 22. Normalized conversion from zonal to eddy available potential energy as a function of pressure for January. 1963 — adiabatic model at n = 1, 2, and 3.

0 g~E-~~~ i^~ C NO(Az, AE) (xlO2) Net:2~^ i'~') S n= 1 36.6.2 n=2: 18.3.3 n = 3: 12.2.4 c*.5.6.70 l =2 n =3" - 1.0 \ __ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 CNdAZ, AE) (xIO2) Figure 23. Normalized conversion from zonal to eddy available potential energy as a function of pressure for January, 1963 —Newtonian model at n = 1, 2, and 3.

92 requires that eddy meridional heat transport be inversely proportional to the wave number, since a and vE are proportional to derivatives of (E (or' ), while vE, from Eqs. (2.2) and (2.21), is also directly proportional to n. For all wave lengths the summation of CNO(Az,AE) through all layers is < 0 for adiabatic flow, indicating transfer from AE to Az. The reverse is true in the Newtonian model, where the magnitude of CNO(Az,AE) approximates that for adiabatic flow, but the energy flow is from Az to AE. Figures 24 and 25 portray CNO(Az,AE) vs. p* for July in the two models at n = 1, 2, and 3. In both there is negligible conversion above approximately 400 mb and C(Az,AE) < 0 for all lower layers. The difference in shape of the curves is very noticeable, approaching asymptotically a limiting value of CNO(AzAE) with increasing pressure for adiabatic flow, while in the Newtonian model increasingly larger negative values of CNO(AZ,AE) result. In spite of this asymptotic behavior the greater magnitudes of CNO(Az,AE) in the adiabatic model at all pressures result in net coonversions fromAE to A more than an order of magnitude larger than for Newtonian heating. Comparing the adiabatic

.I- CNO(Az,AE) (xlO')Net n= 1:-107 n =2: -53.3 n =3: -35.4 n=..6- / /.. x n2/ 1.08 I. A I -- I- I -22 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 CNo(AZ, AE)(x IO) Figure 24. Normalizei conversion from zonal to eddy available potential energy as a function of pressure for July, 1963 —adiabatic model at n = 1, 2, and 3.

0.1 CNo(Az,AE)(xl2) Net.2 n=1: -5.00 3- n=2: -1.53 n 3: -0.70.4 *.5.6.7 n=1 B _~x....8.s~ " ^-"^ ^. ~~* 00-0 ~ ~ ~.~ n 3 -265 x D- ^ 2' 1.0 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0 CNo(AZAg) (xlO ) Figure 25. Normalized conversion from zonal to eddy available potential energy as a function of pressure for July, 196")~Newtonian model at n = 1) 2, and 3.

95 winter and summer cases in Figures 22 and 24, we see an increase in the magnitude of the total conversion from AE to A of some three to four times from E z January to July, while Figures 23 and 25 for Newtonian heating reveal for each wave number a change in both the direction of energy transfer and in the magnitude. While at all levels conversion in January is from A to AE, in July the z E' reverse is the case, at a rate five to ten times less than its winter counterpart. In each model the northward meridional heat transport and amplitude of a" in January usually exceeds that in July. z Figures 26 and 27 reveal the variation of C NO(A,AE) with pressure for the adiabatic and Newtonian heating models, respectively, when the annual average of U (cp) is considered. Interestingly, there is a strong resemblance between each and Figure 22 for the January adiabatic case. Here, however, the pressure level dividing the regimes is somewhat higher than the p* =.88 in Figure 22. For the annual-adiabatic model CNO(A,AE) < 0 for p* <.83, while in the annual-Newtonian model, CNO(A,AE) < 0 for p* <..76. Significantly, the net conversion in the annual-adiabatic model is from AE to A for all wave numbers E z while, on a scale an order of magnitude smaller, the annual-Newtonian model energy conversion is from A to AE. Although this represents an instantaneous z F energy conversion, it does agree with the time-averaged behavior of the atmosphere and supports the validity of inclusion of Newtonian heating in the numerical modelling of ultra-long waves. Figures 28 and 29 combine for n = 1 the graphs appearing in Figures 22-27 for the adiabatic and Newtonian heating models, respectively. These illustrate more clearly the agreement and disparity between C N(A,AE) for January, July,

NO (A A )(xI02)Net n I: -15 n 2: -7.3 / n 3: -5.4 *'5 a~~~~~ oo ~.6 t Ro -7 6 5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 CNO(AZAE)(KIO of pressure for annual, 1963-adiabatic model at n = 1, 2, and 3.

CNO(AZAE)(xlO2)Net n = I: 1.25.*~~31~~~ C 1n2:.63.4 3 n=3.43.6.7.8 -.9 i s.... -.2 17' 1.0 -.6 -.5 -.4 -.3 -.2 -.1 0.1.2.3.4.5.6.7.8.9 1.0 1.1 1.2 1.3 CNO(A, AE)(xlO2) Figure 27. Normalized conversion from zonal to eddy available potential energy as a function of pressure for annual, 1963-Newtonian model.

.1 /.2 CN(AZA E(xlO)Net /x cNO,A)(x102e XO January: -28 X January/ July:-107 J ana.4 Annual:-15 /... " ct.5 - ~~...~~July 4J.. —X.7~..6 7 -x.8 1 rS 9 ~ ~~~~~~~~~~~~~ ~~Annual 7.0 I I -22 -19 -16 -13 -10 -7 -4 -1 0 2 5 8 CNO(Az AE)(xlO2) Figure 28. Normalized conversion from zonal to eddy available potential energy as a function of pressure for adiabatic model at n = 1.

0. g CNO (A, A )(xIO -2)Net.2 ~~~~~~~~~~~~~~O E.2 t January: 36.6.3 July: -5 < r~~~~~~ 3 ^ Annual: 1.25.4 5- ^ CL.6 ^0 jb^^^January y ~88 Q^9 ^ Annual 1.0 I L -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 CN AE n)(XIO-2 Figure 29. Normalized conversion from zonal to edddy available potential energy as a function of pressure for Newtonian model at n = 1.

100 and the annual average in each model. It is of interest that the annual variation of CN (A,AE) is considerably more similar to the January than the July curve in the adiabatic model while the reverse is true for Newtonian heating. Similar curves for n = 2 and 3 are not constructed, for their shapes would closely resemble their n = 1 counterparts with a diminution in amplitude as can readily be seen in Figures 22-27. 4.3. CONVERSION FROM EDDY AVAILABLE POTENTIAL TO EDDY KINETIC ENERGY, C(A K E) For an hydrostatic atmosphere it can be shown that p2 (4.3) C(AE',K) - f fS dS dp (4-3) for a layer of atmosphere between pressures p and p2 bounded by the latitude limits pL and U' As stated earlier ( )E is equivalent to the primed perturbation, ( )', as expressed in Eqs. (2.11). Because WE and o are physical quantities we take only their real parts in ( )' - ( () e (2.21) for t = O0 ( ) = ( ) + i ( ) It is now an easy matter to show E = 0 from integrating cos nX and sin nk around a latitude circle, implying = E. By expanding w in the z-!coordinate system we find dp _ p + p up (4p dt at ax + y w ( where w - dz/dt, the vertical velocity in the z-system. In these coordinates the geostrophic wind approximation as expressed in (2.1) and (2.2) becomes

101 pf py u = -a 1 ap V = P pf ax where p = density, from which the horizontal pressure advection in (4.4), u ap/ax + v ap/dy = 0. Assuming for geostrophic flow ap << ap an assumption well supported by observations, (4.4) reduces to w z.' (4.5) az' In the z-system the hydrostatic relation has the form ap which, substituted into (4.5), gives C — gpw or, for the perturbation, E - -gPzE (4.6) where p represents the basic state density. z Applying the perfect gas law to the perturbation we have R TE aE = p

102 which, combined with (4.6), gives E aE E - WE TE on an isobaric surface. It follows from (4.3) that C(AE,KE) cJS WE TE dS on an isobaric surface. Thus, a positive correlation between wE and TE, i.e., warm air rising and cold air sinking, causes energy conversion from AE to K. As with C(A,AE) we normalize C(AE,KE) with respect to the maximum of IS| over all p and p and then calculate the other perturbation parameters in terms of the normalized X so that we may compare both sign and magnitude of normalized C(A EKE), denoted CNo(AE,KE), between the adiabatic and Newtonian heating models and their respective seasonal cases. Graphs of CNo(AE,KE) vs. pressure are presented in Figures 30-37, illustrating the dependence upon the model, either adiabatic or Newtonian heating, the season, and the wave length of perturbation. The expression in (4.3) is numerically integrated for a given isobaric surface over the entire range of latitude and, as was done with C NO(A,AE), then computed for ten layers in the vertical, each of 100 mb thickness. Taking t = 0 in all calculations gives us instantaneous conversions as in the former case. Figures 30 and 31 depict the January-adiabatic and Newtonian models for n = 1, 2, and 3. It is immediately apparent that there is a strong shape similarity in the curves between the two models. In each CNO(AE,KE) < 0 for p <.56 while conversion is from AE to KE for p* >.56. The comparable variation in each model of CNO(AE,KE) for n = 1, 2, and 3 follows as it did for CNO(A,AE), from the inverse dependence of $ on wave number as seen in (2.22). In this case,

/ n=l:-1.2 o2 > /d 2n 2:-0.6.3 n =3:-0.4.4 ~.5'-5 -3-2-1 012345678 CL.6t- 1.7.8 2 CNo(AE,KE)( O) Figre 30. Normalized conversion from eddy available potential energy to eddy kinetic energy as a function of pressure for January, 1963-adiabatic model at n = 1, 2, and 3.

0 ~I.a 1" CNO(AE,KE)(xIO2 )Net 0 -.2 n=1:-0.21 3L g:ic I n =2:-0.10 03 3 ( L n =3 -0.07.4 V^ x a 5 —.6- n l.7- \.\n-=2 I.0 L I I.. I.. i..!. i a x.8 3/.9 -.6 -.5 -.4 -.3 -.2 -.1 0 1 2 3 4 5 6 7 CNO(AE KE)(xIO ) Figure 31. Normalized conversion from eddy available potential energy to eddy kinetic energy as a function of pressure for January, 1963-Newtonian model at n = 1, 2, and 3.

105 however, only aE is a function of i in (4.3), giving C (A EKE) C 1/n. This E NO E' E A holds strictly for adiabatic flow, where w is independent of n, but the dependence of cu on n in the Newtonian model implies a more complicated relationship between the energy conversions and n. The amplitude thus decreases as wave length decreases for CN (AE, ) as was evident with CN (A,AE). In contrasting the two models it is apparent that the effect of including Newtonian heating is to reduce the amplitude of the AE-KE conversion by approximately an order of magnitude at each layer and summed over the entire atmosphere. The net transfer for both is from KE to AE, but is small compared with the amplitude of CN (AE,KE) at an arbitrary level. This indicates a type of balance between the removal of KE in the upper atmosphere and production of it below. Figures 32 and 33 illustrate the variation of CN(AE,KE) for the Julyadiabatic and Newtonian models. It is immediately clear that there is not nearly as much similarity between them as during January, the shapes resembling closely the curves for CN(A,AE) in the same models. While there is virtually no energy transfer above 350 mb under adiabatic conditions, the Newtonian model does exhibit non-negligible conversion in the upper atmosphere compared with the lower troposphere, especially for n = 1. Like CNO(A,AE) there is an asympototic approach toward a limiting value at each n in the adiabatic model as the lower surface is approached, while amplitudes under Newtonian heating become increasingly greater. At all levels conversion is from KE to AE, with adiabatic amplitudes between one and two orders of magnitude larger than corresponding ones for Newtonian heating. Figures 54 and 35 represent the annual-adiabatic and Newtonian curves and

~ CNO (AE, KE)(xlO'2)Net 1.I n= -9.2 n =2:-4.6 *2 n=3 -3.1.3 0.4.5.6-.7 =2 in=3 o.,.8 1 0 t.9 1.0 I. J.. -1.8 -1.6 -1,4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0 CNO(AEKE )(x 2) Figure 32. Normalized conversion from eddy available potential energy to eddy kinetic energy as a function of pressure for July, 1963-adiabatic model at n = 1, 2, and.

O ~.I CNO(A, KE)(xlO-2)Net n = 1:-0.30 n =2:-0.09.3 n=3: -0.04.45 (3..5-.6.7.8 n=2 X00 0~'n:/.9-15 -.13 -.1 -.09 -.07 -.05 -.03 -.01 0 CNO(AE, KE)(xlCT2) Figure 33. Normalized conversion from eddy available potential energy to eddy kinetic energy as a function of pressure for July, 1963-Newtonian model at n = 1, 2, and 3.

0.1 r ~ CNO(AE. KE)(xlO')Net.2 = / n =1:0.5 3^ ^ ^ ^ n 2: 0.3 ({x ^ Qf3n: 0.2.4 0. O0.6 n=3 Ef IWOOOO "" ^^^~ ^^o^^n=l I.7.8 fl3..E /X.9 -, 1.0 I I -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 1820 CN0 (AE., KE)(xI0y2) Figure 34. Normalized conversion from eddy available potential energy to eddy kinetic energy as a function of pressure for annual, 1963-adiabatic model at n = 1, 2, and 3.,

0 )(3.1- ^^^^~ CNO(AE,,KE)(xlO2)Net ^^.2 n=I:0.01 ^ ^ ^ n 2: 0.004.3~ I~j n =3: 0.0025.4 n.5~~~~~~~~~~~ n=.6 \ Jn "X~~~.7 n ~8~~~~~~~~~ 1 n=~~~~~~~~~~l-1011.9'Oe 1.0 I.-__ _ _ _ _ __ __ __ _ 6 -5 -.4 3 -.2 -.1 0.1.2.3.4.5.6.7.8 CN(AEl KE)( X1O2) Figure 55. Normnalizea conversion from eddy available potential energy to eddy kinetic energy as a function of pressue for annual, 1963-Nevtonian model at n =, 2, and 3.

110 appear very similar to their January counterparts in Figures 30 and 31. In each CNo(AE, ) < 0 for p* <.56 while the reverse is the case below. Again the adiabatic amplitude considerably exceeds the Newtonian, by approximately a factor of 25. Despite the likeness with the January curves there is a very important disparity, in that the conversions for the annual average are of opposite sign summed over the entire atmosphere. In both the annual-adiabatic and Newtonian cases, there is an instantaneous net transfer from AE to KE, agreeing with the time-averaged behavior of the atmosphere. Summarizing our results for the annual-adiabatic and Newtonian models over the whole atmosphere, we have as instantaneous energy diagrams the following; Annual-Adiabatic Annual-Newtonian ~ CNO(AzAE) 0 CNO(AzAE) > 0 AE | E~ CNO(AEKE) > 0 CNO(AE,KE) > 0 -2 -2 Annual-Adiabatic (x 10 ) Annual-Newtonian (x 10 ) n C z (AA) C N(AE KE) dAE/dt n CO(AAE) CN (AEKE) dAE/dt NO z'E NO E'E E NO z(E NO'E E) 1 -15 0.5 -15.5 1 1.25 0.01 1.24 2 - 7 0.3 - 7.3 2 o.63 o.oo4 o.626 5 - 5 0.2 - 5.2 5 0.45 0.0025 0.4275 where dAE/dt = CNO(AzAE) - CNO(AE,KE). The dominance in the magnitude of CNO(Az,AE) over CNO(AE,KE) agrees with Murakami's (1962) results from scale analysis that only the interaction between zonal and eddy available potential energy is predominant for ultra-long waves, exceeding all other energy conversions by two orders of magnitude. We thus conclude that for each wave number there is an instantaneous loss of Am in the annual-adiabatic model which exceeds by approximately an order of -hi

111 magnitude the instantaneous gain of A in the annual-Newtonian model. ComE paring the two diagrams we immediately arrive at the conclusion that the instantaneous energetics of the Newtonian model conforms much closer to the time-averaged state of the atmosphere depicted in Section 4.1 than its adiabatic counterpart. This is strong support for the inclusion of Newtonian heating in an- ultra-long wave model. It should be noted, however, that the energetics of the adiabatic model cannot be considered inconsistent with the time-averaged state of the atmosphere, for it refers to a certain instant of time and there is no requirement for a steady-state condition. Finally, Figures 36 and 37 combine for n = 1 the curves appearing in Figures 30-35 for the adiabatic and Newtonian models, respectively. These illustrate more clearly the agreement and disparity between CN(AE,K ) for NO E E January, July, and the annual average in each model. While we have already alluded to the greater similarity of the annual to the January rather than the July curves, it is evident comparing Figures 36 and 37 how very much closer this agreement is in the Newtonian model. In contrast the fact that there is such a divergence between the annual- and January-Newtonian cases in C (A,AE) for n = 1, as seen in Figure 29, only emphasizes that the two energy conversion parameters, C(A,AE) and C(AE,KE), describe very different physical mechanisms. Similar curves for n = 2 and 3 are not constructed for the same reason as given in presentation of CNo(A,AE). NO z E 4.4. COMPARISON OF ENERGY STRIP AND FULL ENERGY DIAGRAMS For the purposes of this study we define a full energy diagram as one

F0 JX/ I CNO(AE,KE)(xlO2) Net'X/.2 - / January: -1.2.3! July: -9.2 Annual: 0.5.4- * H.6I 6.| iAnnual.7 - 7 i January.8 July /.9 1.0 I I I I -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 CNo(AE, KE)(xl2) Figure 36. Normalized conversion from eddy available potential energy to eddy kinetic energy as a function of pressure for adiabatic model at n = 1.

O. 1-,' CNO(AE,KE)(x1O )Net.2 f January: -0.21 July: -0.30 Annual: 0.01 5.1P hi 5. Annual.7.8 July/ January/,' 1.0 I -.6 -5 -.4 -3 -2 -.1 0.1 2.3.4.5.6.7.8 CN(AE KE )(x102) Figure 37. Normalized conversion from eddy available potential energy to eddy kinetic energy as a function of pressure for Newtonian model at n u 1.

114 having the three energy quantities, A,AE and KE, and consider the instanz E E taneous net (integrated over pressure) conversions CNo(A,AE) and CN(AE,' NO z %E NO E JE over the entire range of latitude as has been undertaken in Sections 4.2 and 4.3. Figures 22-37 pertain to these diagrams. While it is of interest to know the integrated energetics over all latitudes predicted by our models, because of the fact as mentioned earlier that our models become less valid in the higher latitude regime, especially for wave number 3, it is desirable to calculate the net CNO(A,AE) and CN (A K ) in thin latitudinal strips around the NO z'E NO E'E latitude of maximum instability for the January, July, and annual average cases and define these as energy strip diagrams. We recall that these latitudes are 250, 40", and 35~N, respectively, assuredly well within the region of the models' applicability. We consider strips of width 5~ on each side of these latitudes. A further motivation in performing these calculations is to determine the instantaneous energetics in the region of maximum instability in each case and to compare with the integrated energetics. The comparative diagrams for both models follow for wave numbers 1, 2, or 3. Strip Energy Diagram Full Energy Diagram I. Adiabatic model A. January A j IA D z I E g y \ E I~I E ~ ~ KE.04

Strip Energy Diagram Full Energy Diagram B. July A I IA zl z~~~~i z AwK A K iE ij Ly~1j E i i E i ~ E 0o C. Annual average Az A 13 EAE [^M^] E E E.45 II. Newtonian model A. January Ai Az 2 o015 B. July A I IA.02 0

116 Strip Energy Diagram Full Energy Diagram C. Annual average Z Z Here, arrows indicate the direction of energy transfer, while the numbers on the strip energy diagrams represent the magnitude of the particular strip energy conversion relative to a unit value taken for the analogous full energy conversion quantity. The number, 0, represents a negligible magnitude. With the sole exception of C (A,K ) 0 in the January-Newtonian model, this NO E E supports the conclusion that the instantaneous energetics in both models in the region of maximum instability is in agreement with the time-averaged energetics as observed in the atmosphere and predicted by Lorenz (1955). In resolving the apparent discrepancies between the strip and full energy diagrams it is helpful to discuss the energetics of an arbitrary unstable wave perturbation. From the nature of instability the total energy, AE + KE, of the wave must increase, or, equivalently, C(A,AE) > 0. While this is true for all strip calculations, it is false for all the adiabatic as well as the July-Newtonian full energy diagrams. The fact that these diagrams include the entire latitude range and that, at certain low latitudes, there is no instability, i.e., C(A,AE) < 0, invites the conclusion that the negative contributions of C(A,AE) from the stable regions exceed the positive contributions from the unstable regions, even though the latter includes most of the latitude range.

117 Interestingly, vertical distribution of each strip energy conversion is independent of the model or case considered. It is found that CNO(A,AE) is positive at every level and increases monotonically with increasing pressure as shown in the sketch, P. O + CNO(AzAE) On the other hand the shape of CN (A E,KE ) in the strip calculations closely resembles its full" energy counterparts in the adiabatic and Newtonian January or annual average cases, being negative for pK <.55 and positive below as seen in the diagram, p * 55 ~75 CNo(AEIKE) +

CHAPTER V CONCLUSION 5.1. SUMMARY OF RESULTS We have examined the internal dynamics with regards to stability, structure, and energetics of ultra-long, transient, baroclinic wave perturbations superimposed on a meridionally-varying zonal current, assumed linear in pressure, in adiabatic, frictionless, geostrophic flow and how the addition of Newtonian heating modifies the results. Three separate cases are investigated based upon input data~ of the zonal wind between 20~N and 85~N during January, July, and the annual average, 1963. It should be emphasized that these data represent atmospheric behavior during a single year, 1963, and do not necessarily depict climatology. We cannot, thus, assume that the detailed results of this study apply for the normal January, July, or annual atmosphere. Applying the perturbation theory a single second-order ordinary differential equation for adiabatic flow is derived and numerically solved by the "shooting method" in a 40-layer model for the normally complex perturbation phase velocity, c. From this the corresponding quantity for Newtonian heating is easily determined. Knowledge of the eigenvalue, phase velocity, suffices for complete analysis of the structure and energetics of the disturbance in the two models. A discussion of the uniqueness of the unstable eigenvalues is included. Investigation of instability as measured by the imaginary component of 118

119 phase velocity, ci, reveals no dependence upon wave number, n, for adiabatic flow, while instability increases as wave number increases for Newtonian heating, which itself exerts a stabilizing effect. This effect is maximum at maximum instability. For the same input zonal wind distribution adiabatic flow is always more unstable than its Newtonian counterpart. Another measure of instability, e-folding time, is inversely proportional to wave number for either model. At all wave lengths in adiabatic and Newtonian flow greater instability is found over most latitudes in January than in July or for the annual average, having a maximum at 25~N in January, at 40~N in July, and at 35~N for the annual average. In the adiabatic model this corresponds to e-folding times for the third harmonic of 4, 7 1/2, and 6 days, respectively. This is not surprising in view of the fact that the zonal wind in January exceeds that in July and for the annual average at all latitudes except 45~ < c < 60~ and that vertical wind shear tends to make waves unstable. The absence of instability in adiabatic flow occurs only at low latitudes, in particular 20~N in January, 30~N in July, and 20~and 250N for the annual average. Slight instability is associated with zonal easterlies which occur only in July at 20~N, 25~N, and 85~N. The possibility of negative values for ci, which indicates damping, arises only in Newtonian heating where it is maximum at low latitudes for n = 1. The phase speed, defined as the real part of the phase velocity of the perturbation wave, is equivalent for both adiabatic and Newtonian flow and is

120 independent of wave number. Its variation with latitude closely resembles the input distribution of zonal wind with latitude. Maximum values of phase speed are 13.9 m/sec in January at 35~N, 6.7 m/sec in July at 45~N, and 8.2 m/sec for the annual average at 40~N, while the most unstable waves travel at 8.6 m/sec in January (25~N), 5.2 m/sec in July (40~N), and 6.3 m/sec for the annual average (35~N). While most of the waves move eastward relative to the earth, some do not, i.e., the waves in January at 20~N, in July at 20~, 25~, 30~, ard 85~N, and for the annual average at 25~N. The maximum speed to the west is 5 m/sec. The vertical wave structure of only the most unstable waves is investigated and is found strikingly similar for the adiabatic and Newtonian models at n = 1, 2, and 3 for all seasonal cases. Wave number and seasonal dependence is negligible in every case as well as amplitude variation between the models. Only phase angle, 5, which tilts westward with decreasing pressure for every parameter, reveals any distinction, those for Newtonian heating being approximately 5~ west of those for adiabatic flow at each level. The most striking feature in the vertical phase angle variation of perturbation geopotential is the maximum east-west tilt between 300 and 500 mb. The range of 5 over the entire atmosphere is approximately 220~. Examination of the normalized amplitude of geopotential reveals a maximum at the top of the atmosphere and a secondary maximum at the ground, both where the tilt is least, and minimum around 400 mb, where the tilt is greatest. Phase angle variation of the vertical pressure velocity shows, like the geopotential, maximum tilt between 300 and 500 mb

121 but with a much reduced range of b -approximately 400-indicating a slight east-west slope of vertical pressure velocity compared with geopotential. Variation in normalized amplitude of vertical pressure velocity is very nearly opposite to that of geopotential. Minima appear at the bottom and top of the atmosphere to satisfy the boundary conditions, while the maximum occurs at 400 mb. Phase angle of perturbation temperature, which has a range of 85-90~, displays maximum tilt below that of the other parameters-in the interval, 500-700 mb. Normalized amplitude of temperature, owing to the proportionality between temperature and pressure from the ideal gas law, has a maximum at the bottom boundary and minimum at the top of the atmosphere. A secondary maximum appears around 400 mb. Instantaneous.energetics of the adiabatic and Newtonian models is stud ied from two standpoints-(l) by calculation of vertical variation and total over the entire range of pressure and latitude of the normalized conversion from zonal to eddy available potential energy, C O(A,AE), and normalized conversion from eddy available potential energy to eddy kinetic energy, CNo(AE, KE), and (2) by calculation of vertical variation and total over all pressures in a latitudinal strip 100 wide centered on the latitude of maximum instability of the same energy conversions. The A - A -K diagrams integrated over pressure z E E pertaining to the former we refer to as full energy diagrams while that of the. latter we call strip energy diagrams. A comparison between these diagrams for the adiabatic and Newtonian heating models at each wave number may be found in Section 4.4. With the sole exception of CNO(AE,KE)< 0 in the January-Newtonian model, the strip diagrams support the conclusion that the instantaneous energetics

122 in both models in the region of maximum instability is in agreement with the time-averaged energetics as observed in the atmosphere and predicted by Lorenz (1955). In resolving the apparent discrepancies between the strip and full energy diagrams-where C(A,AE) > 0 for an unstable wave in order for its total energy, AE + KE, to increase —we conclude that the negative contributions of C(A, AE) from the stable regions in the full energy diagrams exceed, for all the E adiabatic cases as well as the July-Newtonian case, the positive contributions from the unstable regions, even though the latter includes most of the latitude range. Comparing the full energy diagrams of the annual adiabatic and Newtonian models, it is clear that the addition of Newtonian heating causes the instantaneous energetics to conform more closely to the time-averaged observed energetics over the entire atmosphere. Though this distinction does not appear in the strip energy diagrams it, nonetheless, does provide some support for the inclusion of Newtonian heating in ultra-long wave models. While there is definite dependence over all latitudes of the vertical distribution of each energy conversion quantity with the model and case considered, the strip energy calculations reveal virtual independence around the latitude of maximum instability. It is found for these strip conversions that CN (A,AE) is positive at every level and increases monotonically with increasing pressure, while CNo(AE,KE) is negative for p <.55 and positive below for each wave number in the January, July, and annual average adiabatic and

123 Newtonian heating models. Variation of the energy conversions with wave number-CN(A,AE) 1 /n and CNo(A,KE)' 1/n is shown analytically, holding explicitly for adiabatic flow and approximately for Newtonian. 5.2. CONCLUDING REMARKS AND SUGGESTIONS FOR FUTURE WORK This study pertains solely to the transient baroclinic mode of the ultralong waves. The existence, in addition, of a stationary mode and a transient barotropic mode suggests that an empirical investigation be conducted in the ultra-long wave regime to determine the division between the three modes. The discussion of the uniqueness of the unstable eigenvalues is not a rigorous mathematical proof, for the size of one of the neighborhoods in which uniqueness is valid is indeterminable from the information available, while the size of the other neighborhood is estimated from results of the numerical searches. Thus, while Markushevich's (1965) implicit function theorem appears suitable as a vehicle for proving the desired uniqueness, the problem of rigorously determining the size of the neighborhoods in which uniqueness is valid was not solved. Accomplishing this would be most desirable and would plant the results of this study on a more solid footing. From the conclusions in this study one is encouraged to include a heating term in a model for ultra-long waves. Certainly a more realistic heating mechanism than Newtonian heating should improve our study in this regime. Inclusion of large-scale external heat sources and sinks as pioneered by Smagorinsky (1953) would be a long step in this direction.

124 Charney and Eliassen (1949), Smagorinsky (1953), and Derome (1968), in their studies of the effects of large-scale external forces on the existence of stationary long waves, showed the importance of friction. In the simplest inclusion of friction one would start with a two-, three-, or four-level model. This could be extended to a continuous friction model where both surface and internal friction are included, at the expense of a more difficult lower boundary condition. A possible technique of handling this would be to shrink the friction layer such that we assume its top is at p. = 1 and begin iterating upwards from this point rather than starting at p* = 0 and iterating downwards as was done in this study. In this way we would hope to use the "shooting method" in eigenvalue searches. A more efficient search procedure than the "shooting method" such as Newton's method or an interpolation scheme could be investigated. Further, large-scale topographical effects, though they might require empirical expressions, would contribute to a more realistic model. Other generalizations would be inclusion of a meridional scale for the perturbations and a perturbation expaions of the zonal velocity in powers of the Rossby number which would enable us to make corrections to the zerothorder approximations.

APPENDIX A.1. THE BAROTROPIC LIMITING CASE, A = (cp),a = = CONSTANT Assuming a barotropic basic state and constant static stability reduces Eq. (2.24) to the standard form 2A 2 A a (+ k = 0 2 ap where 2 a p 2 o00 k 2 (A. 1) 2Qa sin cp(A -c) which has the general solution A L = A cos k p + B sin k p for arbitrary constants A and B. Satisfaction of the upper boundary condition, (p O) = 0, requires A = 0 while the lower boundary condition, (p = 1) = 0, specifies the solution to be A a = B sin k p where k = m, m = + 1, + 2,... (A.2) 125

126 Combining (A. 1) and (A. 2), 2 a p ( ^ -1^ ^o o c(rad sec ) A= A- o 222 2 20a m t sin cp Multiplying by a cos cp, we obtain 2 C P c(m sec ) = U - 2 2 s U o 1 2a am ti sin cp which indicates that the waves are always stable and move at a slower eastward (or faster westward) speed than the basic current at any latitude. The difference, cl depends upon the latitude and the mode number, m, as shown in the following table. TABLE Al DIFFERENCE BETWEEN BASIC CURRENT AND WAVE SPEED FOR BAROTROPIC LIMITING CASE, a = a CONSTANT z o cp (deg) 0 15 30 45 60 75 90 cl m=l oo 15.72 3.77 1.54 0.73 0.31 0 cl, m=2 oo 3.93.94 0.38 0.08 8 0 c, m=3 o 1.75 0.42 0.17 0.08 0.03 0 It is clear from this table that the tendency toward westward motion increases with decreasing latitude and decreasing mode number. In high latitudes there is little difference from the basic current regardless of mode number.

127 A.2. THE BAROTROPIC LIMITING CASE, A = A (p), a =p O Z Assuming a barotropic basic state and static stability inversely dependent upon pressure reduces eq. (2.24) to the form 2A 2 a cu k A 2+ - = (A.3) 2 where Eq. (A. 1) specifies k Equation (A. 3) has the same form as the canonical equation 2q-2 y" - cx = o which has the solution, from Kamke (1961), pg. 440, y= I l/2q ( where B is an arbitrary constant and Jl/2q is the Bessel function of first A 2 kind and order 1/2 q. Letting y =, x = p, c = -k, and q = 1/2 we obtain the solution to (A. 3). = B J1 (-2k P*). J being antisymmetric about the origin the solution becomes 1 A =W -jiB P J2k^P) In order o satisfy the lowers boundary condition, (p* = 1) = 0, we must require Jl(2k) = 0

128 The zeros of J1, from Jahnke and Emde (1909), are a = 0, a =+ 3.83, 2 =+ 7.02, a = +10.17,... giving k = a /2, m = mode number = 0, 1, 2, 3,...'5 m Thus (A. 1) may be expressed 2 p2 -1\ T o o COscpI c(m sec ) = U _ - cosi U= - c. o a2 2 o 2 QaaC sin cp m As in A. 1, where static stability was constant, this indicates that the waves are always stable and move at a slower eastward speed (or faster westward speed) than the basic current at any latitude. This discrepancy, c2, depends upon the latitude and mode number as shown in the following table. TABLE A2 DIFFERENCE BETWEEN BASIC CURRENT AND WAVE SPEED FOR BAROTROPIC LIMITING CASE, Cz = co/p* qp (deg) 0 15 30 460 75 90 cl, m=l 0 42.5 10.2 4.16 1.97 0.82 0 c,' m=2 oo 12.6 3.03 1.24 0.58 0.24 0 c, m=3 oo 6.02 1.45 0.59 0.28 0.12 0 It is clear from this table that, as in A. 1, the tendency toward westward motion increases with decreasing latitude and decreasing mode number. Similarly in high latitudes there is little discrepancy from the basic current regardless of mode number. Comparing Tables Al and A2 we conclude for these barotropic cases that inclusion of inverse pressure dependence of static stability increases the tendency toward westward motion by a factor of about 3.

129 A.3. COMPARISON OF THE QUASI-GEOSTROPHIC AND GEOSTROPHIC MODELS FOR A BASIC STATE OF NO MOTION AND CONSTANT STATIC STABILITY We assume in Eqs. (2.11) that U(p,p) = 0 and a = C = constant. z o I. Quasi-Geostrophic Model. Expressing the first two perturbation equations of motion in Cartesian geometry, the geostrophic equations, (2.1) and (2.2), generalized to include local accelerations become -= + -f Vi' (A. 4) 6v a t f u' (A.5) at y where f = 2Qsinq, while the corresponding perturbation forms of the continuity eq., (2.4), and energy equation, (2.5), are u -v' c' u + V~ + ~ 0 (A.6) ax by ap a a + F' = 0 (A.7) at ap o Introducing perturbations of the form, analogous to (2.21), ( ) =( ) e ik(x - ct) where ( ) is assumed to be a function of y and p and k = 2T/(wave length), (A.4) - (A.7) become -ikcu -ik + f v (A.8) = v (a. 8)

130 A 3 A -ikcv = - - u (A 9) A a iku + + = 0 (A. 10) -ikc +ct o = o (A. 11) op o Solving (A.8) and (A. 9) for u and vwe obtain.^ (k2ct + f ) A ay A ik k t v = f( + c 2 2 2 where A f k c which, substituted into (A. 10) gives p -, (k c f r ) + i k (f + c -/ ) 6p A Cy dy A Solving for a in (A. 11) and substituting yields the single equation c a ^ ) + ( ^- ^ = (A. 12) c ^ - Ay A AJ o subject to the boundary conditions A c) = 0 for p = 0, p or, equivalently, from (A. 11), a = o for p = o, p (A. 13) ap

131 Equation (A. 12) is in general very difficult to solve because of the variation of A and f with y. One may, however, employ a beta-plane approximation which, according to the usual practice, is f = f + y (A. 14) o where f is the value of the Coriolis parameter at a central latitude,o, on o 0 -4 -1 the spherical earth and is usually taken as 10 sec and D is the Rossby -12m-1 -1 parameter, 20 cos cp/a, and is commonly given the value 16 x 10 m sec. It is ignored in (A.12) unless it appears with constant factors. Performing the differentiations with respect to y in (A.12) after inserting (A.14) one gets,% 2Pf~$ l[ (f + k22) k2 0 pa 1 af p 2 A 6y -c L A o ay Ignoring P because it appears with nonconstant factors eliminates the term containing ~, while restricting our considerations to wave motions for which 22 2 k c < < f, i.e., to waves whose speed is much smaller than the speed of inertia waves, reduces the equation to 2 1 6 62 B 2A f2 ( - -) + (+ = 0 o ap o 6p 6y2 c where A is approximated by f. o This equation may be solved by separation of variables, i.e., (y, p) = G (y) F (p) (A.15) which transforms the equation to

132 d 2 (d dd 2 G dpo dp_ dyd 1 F 2 f G 0 The first term being solely a function of p and the second of y implies that each is equal to a separation constant, l/gh, where h will be defined shortly and the form of the constant clarified. Setting d( 1 dF dp\Go dp 1 F gh and 2 dy2 (c 1 f2 G gh 0 we find - -- + - F = 0 (A.16a) dp V\ dpj gh

133 and 2 f2 d' ( + k2 + ) G = 0 (A. 16b) dy2 c gh Non-dimensionalizing pressure we obtain from (A.16a) 2 dF + F 0 2 gh dp* which has the solution F = A cos a p* + B sin O p. (A. 17) where 2 op 2 00 (A. 18) gh and (A.15) into (A.513) gives the boundary conditions dF = 0 for p. =, 1 dp' Differentiating (A. 17) with respect to p*, we obtain dF ~d = a (B cos oqp - A sin a p.) dp. which, combined with the upper boundary condition, gives B = O. Repeating for the lower boundary condition we have = m, m = 1, +2,... which, substituted into (A.18) requires that

134 2 o p h = 2 2 (A. 19) gm if This relation explains why the separation constant was expressed as l/gh, for h has the dimensions of a length or a depth. It is normally called the equivalent depth. f2m2~2 2 2 2 fom Letting r = - k + ) (A.20) - p o o after applying (A.19) the solution of (A.16b) becomes G = D cos r y + E sin r y (A.21) subject to the boundary conditions A v = 0 for y = 0, d, equivalent to vertical walls around the Equator and latitude circle specified by y = d. These may be reexpressed, using (A.8) and (A.9) as A 14 u = - c and A o ay f u + = 0 which combine to give f y _ _ o 0 for y = 0, d c

135 or, equivalently, after inserting (A. 15), IG 0 0 dG= 0 dy c Substituting (A. 21) requires (1 D + r Ecos r y + ( E -r D sin r y = 0 for y = 0, d from which the condition at y = 0 gives us f -D +r E 0 or D =- r E c f o and that at y = d yields r d = nr, n = + 1, + 2.. Inserting this condition into (A. 20) and solving for phase speed gives us at last 2 222 1 22 k +fmJT + -2 n it o d 2 cp o o Taking -6 2 _ 2T x 10 3 k = ~L = 10 where 1 is in 10 km., = 2 o p = 100 o m = 1 = n d = 107 in MTS units, we find

136 16 9.4784+.0.34~+-.o 5 0335 2 1 which gives TABLE A3 WAVE SPEED IN QUASI-GEOSTROPHIC MODEL FOR BASIC STATE OF NO MOTION AND CONSTANT STATIC STABILITY 1(103 km) c(m sec-1) 1(103 km) c(m sec-1) 2 -1.07 16 -3.08 4 -2.13 18 -3.10 6 -2.61 20 -3.12 8 -2.83 22 -3.13 10 -2.95 24 -3.14 12 -3.01 26 -3.14 14 -3.06 28 -3.15 II. Geostrophic Model. Neglecting the local accelerations in (A.4) and (A.5) we obtain the perturbation geostrophic relations in Cartesian geometry, fv' = ax fu= - by which transform to A 1 a U - f by A ik v -= - f after assuming perturbations wavelike in x and t.

137 Substituting into the continuity eq., (A. 10), we obtain a: ik a. a S a~ -~ -+ ik T (~) = 0 A ap f by ky f which, after solving for w in (A. 11) and applying (A.14) yields the single equation a 1 a 1 at fmz - at 1 6 - - 0 = o f p CY0 p f by f2 This reduces to a (1 a) -o P oo cf 0 where f f ~- o or, equivalently, 2 p2 U* P 2 00 where y = - (A. 22) c f 0 subject to the boundary conditions, from (A. 13) 6A =- 0 for p = 0, 1 This problem is identical to the F portion of the ageostrophic model in I from which we extract the result, y = mc, m = + 1, + 2,... Solving for c in (A.22) we have

138 P 2 o p. 0 0 C 222 f m it 0 For = 16 x 10 a = 2 0 p = 100 0 -4 f = 10 0 m = 1 -1 in MTS units, we obtain c - 3.2 m sec Comparing this value with those in Table A3, and recalling that the length of the latitude circle at 45~N is about 28 x 103 km, we conclude that for a basic state of no motion at 45~, the percentage errors in assuming geostrophy for wave numbers 1 (1 = 28), 2 (1 = 14), and 3 (1 = 9-3) are 1.6%, 4.6% and 10% respectively. A.4. DISCUSSION OF THE UNIQUENESS OF THE UNSTABLE EIGENVALUE SOLUTIONS We desire to display the uniqueness of the unstable eigenvalue, c*, obtained using the "shooting method" in the solution of Eq. (2.29) and its associated boundary conditions for given cp and U (p'), p'< cp. It is shown in Section 2.1.3 that, with the introduction of a new independent variable, P* 1= -sc (2.29) is transformed to 2A A e(-i) 2 f-f u (cp') = 0 (2.31) 6 L o

139 where <P CP cp ^' L ^ IPL - and 2 CTo 0 E COS p c 0f -+ C, U sin cp' d p' (A. 23) 2MA a sin qp U sin cp L 0 0 0 A 1 under the boundary conditions = 0 at = 0, i-sc^ Eq. (2.31) is a special case of Gauss' differential equation, known also as the hypergeometric equation, which has the canonical form ( 2A 2 Equating (2.31) and (232), we arrive at y'=0 7 t -I = -1 = -f from which we easily obtain = -l+r (2.36a) -= -1 -r (2. 6b) 7 =0 (2. 56c) where r = 1 f (A. 24) From Kamnke (1961), pg. 1i67, the solution to (2.52) subject to the limiting

140 conditions specified by (2.36) is, for Igl < I, A A A = C. Al + C2 (A (2.38) where = 2F1 (r, -r; 2; () (2.39a) A 2 2 = (l-r ) ~~nf2F1 (r,-r; 2; t) + power series in 5 starting with "1" (2.39b) Here, Gauss' function is defined as (a) (b) n F (a,b; d; x) 1 + (2.40) 2 1 n=l (d) nI where (a) = a (a + 1) (a + 2)...(a + n-l) and a, b, d, and x can be real or complex, but d cannot be zero or a negative integer. Slater (1966) shows on pg. 4 the Gauss function to converge for IxI < 1 and diverge for Ixl > 1. For x = 1, convergence occurs if Re (d-a-b) > 0 and divergence for Re (d-a-b) < 0. When Ixi = 1 but x 1, the series is absolutely convergent when Re (d-a-b) > 0, convergent but not absolutely so when -1 < Re(d-a-b) < 0, and divergent when Re(d-a-b) < -1. For Re(d-a-b) = -1 convergence occurs if Re(a + b) > Re(ab) and divergence if Re(a + b) < Re(ab). Copson (1962) shows on pp. 251 and 252 that there exists an analytical continuation of F (a,b;d;x) into I 1 - xl < 1 and that this function has a branch point 21 at x = 1. Eqs. (2.39)'at the upper boundary, O = 0, give = 0 and a2 = 1.

141 A Since c = 0 here, C2 must be 0. To satisfy the lower boundary condition we seek e = such that 0 F (r,-r;2; ) = 0 (2.41) Based upon the numerical unstable solutions of (2.29) and its associated boundary conditions, together with (A.23), (A.24), and the identity, 1 o 1-sc* we display in Table A4 for the January, July, and annual average cases, the parameters r and e which must necessarily satisfy (2.41). It is clear that these unstable numerical solutions lie outside the g unit circle of convergence of F (r,-r; 2; e ). Thus, in seeking to show 21 o their uniqueness we require an analytic continuation of 2F(r,-r; 2; o ) beyond 1| I < 1. One possible approach is to express 2F(r,-r; 2; 1 ) as a combination of two linearly independent solutions of (2,32), valid in a region having some part in common with the 0 unit circle, but extending beyond it, i.e., Re ( o) > 1/2. The presence, however, of the parameter, d = 2, in the Gauss function precludes usage of more than one of Kummer's standard 24 solutions, all of which have the rather simple form of the product of functions of e and/or (1-g ) and a Gauss function. The second solution takes a conO O siderably more complex and unwieldly form, causing the calculation of the two constants of proportionality to be no mean task. A second' much less burdensome approach is to express the Gauss function

TABLE A4 r, ^ CALCULATED FROM NUMERICAL UNSTABLE SOLUTIONS (a) January, 1963 (b) July, 1965 (c) Annual Average, 1963 r 5 r 5 r 5, 250 1.58 (1.238,.235) 200 3.121i (2.309, -.096) 300 1.72 (1.11.21) 30 1.41 (1.395,.248) 25 2.80i (2.397, -.103) 35 1.50 (1.308,.246) 35 1.36 (1.471,.246) 35 1.83 (1.079,.163) 40 1.39 (1.418,.246) 40 1.39 (1.418,.246) 40 1.46 (1.34o,.247) 45 1.35 (1.455,.242) 45 1.48 (1.328,.246) 45 1.35 (1.460,.243) 50 1.32 (1.489,.236) 50 1.63 (1.197,.224) 50 1.33 (1.482,.239) 55 1.33 (1.4835.237) 55 1.63 (1.197,.224) 55 1.33 (1.481,.240) 60 1.34 (1.466,.241) 60 1.41 (1.395,.248) 60 1.33 (1.481,.238) 63 1.38 (1.424,.244) 65 1.28 (1.541,.228) 65 1.28 (1.543,.226) 70 1.30 (1.513,.234) 70 1.21 (1.644,.197) 70 1.22 (1.623,.206) 75 1.28 (1.540,.228) 75 1.15 (1.731.,.166) 75 1.24 (1.596,.212) 80 1.25 (1.580,.219) 80 1.11 (1.741,.140) 80 1.52 (1.502,.081) 85 1.07 (1.741,.134) 85 0.35 (1.612, -.091)

143 as a Barnes-type integral. From Slater (1966), pg. 25, F ( 2. )-(2) 4_ r(r+s) P(-r+s) (-s)(,)S d 2F1 (ro 29 ir(r)r(-r) - o r(2 + s) ( s (A. 25) where larg (- )I< it. The convergence of the integral implies the analyticity of 2F(r,-r; 2; g ), when represented by (A.25), in the entire ~ -plane cut along the real axis from 0 to oo. Thus, when selecting a ~ -domain of analyticity, the only restriction is to confine it above or below the real axis. From Rainville (1960) F (z) is analytic except at z = 0 and all negative integers, where it has a simple pole, and z = oo, where an essential singularity occurs. For positive integer x, r (x + 1) = x, giving F (2) = 1, while Whittaker and Watson (1961) prove r (z) r (-z) - z We thus conclude from (A.25) that, in the r-plane, 2Fl(r,-r; 2; g ) is analytic for all r whose real components differ from 0 or a positive or negative integer. Selection of an r-domain of analyticity must accord with this condition. In proving uniqueness of the unstable numerical solutions we apply Markushevich's (1965) implicit function theorem, pg. 109, reexpressed in our notation: Let F(r, o) be a function of two complex variables which is analytic in a neighborhood Ir-r 0 < r, { - I <p of the point (r, ), and suppose that

144 F( ) = o, 1 o. (r0 00 Then there are neighborhoods N(r ) and N(l ) such that the equation 0 00 F(r,S )=0 has a unique root 0 = g (r) in N(O ) for any given r e N(r ). Moreover, the function t = 0 (r) is single-valued and analytic on N(r ), and satisfies the condition o (r) = 0." Clearly, we identify F(r,~ ) with the Gauss function, 2F(r,-r; 2; o). The point, (r, I ), may be selected, arbitrarily, from Table A4, from which it follows that F(r, ) = 0. The condition at L) 0 necessarily follows for non-trivial (r OO ) solutions of the second-order equation, (2.32), where F(r,i ) = 0 is 0 oo00 specified. If we select from Table A4(c) the solution at p = 35~ as r 1.50, (1.308,.246) and set r =.50, p =.246 to define the 00 domain of analyticity about (r, ), this domain will include all the 0 00 numerical solutions in Tables A4(a) and A4(c) and all but those at cp = 20~, 25~, and 85~ in Table A4(b), i. e., where the basic state zonal wind is easterly. The first two of these special cases cannot be treated by Markushevich, for Re(r) = 0, where F(r, 5 ) is not analytic. The case at (p = 85~ can, however, be taken separately by specifying r =.35, = (1.612,-.091), 3 =.35, and p =.091. Although the domain of 00 analyticity contains no numerical solutions other than (r, 0), Markushevich can still be applied. In determining the neighborhoods N(r ) and N(t ) we refer to Markushevich's proof of Weierstrass' preparation theorem, pp. 105109, where they are constructed as interiors of circles, I r-r I <' < f and 0 { o - 1oo < p' < p, subject to the conditions: 0o 00

145 (1) p': F(r,g ) has no zeros within {I - Io = p' except at the point g itself. 00 (2) d': Let p > 0 be the minimum of IF(r, e )I on the circle, I -5 1 = p'. Choose *' such that F(r,B ) - F(r, 0)I < o for all r within [r-r = r' and o on { - i, 0 0 0o = p'. We specify p' = p -.001 with the justification that, in computing, given r, 1 (or its equivalent) numerically, the search procedure encompassed the e -domain of analyticity, but yielded only the zero at =. Unfortu-domain of analyticity, but yielded only the zero at o = Unfortu0 0 00 nately, given the information available, Lt cannot be computed, making *' indeterminable. The size of N(r ) must therefore remain unknown and Markushevich 0 applied with that reservation and the hope that Ir-r i < s' does indeed encompass the other numerical solutions. The fact that, during the entire period of computation, in which the searches were often begun from different starting points, only one g was obtained for each r lends credence to the above ex00 O pectation. A.5. RELATIONSHIP BETWEEN WESTWARD WAVE TILT AND NORTHWARD SENSIBLE HEAT TRANSPORT We desire to show that, for geostrophic, hydrostatic flow there exists a necessary and sufficient relationship between westward wave tilt and northward sensible heat transport. For simplicity we consider a perturbation having wave number 1. On an arbitrary isobaric surface, p = p1, we select the

146 origin of a Cartesian coordinate system to coincide with the ridge of geopotential, giving us 1= cos k x (A.26a) 1 1 T = cos (kx + ) (A.26b) where l1 are constant, > 0. k = 2Jt/L L = wave length = length of latitude circle at latitude (p. C = phase angle of temperature Graphically, we have, assuming 0 <x < </2, N atr< P Pp \I\ // i \ 1\ W II 1 IE S -y-= \ / For H = northward sensible heat transport averaged over the latitude circle at p = pl

147 H c (T )z where v = S-N wind component at p! L ( )z = L o ( ) dx The geostrophic approximation on p = p1 requires that a1 1 1 V = f 6x giving H c (T x Applying (A.?6) we obtain H c - T (cos c0 cos kx - sin o sin kx) k sin kx which is equivalent to H c J k 1 sin a 2f 1 1 It follows that, on an isobaric surface, if geopotential leads temperature by no more than T/2 (more generally, t) as we have assumed and as is usually the case in the real atmosphere, there will result a northward sensible heat transport. Conversely, if geopotential lags temperature by no more than T/2 (more generally, aT), a southward sensible heat transport will ensue. Assuming hydrostatic equilibrium on p = p1, Eq. (2.3) becomes

148 RT1 pJ 1 P1 For isobaric surface p = P2 slightly above p ^ P, i.e., p2< P1 we take the forward finite difference to obtain = _ + R - T 2 1 P1 1 where Ap = P1 - P2 > 0 Assuming 2 = $ cos (kx + B) and replacing 1 and T from (A.26), the 2 2 1 equation becomes 2 cos (kx + 3) =' cos kx + PE cos (kx +). 2 1 p1 1 Expanding the cosines and equating coefficients of cos kx and sin kx, we now have 2 cos = - + Ra cos a 2 1 Pl 1 and 2 sin - R sin Dividing, it follows that RAp sin a tan P = 1 Pi 1

149 Thus, for 0 < ca < /2, tan 3 = +/+, implying 0 < 3 < c / 2 which indicates a westward tilt of geopotential. For -T/2 < a < O, tan D = - / +, implying -iT/2 < 3 < 0 which gives an eastward tilt of geopotential. We have thus shown that northward (southward) sensible heat transport and westward (eastward) wave tilt are, under geostrophic, hydrostatic flow, inextricably linked.

BIBLIOGRAPHY Bolin, B., 1950: On the Influence of the Earth's Orography on the General Character of the Westerlies. Tellus, 2, No. 3, 184-195. Bradley, J.H.S. and A. Wiin-Nielsen, 1968: On the Transient Part of the Atmospheric Planetary Waves. Tellus, 20, No. 3, 533-544. Burger, A. P., 1958: Scale Consideration of Planetary Motions of the Atmosphere. Tellus, 10, No. 2, 195-205. Burger, A. P., 1962: On the Non-Existence of Critical Wavelengths in a Continuous Baroclinic Stability Problem. J. Atmos. Sci., 19, No. 1, 30-38. Charney, J. G., 1947: The Dynamics of Long Waves in a Baroclinic Westerly Current. J. Meteor., 4, No. 5, 135-162. Charney, J. G. and A. Eliassen, 1949: A Numerical Method for Predicting the Perturbations of the Middle Latitude Westerlies. Tellus, 1, No. 1, 38-54. Charney, J. G. and M. E. Stern, 1962: On the Stability of Internal Baroclinic Jets in a Rotating Atmosphere. J. Atmos. Sci., 19, No. 2, 159-172. Copson, E. T., 1962: Theory of Functions of a Complex Variable, University Press, Oxford, 448 pp. Cowling, T. G., 1957: Magnetohydrodynamics, Interscience Publishers, New York, 115 pp. Deland, Raymond J., 1964: Travelling Planetary Waves. Tellus, 16, No. 2, 271-273. Deland, Raymond J., 1965: Some Observations of the Behavior of Spherical Harmonic Waves. Mon. Wea. Rev., 93, No. 5, 307-312 Derome, Jacques F., 1968: The Maintenance of the Time-Averaged State of the Atmosphere. Technical Report 08759-2-T, Department of Meteorology and Oceanography, The University of Michigan, 129 pp. Derome, J. F. and A. Wiin-Nielsen, 1966: On the Baroclinic Stability of Zonal Flow in Simple Model Atmospheres. Technical Report No. 2 06372-2-T, Department of Meteorology and Oceanography, The University of Michigan, 93 PP. 150

151 BIBLIOGRAPHY (Continued) DoYs, B. R., 1962: The Influence of Exchange of Sensible Heat with the Earth's Surface on the Planetary Flow. Tellus, 14, No. 2, 133-147. Eady, E. T., 1949: Long Waves and Cyclone Waves, Tellus, 1, No. 3, 33-52. Eckart, Carl, 1963: Extension of Howard's Circle Theorem to Adiabatic Jets. Phys. Fluids, 6, No. 8, 1042-1047. Eliasen, E., 1958: A Study of the Long Atmospheric Waves on the Basis of Zonal Harmonic Analysis. Tellus, 10, No. 2, 206-215. Eliasen, Erik and Bennert Machenhauer, 1965: A Study of the Fluctuations of the Atmospheric Planetary Flow Pattern Represented by Spherical Harmonics, Tellus, 17, No. 2, 220-238. Eliasen, Erik and Bennert Machenhauer, 1969: On the Observed Large-Scale Atmospheric Wave Motions. Tellus, 21, No. 2, 149-165. Fjortoft, R., 1959: Results Concerning the Distribution and Total Amount of Kinetic Energy in the Atmosphere as a Function of External Heat Sources and Ground Friction. Rossby Memorial Volume, Rockefeller Press, 194-211. Green, J.S.A., 1960: A Problem in Baroclinic Stability. Quart. J. Roy. Meteor. Soc., 86, No. 368, 237-251. Haurwitz, B., 1940: The Motion of Atmospheric Disturbances on the Spherical Earth. J. Mar. Res., 3, No. 3, 254-267. Haurwitz, B., 1941: Dynamic Meteorology, McGraw Hill, New York, 365 pp. Hirota, I., 1968: On the Dynamics of Long and Ultra-Long Waves in a Baroclinic Zonal Current. J. Meteor. Soc. Japan, 46, No. 3, 234-249. Howard, L. N., 1961: Note on a Paper of John W. Miles. J. Fluid Mech., 10, No. 4, 509-512. Jahnke, Eugen and Fritz Emde, 1909: Funktionentafeln mit Formeln und Kurven, B. G. Teubner, Leipzig, 176 pp. Kamke, E., 1961: Differentialgleichungen, Losungsmethoden und Losungen, I, Akademische, Verlagagesellschafg, Leipzig, 666 pp. Kuo, H.-L., 1952: Three-Dimensional Disturbances in a Baroclinic Zonal Current. J. Meteor., %, No. 4, 260-278.

152 BIBLIOGRAPHY (Continued) Lorenz, E. N., 1955: Available Potential Energy and the Maintenance of the General Circulation. Tellus, 7, No. 2, 157-167. Markushevich, A. I., 1965: Theory of Functions of a Complex Variable, II, Prentice-Hall, Englewood Cliffs, N. J., 333 PP. Martin, D. E., 1958: An Investigation of Systematic Errors in the Barotropic Forecasts. Tellus, 10, No. 4, 451-465. Miles, J. W., 1964: Baroclinic Instability of the Zonal Wind. Rev. Geophys., 2, No. 1, 155-176. Miles, J. W., 1965: Instability of Very Long Waves in a Zonal Flow. Tellus, _7, No. 3, 302-305. Murakami, T., 1962: On the Energetics of the Various Large-Scale Atmospheric Disturbances. Papers in Meteor. and Geophys., 8, 103-130. Peixoto, Jose P., 1960: Hemispheric Temperature Conditions During the Year 1950: Scientific Report No. 4, Planetary Circulations Project, Mass. Inst. of Technology, 211 pp. Phillips, N. A., 1963: Geostrophic Motion. Rev. Geophys., 1, No. 2, 123-176. Rainville, Earl D., 1960: Special Functions, The Macmillan Company, New York, 365 pp. Rossby, C.-G., and collaborators, 1939: Relation Between Variations in the Intensity of the Zonal Circulation of the Atmosphere and the Displacements of the Semi-Permanent Centers of Action. J. Marine Res., 2, No. 1, 38-55. Saltzman, B., 1959: On the Maintenance of the Large-Scale Quasi-Permanent Disturbances in the Atmosphere. Tellus, 11, No. 4, 425-431. Saltzman, B., 1965: On the Theory of the Winter-Average Perturbations in the Troposphere and Stratosphere. Mon. Wea. Rev., 93, No. 4, 195-211. Slater, Lucy Joan, 1966: Generalized Hypergeometric Functions, University Press, Cambridge, 273 pp. Smagorinsky, J., 1953: The Dynamical Influence of Large-Scale Heat Sources and Sinks on the Quasi-Stationary Mean Motions of the Atmosphere. Quart. J. Roy. Meteor. Soc., 79, No. 341, 342-266.

153 BIBLIOGRAPHY (Concluded) Welander, P., 1961: Theory of Very Long Waves in a Zonal Atmospheric Flow. Tellus, 13, No. 2, 140-155. Whittaker, E. T. and G. N. Watson, 1961: A Course of Modern Analysis, University Press, Cambridge, 608 pp. Wiin-Nielsen, A., 1961: A Preliminary Study of the Dynamics of Transient, Planetary Waves in the Atmosphere. Tellus, 13, No. 3, 320-333. Wiin-Nielsen, A., A. Vernekar and C. H. Yang, 1967: On the Development of Baroclinic Waves Influenced by Friction and Heating. Pure and Appl. Geophys., 68, No. 3, 131-161.

UNIVERSITY OF MICHIGAN 3 901111 111111 1111115 02826 5786 3 9015 02826 5786