'11 1',; lU N 1. V I': I. 1'' Y 0 I,' M I (C I I (1 \A. N C OLLE'GE OF' ENGIN E'PI1'R ING Department of Atmospheric and Ocennic Science Technical Report ON THE- SIMULATION OF THE AXISYMMETRIC CIRCULATION OF THE ATMOSPHERE Humberto A. Fuenzalida Aksel A.. Wiin-Nielsen Project Director DRDA Project 00260. supported by: NATIONAL SCIENCE FOUNDATION GRANT NO. GA-16166 WASHINGTON, D.C. administered through: DIVISION OF RESEARCH DEVELOPMENT AND ADMINISTRATION ANN ARBOR August 1973

.; K t'It "IV

TABLE[ OF CONTENT'S Page LIST OF TABLES v LIST OF FIGURES vi LIST OF APPENDICES ix LIST OF SYMBOLS x ABSTRACT xv CHA PTER I. IINTRODUCTION 1 1. 1 Motivation of the Study 1 1. 2 Brief Review of Most Recent Previous Work 3 1. 3 A Scope and Outline of the Study 5 II. THE MODEL 8 2. 1 Basic Equations 8 2, 1. 1 Quasi-geostrophic Equations 8 2. 1. 2 Two-level Approximations 15 2. 1. 3 Zonally Averaged Equations 18 2. 1. 4 Subgrid Parameterizations 21 2. 1. 5 Additional Diagnostic Relations 2. 2 Boundary Contitions 25 2. 2. 1 Boundary Conditions at the Poles 26 2. 2. 2 Boundary Conditions at the Equator 32 2. 3 Energetics of the Model 36 2. 3. 1 Momentum Transport 36 2. 3. 2 Heat Transport 42 2. 3. 3 Energy Balance 42 III. DIABATIC HEATING 49 3.1 Introduction 49 3. 2 Heat Balance for an Atmospheric Column 50 3. 3 Surface Heat Balance 53 3. 4 Vertical Transfer Processes and the Water Balance 55 3. 5 Computation of the Heating Function 58 3. 6 Tuning of the Heating Model 60 3. 6. 1 Numerical Procedure 61 3. 6. 2 Determination of Parameters 62 3. 6. 3 Characteristics of the Tuned Heating Function 68 3. 6. 4 Diabatic Heating in the Simulation 78 iii

TABLE OF CONTENTS (continued) Page IV. LINEAR SIMULATIONS 83 4. 1 Northern Hemisphere Simulation 85 4. 1. 1 Temperature Field 85 4. 1. 2 Zonal Motion 88 4. 1. 3 Meridional Circulation 90 4. 1. 4 Transport of Heat and Momentum 92 4.1. 5 Energetics 92 4. 2 Southern Hemisphere Simulation 100 4. 2. 1 Temperature Field 100 4. 2. 2 Zonal Motion 103 4. 2. 3 Meridional Circulations 105 4. 2. 4 Transport of Heat and Momentum 105 4. 2. 5 Energetics 105 V. SENSITIVITY STUDY 111 5. 1 Streamfunction Constant 111 5.2 Static Stability 114 5. 3 Boundary Layer Friction 5.4 Internal Friction 117 5. 5 Eddy Coefficients 118 5. 5. 1 Heat Eddy Exchange Coefficient 120 5. 5. 2 Potential Vorticity Eddy Exchange Coefficients 124 5.6 An Improved Case 127 VI. EDDY EXCHANGE PROCESSES 133 6. 1 Nature of the Eddy Exchange 6. 1. 1 Exchange Properties of Eddies 136 6. 1. 2 Generation of Eddies: Baroclinic Instability 140 6. 2 Recent Models of the Heat Transport 143 6. 3 A Parameterization for Momentum Transport 145 6. 4 Results of a Nonlinear Simulation 148 VII. CONCLUSIONS 160 APPENDIX I 165 APPENDIX II 173 APPENDIX III 176 A'PPENDIX IV 182 REFERENCES 184 iv

LIST OF TABLES Table Page 3. 6. 1 Toulndary layer vertical transfer of water vapor and sensible heat 63 3. 6. 2 Latent heat release normalization factor:m 65 3. 6. 3 Radiation parameters after tuning 69 3. 6. 4 Energy balance after tuning 71 5. 0. 1 List of sensitivity experiments 112 5. 5. 1 Eddy exchange coefficients from Sela & Wiin-Nielsen (1971) 119 v

LIST OF FIGURES Figure Page 3. 2. 1 Energy fluxes related to the diabatic heating of the atmosphere 52 3. 3. 1 Energy fluxes related to the surface layer 52 3. 6. 1 Energy balance of the earth-atmosphere system 72 3. 6. 2 Heat transport by the ocean-atmosphere system and by the oceans 73 3. 6. 3 Heat released by the oceans 75 3. 6. 4 Mean annual heating function 76 3. 6. 5 Atmospheric heat budget and its components 77 3. 6. 6 Seasonal distribution of the atmospheric diabatic heating 79 3. 6. 7 Diabatic heating of an atmospheric column 80 3. 6. 8 Comparison of mean annual heating obtained in the simulation and the tuning 81 4. 0. 1 Flow chart of computation of linear simulations 84 4. 1.1 Thermal field at 50 cb for the Northern Hemisphere in the linear simulation 86 4. 1.2 Thermal field at 100 cb for the Northern Hemisphere in the linear simulation 86 4. 1. 3 Meridional temperature profile for January and July at 50 cb, Northern Hemisphere linear simulation 87 4. 1.4 Meridional temperature profile for January and July at 100 cb, Northern Hemisphere linear simulation 87 4. 1.5 Zonal flow at 25 cb for the Northern Hemisphere in the linear simulation 89 4. 1. 6 Zonal flow at 75 cb for the Northern Hemisphere in the linear simulation 89 4. 1. 7 Vertical velocity at 50 cb for the Northern -Temisphei - re in the linear simulation 91 4. 1. 8 Heat transport by the atmosphere for the Northern HTemisphere in the linear simulation 91 vi

IST ()F FIGURES (continued) Figure Page 4. 1. 9 Angular momentum transport for the Northern Hemisphere January in the linear simulation 93 4. 1.10 Simulated variations of zonal available potential energy and zonal kinetic energy for the Northern Hemisphere in the linear simulation 95 4. 1.11 Simulated variations of the generation of zonal available potential energy and dissipation of zonal kinetic energy for the Northern Hemisphere in the linear simulation 95 4.1.12 Simulated variations of energy conversions for the Northern Hemisphere in the linear simulation 98 4. 1.13 Simulated energy box diagram in a typical year for the Northern Hemisphere in the linear simulation 98 4.2. 1 Same as Figure 4. 1. 1 but for Southern Hemisphere 101 4.2.2 Same as Figure 4. 1.2 but for Southern Hemisphere 101 4. 2. 3 Same as Figure 4. 1. 3 but for Southern Hemisphere 102 4.2.4 Same as Figure 4. 1.4 but for Southern Hemisphere 102 4.2. 5 Same as Figure 4. 1.5 but for Southern Hemisphere 104 4.2. 6 Same as Figure 4. 1. 6 but for Southern Hemisphere 104 4.2. 7 Same as Figure 4. 1. 7 but for Southern Hemisphere 106 4.2. 8 Same as Figure 4. 1.8 but for Southern Hemisphere 107 4.2. 9 Same as Figure 4. 1. 9 but for Southern Hemisphere and showing corrected values for January and July 107 4.2. 10 Same as Figure 4. 1. 10 but for Southern Hemisphere 109 4.2.11 Same as Figure 4. 1.11 but for Southern Hemisphere 109 4.2. 12 Same as Figure 4. 1. 12 but for Southern Hemisphere 110 4.2. 13 Same as Figure 4. 1.13 but for Southern Hemisphere 110 5. 5.1 Profiles of K2 in experiments #7, #8 and #9 122 vii

LIST OF FIGURES (continued) Figure Page 5. 5. 2 Temperature profile for January at 50 cb in experiments #10, #11 and #12 122 5. 6.1 Same as Figure 4.1.5 for an improved case 129 5. 6.2 Same as Figure 4.1.6 for an improved case 129 5. 6. 3' Same as Figure 4. 1. 7 for an improved case 130 5. 6. 4 Same as Figure 4.1.8 for an improved case 130 5. 6. 5 Same as Figure 4.1.13 for an improved case 1 32 6. 4. 1 Seasonal variations of eddy coefficients in the nonlinear simulation of the Northern Hemisphere 149 6.4.2 Same as Figure 4. 1.1 for nonlinear simulation 151 6.4. 3 Same as Figure 4.1.2 for nonlinear simulation 151 6.4.4 Same as Figure 4. 1. 3 for nonlinear simulation 152 6. 4. 5 Same as Figure 4.1.4 for nonlinear simulation 152 6.4. 6 Same as Figure 4.1. 5 for nonlinear simulation 153 6.4. 7 Same as Figure 4. 1. 6 for nonlinear simulation 153 6. 4.8 Same as Figure 4.1. 7 for nonlinear simulation 155 6.4. 9 Same as Figure 4. 1.8 for nonlinear simulation 156 6. 4.10 Same as Figure 4. 2. 9 for nonlinear simulation 156 6. 4.. 11 Same as Figure 4. 1.10 for nonlinear simulation 158 6.4. 12 Same as Figure 4. 1. 11 for nonlinear simulation 158 6. 4. 13 Same as Figure 4. 1. 12 for nonlinear simulation 159 6.4. 14 Same as Figure 4. 1.13 for nonlinear simulation 159 vi. i. i

LIST OF SYMBOLS A internal friction factor, definedon page 22 Az zonal available potential energy a radius of the earth b turbulent transfer of sensible heat across the boundary layer C(MN) conversion of energy from form M to form N c complex phase velocity ( = cr + i c) r i c specific heat at constant pressure for dry air cd drag coefficient E evaporation f Coriolis parameter ( = 2 l sin g ) f standard value of Coriolis parameter f(, t) shape factor defined on page 182'iT ~horizontal friction stress vector at level i $M., ~ horizontal friction stress x component at level i 1, x G generation of available potential energy or a quasi-conservative quanti ty g acceleration of gravity H diabatic heating H( downward flux of energy by process i related to the diabatic heating of an atmospheric column H(i) downward flux of energy by process i across the surface HT poleward heat transport by the atmosphere (HT)2 poleward heat transport by the atmosphere at level 2 HI, scale height I integral defined in page 56 x

Ptage APPENDIX I:.FINITE DIFFERENCES EQUATIONS 165 APPENDIX II: SPLINE INTERPOLATION 173 APPENDIX III: INSTABILITY OF A TWO-LEVEL MODEL 1 76 APPENDIX IV: TIME DEPENDENT COEFFICIENTS 182 ix

LIST OF SYMBOLS (continued) ji fraction of a latitude circle covered with ith kind of surface K kinetic energy Ki eddy exchange coefficient at level i k coefficient defined on page 56 k, constant used in expression (3.4.2) k vertical unit vector L typical horizontal scale or latent heat for water M total angular momentum transport Mi angular momentum transport at level i n parameter of beta distribution P (x) Legendre polynomial of order n p independent variable pressure Qji potential vorticity at level i q constant defined on page 16 R gas constant for dry air Ro Rossby number defined on page 1 Ri Richardson number defined on page 12 Ro incoming solar radiation at the top of the atmosphere r rainfall or parameter of beta distribution defined on page 182 r atmospheric albedo a r surface albedo Ti temperature at level i TD oceanic subsurface temperature t independent variable time U typical wind speed or velocity of a basic flow xi

LIST OF SYMBOLS (continued) u zonal velocity V/ wind vector V component of V v meridional velocity x zonal coordinate or independent variable (= sing ) y meridional coordinate or dependent variable z height of an isobaric surface a specific volume of the dry air P Rossby parameter or beta probability density function r atmospheric absorption to longwave radiation or spectral density function y ratio of specific heats for dry air 8 aspect ratio or factor defined on page 45 E ratio of 100 cb temperature to sea surface temperature ~ boundary layer friction constant defined on page 21 or dimensionless number defined on page 12 Ef ~relative vorticity r~) meridional displacement of a parcel 0 potential temperature K typical magnitude of the static stability of the atmosphere ( =ln 0/Iz) X independent variable longitude or constant defined on page 16 An nth eigenvalue 4 netmsoscalet dyrla.lmic viscosity aV mesoscale kinematic viscosity ( = / / p ) ^J 1d

LIST OF SYMBOLS (continued) 91 atmospheric downward emissivity f2 iatmospheric upward emissivity P density of dry air ur velocity potential or Stefan-Boltzmann constant 0" ~static stability T time scale perturbation streamfunction ) anamplitude of P) or geopotential of an isobaric surface? independent variable latitude X atmospheric opacity to solar radiation qi streamfunction fl angular velocity of the earth W vertical velocity in the isobaric system X ( =R/cp) ( ) irrotational component ( ), solenoidal component ( )A zonal component ( )ej meridional component ( )r radial component ( )s standard or surface value s ( ) constant value ( )p isobaric component or particular solution ( )H horizontal component )H ( )0 variable at 0 cb level xiii

LIST OF SYMBOLS (continued) )I( variable at 25 cb level )2 variable at 50 cb level )( ) variable at 75 cb level ( )4 variable at 100 cb level ( )z zonal mean ( ) eddy component, deviation from zonal mean ( )' isobaric deviation ( ) time mean or horizontal mean ( )* complex conjugate T ( ) transpose matrix xi.v

A I STJ:ACT ( r ON THE SIMULATTON OF TT-TE AXISYMMErTRIC CIRCULA7TION OF THE ATMOSPHERE by Humberto Adrian Fuenzalida Chairman: Aksel C. Wiin-Nielsen The governing equations of a two-layer quasi-geostrophic model are cast in terms of potential vorticities and zonally averaged. Using suitable boundary conditions at the pole and equator they are integrated in time for three years. The thermal forcing is provided by an ad hoc diabatic heating function which includes radiative fluxes, latent heat released by condensation, turbulent transport in the boundary layer and thermal interaction with the oceans. The poleward eddy fluxes of sensible heat and angular momentum are parameterized first by means of diagnostically determined exchange coefficients for heat and potential vorticity, which depend on latitude and pressure, but do not vary in time. The integration results include the temperature field at 100 and 50 cb, the zonal flow at 75 and 25 cb, and the vertical motion at 50 cb, all of which are evaluated every third day. The kinetic and available potential energies and their conversions, generation and dissipation are also evaluated with the same frequency. Starting from isothermal and motionless conditions the model reaches a quasi-steady state after a few months. Separate integrations are performed for each hemisphere using the same set of eddy coefficients determined from Northern Hemisphere data by Wiin-Nielsen and Sela (1971). By comparison with empirical data the outcome of the simulation shows that substantial improvement is obtained with respect to previous results xv

(Sela and Wiin-Nielsen, 1971) by virtue of a more realistic thermal forcing. A sensitivity analysis of the model reveals that some results can be modified by a change of the value of the streamfunction constant and the boundary layer friction and that assuming an exchange coefficient independent of latitude and pressure leads to unrealistic results. An analysis of the mixing process for quasi-conservative properties together with results by Green(l 970), Bretherton(l 966), and Pedlosky(1 964), and the hypothesis that unstable waves are the responsible agents for the eddy transport, leads to a parameterization of the eddy fluxes such that they depend on the mean gradient of potential vorticity and vary with time. This is a generalization of a result by Stone(1972). An integration of the model using this parameterization provides results whi"ch are in better agreement with observations than those obtained in the former simulations. xvi.

CHAPTER I INTRODUCTION 1.1 MOTIVATION OF THE STUDY Today there exists a complete hierarchy of atmospheric models. In terms of spatial dimensionality they range from the Global Average Models (GAM), which due to a very detailed vertical specification allow a description of photochemical and radiative processes, to the so-called General Circulation Models (GCM) which cover the globe with typically ten levels of a 300-km mesh size grid. Between the onedimensional GAM and the three-dimensional GCM stand the Models in the Meridional Plane (MPM) with a two-dimensional grid and a moderate cost of operation. They were not created as an economic compromise, but rather are the natural outcome of the symmetry imparted to the atmospheric motions by the rotation of the earth and the thermal forcing imposed by the sun. Atmospheric motions can be decomposed in a zonal (longitudinal ) mean and deviations from it. In extratropical latitudes and for motions in the meridional plane, the temporal mean of the deviations over, say a month, is orders of magnitude smaller than the instantaneous values revealing a turbulent behavior. This turbulent motion is able to diffuse any entity, in particular those with Lagrangian variations smaller than local changes, i. e., a quantity which is conservative to some extent. Due to the nonlinearity of the governing equations the isolated study of the zonal component of the atmospheric dependent variables is not possible. The zonally averaged equations of 1

2 momentum and energy will contain mean products of the variables, whose behavior cannot be predicted by the system formed by these equations only. Thus, starting with a complete set of relations, we end up with more unknowns than equations. This closure problem is the basic difficulty faced in the formulation of the MPM. It is not present in the GAM because the nonlinear terms, which represents interaction within the fluid, vanish when integrated over a closed surface. Neither is it present in the GCM which deals explicitly with the zonal deviations. The MPM have not enjoyed much popularity among scientists dedicated to climate study probably because the parameterizations of heat and momentum transports have never been done properly. In the present state of climatic simulation they could be extremely helpful since the GCM have a prohibitive cost in experiments comprising several years and the GAM are unable to consider interactions between low and high latitudes. Lately, the interest in the MPM has been revived by the immediate need of answers to questions related to the impact of present and future technology on the environment. Such answers are to be obtained through climate simulations covering several decades if the interactions between the oceans and the atmosphere are to be included. Even with the next generation of computers this is a prohibitive task for GCM, (SMIC, 1971). Therefore, it is not a mere academic exercise to attempt new formulations of the eddy transfers in a MPM. On the contrary their applications are presently needed.

3 From a less applied point of view, the ability to simulate the zonal part of the atmospheric phenomena provides a frame within which the deviations or eddies must evolve and possesses a substantial intrinsic interest since the zonal motion constitutes a large fraction of the total atmospheric energy content. 1.2 BRIEF REVIEW OF MOST RECENT PREVIOUS WORK The kind of turbulent motion we are attempting to parameterize differs in several important aspects from the most commonly studied turbulence. First of all it is highly anisotropic. Typical eddies range from 10, 000 to 3, 000 km in horizontal dimensions while their vertical extent hardly exceeds 10 km. This configuration forces the particles trajectory to be quasihorizontal. Secondly, the turbulence is nonhomogeneous. Mechanisms that generate this turbulent motion (baroclinic instability, large scale orography and heat sources) are active over most of the wave number range so that eddy motions generate mean values (and other statistics) with important variations in latitude and longitude. Finally, because of the seasonal variations of the acting mechanisms the process in nonstationary. Therefore, the extension of results of turbulence theory to our case is not permissible unless specifically justified. This parameterization problem is an old one. Lately the heat transport has been modelled successfully (Stone, 1972b), but the momentum transport is still waiting for a solution. This may be due to the fact that momentum as opposed to temperature is not a quasiconservative quantity (Green, 1970).

4 Following the completion of a diagnostic study of the annual variations of the atmospheric energy cycle, Wiin-Nielsen (1967), this author attempted to simulate the part related to the zonal flow using a simple dynamic model forced by a Newtonian heating function WiinNielsen, 1969, 1970). In these investigations the heat transport was parameterized through an austauch coefficient and the momentum transport was neglected in a first order approximation. Taking Green's (1970) suggestion, that exchange coefficients should be considered in relation to conservative quantities only, WiinNielsen and Sela (1971) evaluated annual diagnostic coefficients for the exchange of quasigeostrophic potential vorticity at several levels and every 21 degrees of latitude. Casting the two-level quasigeostrophic model in terms of potential vorticities Sela and Wiin-Nielsen (1971) made another attempt to simulate the axisymmetric component of the general circulation. The heat transport was parameterized by an eddy coefficient dependent on latitude and the thermal forcing was Newtonian. The results of this investigation gave the impression that much could be improved by a more realistic specification of the heating function. In particular, the release of latent heat by condensation and the exchange of sensible heat with the underlying surface of the earth were considered to be the most important mechanisms ignored in the simulation. The eddy exchange coefficients evaluated by Wiin-Nielsen and Sela (loc. cit. ) constitute a first step to prove the usefulness of parameterizations of conservative quantities and are able to generate good simulations of the annual energy cycle, but even if such evaluations were developed to the extent of providing a time dependence for the coefficients these diagnostic values are closely linked to the present

5 climatic conditions and are of little use in predicting climatic evolutions under conditions significantly different from the current ones. One has to realize that parameterizations,because of their very nature are approximations to the actual mechanisms. They can never attain complete fidelity and their success depends on the fraction of the variance of the process that can be accounted for in terms of dependent variables. Ideally a parameterization must depend to a large extent on the internal variables of the model and include coefficients or functions with a much smaller variability than that of the process to be modelled. Both Green (1970) and Stone (1972b) have designed parameterizations for the heat transport using an eddy exchange coefficient proportional to the meridional gradient of temperature. The proportionality factor depends on quantities which have a more restricted variation than the heat flux itself. We will attempt a time dependent formulation which can be considered an extension of their results and because it introduces a feedback mechanism it is more useful in the study of climatic variations. 1. 3 A SCOPE AND OUTLINE OF THE STUDY This study endeavours first to show that a crude dynamic model can account for most of the energy cycle and reproduce its annual variations to a larger extent than that obtained by Sela and Wiin-Nielsen (1971). At the same time these results provide more evidence to support the parameterizations in terms of quasi-conservative quantities. A second related goal is an improvement in the modelling of the momentum transport.

6 To attain the first goal the model used by Sela and WiinNielsen (loc. cit. ) was critically reviewed, its boundary conditions were studied from an analytical point of view and the energetics formulated without referring to balance considerations. All this is done in Chapter 2. Chapter 3 contains the description of the thermal forcing which now includes latent heat release and exchange of sensible heat with the oceans. This formulation also allows the evaluation of surface temperatures. The heating function is tuned to satisfy the present seasonal variations and its behavior in the simulations is quite satisfactory. Results from the simulations of both hemispheres separately are reported in Chapter 4. These results are directly comparable with those of Sela and Wiin-Nielsen (loc. cit. ) since values of constants and the eddy exchange coefficients are the same. The only difference is a larger time step and grid size. The validation method used to judge the performance of the model is a comparison with observed values. At this point we hope to have shown that substantial improvement was obtained by the use of a more realistic thermal forcing. The main deficiencies of the model are due to its quasigeostrophic character and the time independent value assigned to the eddy exchange coefficients. It is felt that the advantage of parameterizing the transport of a quasi-conservative quantity such as potential vorticity rather than momentum, which has a very active source in the Coriolis force, is confirmed by this result.

7 Chapter 5 contains the outcome of sensitivity studies. Various constants of the model are modified to find to what extent the results depend on their particular values. A numerical study is made of the most relevant features of the latitudinal dependence of the three eddy exchange coefficients included in the model to extimate their importance and identify features that depend on them. One result that follows from this study is the inadequacy of a constant austauch coefficient modelling for the transport under consideration, Finally, the results of an improved case are shown. In Chapter 6 an analysis is made of the different mechanisms that can perform the heat and momentum transport. Under the assumption that the main agents responsible for these fluxes are the waves generated by the baroclinic instability of the zonal flow some results obtained by Green (loc. cit. ) and Stone (loc. cit. ) are generalized in the light of investigations by Bretherton (1966b) and Pedlosky (1964a) to produce a time dependent formulation for the exchange coefficients. Finally, Chapter 7 summarizes the conclusions of the work and contains some ideas as to how to make improvements in further attempts.

CHAPTER II THE MODEL 2.1 BASIC EQUATIONS This study will be done by using the standard two-level quasi-geostrophic model introduced by Phillips (1956). In the next pages we will summarize the most important simplifications introduced in its formulation and indicate the limitations they convey. 2. 1. 1 Quasi-geostrophic Equations The dynamic model to be used is based on the quasigeostrophic equation of vorticity. (l-at +/V)+ f -oir pkV(2.1. 1) where the independent variables are time, t, and the space coordinates AN,, Ap, longitude, latitude and pressure, respectively. Hence o the local derivative a and the two-dimensional operator \V V V+ a ~ aV os aA a C a (2.1.2) are evaluated isobarically. The horizontal velocity has been decomposed in a solenoidal part VkP V'VY, where [ is the vertically upward unit vector and V a streamfunction, and an irrotational component V =Vc, where G is a velocity potential. The vertical component of the relative vorticity is given by the horizontal part of the Laplacian of the streamfunction. 8

i- = -Vo b [ ^^ ~)1? -QoL a^ ^ cow l aA ^^ an aC aCosop? ACOth If ^ L +a__ (2.1.3) As customary C represents the Coriolis parameter, 2 Q Arvi, where 2 is the angular velocity of the earth. The subscript zero denotes a constant value. The divergence operator is given by ~VW^ sa lI a a[V ~I (Vcoi')] (2.1.4) By the same token the vertical component of the rotational operator is AgtJ - c0I a1 A -'e a) where the vector V has components (V%, V.,Vr ) along unit vectors oriented west to east, south to north and upwards respectively. The last term in (2. 1.1) accounts for the frictional force due to the vertical variation of the horizontal stress'. Here represents the acceleration of gravity. The irrotational part of the wind will be eliminated using the equation of continuity V'%+ ab=p (2. 1.5) where 00 is the individual variation of pressure, the vertical "velocity" of the coordinate system x, y, p. A third relation to be used is the first principle of the thermodynamics in its quasi-geostrophic form

10 Vat.apv - ap - V c, (2 (2. 1. 6) Here the subscript s denotes a standard value independent of the horizontal coordinates and time. The potential temperature is related to the temperature T, density e and pressure by 6= wrt= ot' p- (2. 1.7) G =Tp ~con d C~rrst.(;"F (2.1.7) where X =C-, R is the dry air gas constant and Cp the specific heat of the dry air at constant pressure. Finally z is a dependent variable which denotes the height of an isobaric surface and H the diabatic heating rate. In the quasi-geostrophic formulation the balance equation, which completes the basic set, is reduced to "V = f (, where ( = gz. In the process of obtaining this set of relations starting from the complete equations (vorticity, divergence, continuity and the first principle of thermodynamics) scale, filtering and constraintpreserving simplifications that limit the applicability of the equations have been introduced. A brief description of them,based on the discussion of Charney (1960) and Charney and Stern (1962),follows. The first filtering approximation we will be concerned with is the hydrostatic assumption which filters out vertically travelling sound waves. The scaling condition equivalent to the replacement of the vertical equation of motion by the hydrostatic relation is d =-d <( where D) and L are vertical and horizontal distances across which the dependent variables vary in their typical order of magnitude.

11 S is sometimes called the aspect ratio. This approximation upsets the conservations of angular momentum and energy under adiabatic and reversible conditions. In order to restore them the kinetic energy has to be redefined excluding the vertical motion,Coriolis and metric terms including the vertical velocity must be deleted from the horizontal equations of motion, and the radial spherical coordinate must be replaced by the radius of the earth ("shallowness" approximation). The hydrostatic approximation is satisfied with a high degree of accuracy in the atmospheric large scale systems and permits the introduction of the isobaric quasi-Lagrangian coordinate system x, y, p. The next filtering approximations are the geostrophic conditions. These are stated as a small Rossby number,R-= 2,0() where U is a typical value for the horizontal velocity, and by setting a lower bound to the value of the typical time scale, T, (time required for a dependent variable to change by its characteristic magnitude) in terms of the advective time: T a 0 (LO1J). Although the Rossby number just defined does not show it explicitly this approximation fails in the vicinity of the equator. The small Rossby number allows a perturbation expansion which has as zero order approximation the geostrophic relation, and as first order approximation the quasi-geostrophic equations. The replacement of the divergence equation by the geostrophic relation filters out gravity waves and to maintain the energy constraint it is necessary to redefine the kinetic energy as that due to the solenoidal part of the wind only. The third scale consideration relates to the importance of the conversion of available potential into kinetic energy. The term

12 that determinates the conversion in the vorticity equation is the divergence term whose order of magnitude is then assumed to be one or more, e. g. = 40 L > ^ (om ~= sD -i where K is the typical magnitude of the static stability, aLt i/., and LR is the Rossby radius of deformation. Since SRo2 Ri = 1 this condition is sometimes stated as requiring a large value of the Richardson number UiTh The next set of relations are chosen as to guarantee horizontal pressure deviations much smaller than the magnitude of the variable itself. These are )K<< I, D/H, 0(0). Since g0(0) the first relation implies a very small value for the static stability. In fact observations show ~ to be order one, and Eady (1949) concluded that this corresponds to the most efficient systems in converting potential to kinetic energy. In this case )(< 0() is enough. The second relation selects systems whose vertical extension D is similar or smaller than the scale height H*. The smallness of the horizontal baric deviations carries with it a similar bound on the density, specific volume and more important, on the static stability, which then can be replaced by a standard value within the accuracy in use (first order in R ). However, such a substitution breaks the link between total potential and kinetic energy and no stabilization results from the energy conversion, Lorenz (1960). This assumption does not seem of great importance in climate simulation according to some recent results by Stone (1973).

13 The relation 1)/11, ()(1) together with the hydrostatic and geostrophic approximations imiply an irrotational part of the wind, V. of order FRo while the solenoidal part, Vy, is order one. By the same token the relations used to transform derivatives from a quasiLagragian system of reference (x, y, 5 ), isobaric or isentropic, to the fixed (x, y, z) system are of the form ( ) (-a1 - (E(R) where S stands for time or any horizontal space variable, and all variables are dimensionless. Finally, we use a simplification of the governing equations which can be justified through a scale analysis when the beta-plane geometry is being used. This simplification is related to the space variations of the Coriolis parameter in the divergence term of the vorticity equation and in the linear balance equation. The spherical geometry makes invalid the scale argument normally used. But, the simplicity introduced by this approximation is most desirable since it makes the streamfunction proportional to the geopotential. Because this approximation gives good results in the beta-plane geometry and the great simplification it provides we will use it. Once this decision has been taken it is necessary to replace the Coriolis parameter by a constant value in the terms affected by the previous simplification. Such substitution has to be done to avoid fictitious sources of vorticity. However, this is obtained at expense of a discontinuity in the Coriolis parameter at the equator and henceforth on the geostrophic winds. The sole simplification of the governing equations by an order of magnitude criterion can lead to serious inconsistencies when

14 they are integrated over long times. One way to keep some control on this aspect is to modify the equations without altering energy invariants that the original complete set of equations satisfy. Although in some of the filtering approximations care was taken to maintain the energy invariants, at least with good approximation, it is necessary to check the final equations and the energetic processes derived from them. We will do this in 2. 3.3. Equation (2. 1.1) and (2. 1. 5) allow the elimination of V,7 giving (%at +8*) +i) oap-i aid (2. 1. 8) where the last term represents the friction term that includes internal and boundary layer effects, is the horizontal stress vector. Substituting ( = gz and T=-, a standard value depending only on p,, and including the heating term we can write (at +W.v) a H- R 4 P'(2.1.9) Equations (2. 1. 8) and (2. 1. 9) constitute the basic relations which are used in the two-level quasi-geostrophic model, beside XF'. Eliminating the geopotential and substituting the streamfunction where necessary we obtain (at 1~kV 1 V)(Vz v +F) =o a- 2 H (8t +kX5;7YV)8tv, C+c, a, @ =_ At(2.1. 11)

15 These equations determine the solenoidal field of horizontal motion throughY and the vertical circulation by u, when the stress ] has been expressed as a function of them and the diabatic heating H is prescribed. In frictionless adiabatic conditions these equations (by elimination of w) define an invariant under horizontal nondivergent motion (at +v/ v)(FI + a4 (2. 1. 12) this is the pseudo-potential vorticity of Charney (1971), which replaces Ertel's potential vorticity in quasi-geostrophic motion (Charney and Stern, 1962). 2. 1. 2 Two-level Approximations In the two-level model the vertical resolution is kept to a minimum: two levels for s1 and one for w. Using the classical notation we will denote the levels 0, 25, 50, 75 and 100 cb by 0, 1, 2, 3, 4 respectively. 0 - u= 0 O0 cb 1..... r1..25 2. 2 50 3. —- 3 75 4 ~w 4=0 - - 100 the boundary conditions are then applied to w: w = 0 at levels 0 and 4. Orography effects are not considered in this study.

16 Equation (2. 1. 10) applied to levels 1 and 3 gives aa; W' (F + Z) = 210() \Z0/~( @t_2Z-(2 1 13) at 2]s/ (2. 1.14) where )P= 100 cb and,'j/3 are understood to be solenoidal velocities. Equation (2. 1. 11) applied to level 2 gives at (t.-',)+ ~, (,) ~ -F H, To close the system we set v i(v+ Substituting L2 from (2. 1. 15) into (2. 1. 13) we obtain (4I o)[F+C. - ~-v.(~-~o a t T; H- 2 We define a thermal wind by rT= (V,.-) and a corresponding streamfunction jT {(M -V)' Then replacing we 8can wre we can write (at~ +t(F +, -t)=

17 In a similar way we get (2.1.17) These two equations define the potential vorticities at levels 1 and 3. These quantities are going to be of central importance in the eddy flux param eterization. Note that the left hand side of these equations can be expressed in terms of ~, andW only, since`f, = vlt:V2 1P X! =V2V3 3 kV (x) Defining the potential vorticities "Q,-f~~ +vV,~ A^ ^(2.1.18) Q= F + vw+' (2. 1. 19) we can formulate a diagnostic identity as follows vZ tP v-=At(Q1 - Q6) (2.1.20) and write our prognostic equations in the form: at I I)' cal; - V Q1)-XII - ( 2 - (2.1.21) a t = - V (%Q^)+ H -% Ps i 7> (ffi4 ~ ) (2. 1. 22) ^-^^^ ^'

18 The thermal streamfunction relates to the temperature at level 2 by the relation = Fo, which gives o AT T 2 (' )3 Using the hydrostatic relation applied at level 2 we obtain O = 2 ~= 2, z2 ) where T2 is the mean temperature between 25 and 75 cb. Therefore -2FTr ='T, (2.1.23) 2. 1. 3 Zonally Averaged Equations In order to study the axisymmetric circulation we average the governing equations with respect to longitude. Let us define the operator 27Tr ( ) = i ( )d^( )~- o and apply it to equations (2. 1. 20) to (2. 1. 23). Then aQt = acor [c> (Co 1] 2z (2.. 24)a (2. 1.24)

a - - ~ __o (Q ) ^ CO-x J(2.1.25) W^ y CT7. 2 (2.1. 26) 2 fo TIZ= RT21 (2.1. 27) The advective terms in (2. 1. 24) and (2. 1. 25) can be decomposed as: ( )LL [( L )(QL Q + )] LQl \ z JZ = HiZ Q\ + E( Q,), ( eQ) since in a quasigeostrophic model 0L = O. The subscript E denotes deviations from the zonal mean. Considering simple parameterizations of the internal and boundary layer friction described below and an exchange coefficient law for the eddy transport of potential vorticity, whose validity will be analyzed later (,c Q X -Kc )^ (!)' Q' (2. 1. 28)

20 ( ) ( )a0 (2.1.29) the prognostic equations become at Q a~c: sKC^ aaQ2)- JTb o^T~r^^ - c^^^^^^ ^g (K^ccnD a'l -bH,,-n(2. 1. 30) at = acos (K ad Q,7; Z + ( bs) (..Q — ) (2.1, 31) These equations allow the computation of future potential vorticities QMz and Q37 if we know the current values of the heating N2z and the vorticities EL,, Qi, Q3z. With the new values of Q&. andQz equation (2. 1. 26) permits the evaluation of the new \. Then (2. 1. 27) gives the temperature 22 and the identities Vt -: - ^ )T;2 (Q 5) t ^ 7,3 ) (2. 1. 32) give the new vlaues of B, and fz. Finally if we have Z as a function of Tz, for all time and location, a new time integration is possible and this process can be repeated indefinitely.

21 Knowing the vorticities we can compute the zonal flow by integration. The vertical motion can be calculated from the thermodynamic relation (2. 1. 15) or the omega equation of the system. Meridional flow, if necessary, can be obtained by the continuity equation (2. 1. 5). The temperatures at the 100 cb level will also be obtained through the heating function H2z. 2. 1. 4 Subgrid Parameterizations The only subgrid processes that will be included, outside of those considered in the heating function, are the boundary layer friction and the internal dissipation due to vertical shear. Boundary layer friction enter the vorticity equations (2. 1. 24) and (2. 1. 25) in terms of the form F MOM CcM V -I It will be assumed that'=0 and that 3 is entirely due to the vertical flux of momentum in the Prandtl layer. A simple formulation for this process is I4 =Cd,4 V which with an additional linearization, that consists of using V4 and A as constants, gives fs aCo3J as (Cap E) P- C atp 4 4 The drag coefficient Cd can attain values between 10 -2 and 10 depending on the roughness of the underlying surface (Cressman, 1960). Then define p 2a 4Cd ^ ^ ^~V

22 a constant for a given value of C^. In the experiments a value 6 -1 = 3 x 106 sec was chosen corresponding to C =2.5 for V4 =5msec. The unknown value of vorticity at 100 cb will be estimated by linear extrapolation as ~4 = (- )= —2- T Then 1 b aTc 4 X call Outside the surface layer we will use ~ so that, 2 2' =2, In finite difference form ZX R2T Considering A V- as a constant we obtain z 8 135 acost.a+ (%>@)l=-2An; — A(,z-) Therefore, the dissipation term in (2. 1. 24) becomes'PS acorspodn terxi (2.1.5 is while the corresponding term in (2. 1. 25) is 2__ t [ (8 cos(-Yid A ^M ^ (-^ ) vpS;a 2c~ bi L x txz x Z Vs

23 2. 1. 5 Additional Diagnostic Relations Equations (2. 1. 30), (2.1. 31), (2. 1. 26) and (2. 1. 27) together with identities (2. 1. 18) and (2. 1. 19) provide a way to evaluate the vorticities, 5, the thermal streamfunction.T and the temperature at a midtropospheric level, T-, once a procedure to obtain the diabatic heating is designed. This problem will be postponed until Chapter 3, where it will be shown how to obtain also the temperature at 100 cb or a near surface level. Presently we will be concerned with the determination of the zonal flow and the meridional circulations. Zonal velocities at levels 1 and 3 will be computed by integration from the meridional distribution of relative vorticities. This procedure is more accurate and economic than the use of streamfunctions. By definition?T^ acotQp ben (d.{ ) x2 CaT An integration from equator to pole gives 7r/z ui. (0=o)=-| fCofdY 0 Then integrating from equator to any latitude circle one obtains ULZ(^^) )ccn 9 1z cwc'dvP^ ^(2.1.33)

24 This relation can be used also for LLU since it depends on vorticity only. The vertical motion can be calculated from the thermodynamic equation in the geostrophic formulation (2. 1. 15). In terms of the thermal streamfunction this equation can be written as at 4 T^ {\^ (828r) -, 4 f^2= __Cc (2. 1. 34) Taking the zonal average we obtain -t ~t fa( A ) [Z1- u HO _ 4 (2. 1. 35) Since in a quasi-geogtrophic theory tz = 0 \(AII )i (u2E\)Z K/')= - (2. 1.36) The last equality needs some comments. In Chapter 6 we will try to justify a diffusive parametrization for quasi-conservative properties. The thermal streamfunction does not fit into this category in the same sense as potential vorticities, iL e., considering only the advection by the solenoidal part of the horizontal wind. Equation (2. 1. 34) shows that VT is conserved under adiabatic conditions when the contribution of the vertical motion is also taken into account. However, Clapp (1970) has provided observational evidence that a parameterization through an eddy coefficient is adequate for the horizontal heat transfer by transient eddies. We will use it realizing that some inconsistency in our approach is introduced. Equation (2. 1. 35) allows the evaluation

25 of O, if the eddy coefficient 2(p) is available: U.): ~i a } az z- c -. H c (2. 1. 37) " L a as/ M Once the vertical motion is known the meridional velocity is provided by the equation of continuity integrated with respect to latitude: =-B OR CON sI C (2.1.38) 2.2 BOUNDARY CONDITIONS It is our purpose to integrate equations (2. 1. 30) and (2. 1. 31) in time and relation (2.1. 26) in space over each hemisphere. To accomplish such a task we need boundary conditions at both poles and equator and an initial condition which as we shall see bears little importance to the outcome of the experiments. For the time being we will be concerned with the boundary conditions only. The three equations present regular singularities at the poles that owe their existence to the coordinate system. From a physical point of view the solutions sought have to be finite everywhere in the domain | < 2 Furthermore, the parabolic character of the predictive equations for potential vorticity indicates that a kind of insulating conditions L". O z = 15 (2.2.1) are appropriate at any finite boundary. Likewise, to preclude any flux of heat across boundaries one should select

26 ='TZ^ (2.2.2) However, the same insulating conditions could be obtained by imposing zero eddy diffusion coefficients at the boundaries: K, = K - K3 =0 and using some other form of boundary condition. Furthermore, at the poles a limiting process is necessary since the area vanishes there. To decide upon the most appropriate boundary condition we will now consider a more detailed analysis. For the axisymmetric flow the poles are natural (although singular) boundaries, on the contrary the equator has to be chosen a boundary because of limitations of the governing equations and data availability. This is an essential distinction which introduces profound differences in the search for suitable boundary conditions. 2. 2. 1 Boundary Conditions at the Poles Equation (2. 1. 26) is to be integrated over the domain 0< < <-. Having in mind that the thermal streamfunction YTis proportional to the mean temperature between 25 and 75 cb, any boundary condition in terms of z only would be too restrictive. Hence more suitable forms will express conditions on the first derivative at both ends of the interval. At this point it seems appropriate to consider the balance of terms in equation (2. 1. 26) since our main concern is global and seasonal variations, a suitable scaling magnitude for spatial derivatives is the radius of the earth a. Then the Laplacian is of order -2 w a a quantity which compared with q is about 160 times smaller.

27 Hence in (2. 1. 26) the second term on the left hand side balance almost completely the right hand side. To a good approximation we have: TZ= - 2 (Q Q)Q3Z) (2. 2. 4) which is the inner solution. Unless boundary conditions are chosen in such a way that they are compatible with this inner solution thermal boundary layers are bound to appear. Observational data does not support the existence of such layers and care must be exercised when choosing the boundary conditions for KTZ Q1Z and 3Z Returning to the complete differential equation(2. 1. 26) if we use a new independent variable x = sing it can be written as (I - x2> b-2xy' - ()= F (x) (2.2.5) where )() and F(x (QZ The homogeneous equation ~(I ~- X2) - 2 x - a(d >(XC) =O (2. 2. 6) has eigenvalues XAn,= CO4- ~+(n* ) 0,Z,4,... with Legendre polynomials',(x) as eigenfunctions. This choice of h implies an even y(x) which keeps some global invariants on the hemisphere, as will be shown further on. The equation then can be solved by a Fourier analysis technique

28 o0,(,x)-_ C, E(x)',=0 /'' evenl where F (x)=- cP (x) n=o even This procedure suggest that at x= 1 no other condition than finiteness is strictly required once the regular solution has been chosen. In the Fourier technique the choice of the regularly behaved solution at x = ~1 is made when the second solution of the Legendre equation is left aside. The problem is now how to perform this selection in a finite difference formulation. The differential equation (2. 2. 6) has two linearly independent solutions. A power series technique around x=l gives a first solution I, = I+ "-x + while the second is - aL[,- x] + E bn[-x] Hence a boundary condition of the form 2 ~ X=l would select the finite solution of the homogeneous equation. We still have to consider the contributions from the forcing function F (x). The particular solution of (2. 2. 5) can be expressed in terms of a Green's function G (x, x') as

29 ) —, )(x,xx') (x) d ux and since evenl it is seen that G (,'x) and its derivative are finite at x = ~1. Hence'p[(=+I) is finite when F(x) and its derivative exists at x = l~ which is always true, and we can use a condition like = M < oo x- I (2. 2. 7) to choose the regular sblution of (2. 2. 5). Returning to our original independent variable ~ relation (2. 2. 7) corresponds to'O a' - -M = From the finiteness of ^T one concludes then that a useful boundary condition, would be Meln 1 z0 (2. 2. 8) 2 A completely similar analysis gives the same condition as Y ---.. Once this choice is made as the boundary condition for equation (2. 1. 26) relation (2. 2. 4) suggests that the expression

30 (Q,-Q o 0 (2.2.9):an be used as a boundary condition for the prognostic relations 2. 2. 30) and (2. 1. 31) thereby avoiding the formation of thermal )oundary layers near the poles. Equation (2. 1. 26) can be written as y2 coi w a 1T 2 Q aking the limit as ~p —. and using L'Hospital's rule and (2. 2. 2) ve get 2 ^ ^^m ^, Q Qs a atpf 9'$~-^ tt' 2 (2. 2. 10) i relation that replaces (2. 1. 26) at the poles. We now turn to the predictive relations. Equations'2. 1.30) and (2.1.31) share a term of the form I;91[K ~ Q] a^co^^9 B^ [nM, Kco^ D Ahich can be written nd K aiQ e oK Ktean ta 3nd the finiteness of the expression demands that K tan%?-B remain rounded as -- 7/. This requirement and relation (2. 2. 9) gives as boundary conditions em 0=o or K 0 =~2.o,~__.~

31 The first alternative was chosen because it satisfies in a very simple way relation (2. 2. 9), is quite adequate for an axisymmetric solution and because a vanishing eddy coefficient is not a very realistic representation for the polar cyclonic activity which although sporadic and weak still exists. L'Hospital rule then gives and twrn a& =- 2\ ll< iQ t lt Then the two prognostic equations have as limit form over the poles aQ = i P. t -A(d -} ) abt - at P a \ z 3Z/ (2.2.12) aA z 2 Q aQ3z L aK a Q3Z H at -- 0. ay1 al a9 7+A(t,-i z\ 2 (33, 2) (2. 2. 13) The thermodynamic equation, (2. 1. 34) used to evaluate the meridional circulation, also presents this kind of term which by (2. 2. 8) is reduced so that the limit form of the equation for the pole is I Ka. * A2 Ct a., a

32 2. 2. 2 Boundary Conditions at the Equator Since any vorticity is by definition, at any level and time, integrating from pole to pole one gets I cet o d' =O IT/2 If the integration is over one hemisphere only so that to retain a zero total vorticity, at any level and time, over the hemisphere we must set indicating that convenient boundary conditions at the equator are at. _ 3 r = (2. 2.15) This also precludes heat fluxes across the equator, and makes i (Q.+Q3) c>Wo =

33 as can be easily shown from the definitions of Q. and Q3. From the potential vorticity expressions one deduces that Q' =2(' + )' ^a~~ge =2Q + ah~ ~(2.2.16) at the equator. Adding the equations (2. 1. 30) and (2. 1. 31) multiplying by cos and integrating from pole to pole we get bt f (,^n)C= -0 When the integration is performed over one hemisphere only, subject to (2. 2. 15) we obtain [: ^ ^ co+-.1-, ) o^\~~~~~' t^ ^Qc=O and to preserve the constancy of (Q,, + Q5z ) over the hemisphere we have to assume 1oFor arbitrary K, and K( this can be attained only if s o IF- n o (2. 2. 17)

34 which is another convenient boundary condition at the equator. Then by (2. 2.16) 2+ = -o 3- +2 - o E (2.2.18) at the equator. These boundary conditions imply symmetry conditions with respect to the equator as can be seen from the power series solutions of (2. 2. 5) around the equator. Equation (2. 2. 5) has as solutions around L=O,~ = I+,Ax I2+ +x'+.4. 3 +-: and ]: X[. )( + 2' 6 o 42 where CL, is arbitrary. Imposing or inerms ofo or in-terms of IP a =0s tP o

35 we are picking an even solution around the equator for AtZ. An odd solution is not possible because of the physical interpretation of JTU ( proportional to mean atmospheric temperature), and points out that it is impossible to cross the equator with a formulation in terms of the geostrophic streamfunction. Since F1T2 = 2 A ^t and is discontinuous across the equator, the symmetry of'Tz implies a discontinuity of Te or vicerversa. The diagnostic equation for TTZ at the equator becomes --' 2 (2. 2. 19) The potential vorticity equations are for I-'' _ _ _ K_ _ bt ^ o? as # d^^ 2LJ (2. 2. 20) aQor K-3. S^ IL ad the thermonami c (2.eq2. 2) and the thermodynamic equation _ aP4 RL J^ H2i I t 2 2 at CLZ a92 4Co ZZ- 2 CP Fb (2. 2. 22)

36 We have now given mathematical justification for the choice of boundary conditions at the poles. Conservation principles valid for the globe help to define them at the equator. The boundary conditions at the pole include vanishing meridional first derivatives of the two potential vorticities. However to state the finite difference solution in a most economic fashion such conditionsmust be supplemented. Specifically to reduce the integration matrix to a tridiagonal one it is necessary to require a second order symmetry as polar boundary condition: B"Qn= r 1 = o n= 13 (2.2.23) which by (2. 2. 16) also implies =a -0 -ap' -^ pi s n= l,3 (2. 2. 24) 2.3 ENERGETICS OF THE MODEL Under this heading we consider first the meridional transports of angular momentum and heat to continue with that part of the energy cycle linked to the zonal kinetic and available potential energies. 2. 3. 1 Momentum Transport In a quasi-geostrophic model with continuous vertical structure potential vorticity is defined from equation (2. 1. 12) as

:37 Q - fZ + 5 + at ($ at )(2. 3. 1) The transport of potential vorticity is obtained from this definition multiplied by the meridional velocity, v, and zonally averaged. If the transport of relative vorticity is expressed in terms of the convergence of momentum and the temperature is substituted for the streamfunction, by using the hydrostatic relation we obtain AU1 -, Ir~~ at3~ P(2,.3.2) The substitution for the fluxes of heat and potential vorticities by expressions in terms of the eddy coefficient for quasi-conservative quantities, K, whose uniqueness was proved by Bretherton (1966b), provide after some simplification cc_' |____3 F U }(l i.a _ _ _ OT ^Tl a) Li aM a P ^ a 7Z.'~~~~~~ "'(2.3.3) This equation shows that the vertical variation of the exchange coefficient K is linked to the evolution of the dependent variables. It also suggests that momentum transport can be modelled indirectly if parameterizations of the heat and potential vorticity transports are available and that these parameterizations are not completely independent of each other, but have to satisfy the integral constraints of

38 the transport of angular momentum. Assuming steady state conditions, integrating over the complete atmosple ric column and remembering that in steady state or for annual mean conditions the total divergence of angular momentum equals the surface torque exerted by the air on the earth one obtains: -ac, = K [F.] + ^ a 2 ^^~- Jo Q. p p o J (2. 3.4) Under the assumption that deviations from steadiness are small and considering that rs is proportional to the surface wind where surface westerlies prevail, the second integral must be negative and dominate the right hand side of the equation. The first integral is positive because under quasi-geostrophic conditions f >>}. If we now integrate (2. 3. 3) with respect to latitude from 2 pole to pole after multiplying by cos we get: J ^ a[^4+ i] +Po. a3<' aT o,'?ao 0 jI K say + ^o- ap'O.'a so ~,\ (2. 3. 5) This is a condition that has to be satisfied on any isobaric surface and is also valid for the hemispheric case if the boundary conditions imply a vanishing transport across the equator. Since quasi-geostrophy implies a misrepresentation in the low latitude regions and vanishing mean meridional motions,the important transfer of angular momentum by the Hadley cell is underestimated and one has to expect a large negative residual in the hemispheric balance of momentum. Within a quasi-geostrophic framework it is illusory to comply with relation (2. 3. 5).

39 Let us now turn to our two-level model. From the definitions of potential vorticities (2. 1. 18) and (2. 1. 19) one can deduce the relations - c^ { ] ]= [l0Qy]+ [;Yv]2 (2. 3. 6) - Q. coftp {[L 4[jc]C~ =[o~ jl- L;'J2. (2.3.7) that relate the transports of potential vorticities to those of momentum and heat. In obtaining these expressions use has been made of integration by parts and the vanishing advection of temperature by the thermal wind. In view of the parameterizations introduced by equations (2. 1. 28), (2. 1. 29), and (2. 1. 36) these relations can be cast in the form co- f { [ l ]c - K, Q - 2 (2. 3. 8) -2st [L ]i o -I K[ x 1 \ K Qn3 aoT~r yF ^ ^^ {^^ [Y1C*q K(2.3.9) If we multiply by a. COOZ and integrate these equations from pole to pole we obtain J {ax +f K2Fd4= 0 (2.3.10) JQ a<Kc ^Ydp=O J^T;: 1^ 2QN_~K2^t(M~q~ )(2.3.11)

40 These integral constraints show that the use of time independent values for K1, K2 and K3 introduces some inconsistency. In hemispheric models these constraints reduce to [ iK an t iK ea~ - A I~1 0 0 O 0 ^^ ^K lM^ p L 1 Cr% t)K a(f3Z-q2K2^;casPfd L -j [uff]z ISO ~ (2. 3. 13) The last equality is based on boundary conditions which isolate one hemisphere from the other: I, ('=O0)= L3(pO0=) 0. Addition of these relations gives f) - f- {Kg aqua + 5 J to Ccf( =O (2. 3. 14) 0 If we add equations (2. 3. 6) and (2. 3. 7) to obtain the vertically integrated divergence of momentum, by the same procedure we got the relation (2. 3. 4) for the continuous case, the two-level model gives after some arrangements and substitutions -E U2.ati (Ka + st b (Q Qi - ( — eyQton-Qs o) (2. 3. 15) a relation that must be satisfied under steady conditions or over a typical year. As far as the zonal flow is concerned this cannot be attained by adding a constant velocity to the flow at all levels. Such additions implies a change in vorticity so that the corrected solution

41 is riot compatible with the governing differential equations and the boundary conditions. Such a correction was evaluated for each year so as to make the atmospheric torque on the earth vanish, but is not included in the results presented later. The angular momentum transport will be evaluated from an integration, from a given latitude to the pole, of equations (2. 3. 8) and (2. 3. 9): Col [ = i2K FU.-~$1,= C0 p 41(K. aiQgtsZ-2K. arg C> The angular momentum transported poleward across the latitude circle k in the upper half of the atmosphere is, using the thermal wind L4U = - ~ a M,:- -c 5 - K1a^~ ~^ ^- (K -q qKpuc J^ad (2.3. 16) and below 50 cb we obtain a similar relation: T f/2 M1 _ (Kl t adosa- <,a. - o(.P.d1 9 ~~~~~~~~~~~(2. 3. 17)

42 giving a total amount for the whole column: m^^ ^ ^J 1a41 t BQtt +,,, )tXcayd(2. 3. 18) 2. 3. 2 Heat Transport By definition the heat transport is given by [H T] (\Ti2 ~PT= C C ( ) Using the eddy exchange parameterization (2. 1. 36) we get (HT)2- ^Cc or in terms of the zonal thermal wind (H T)c _ For the whole atmospheric column we finally obtain (HT)= 2 a -t c^ ^c' 9 (2. 3. 19) 2. 3. 3 Energy Balance Through 2. 1 we have described briefly the various approximations behind equations (2. 1.30), (2. 1. 31), (2. 1. 26) and (2. 1. 27). Several of these simplifications have been justified from an order of

43 magnitude viewpoint. Such an approach seems justified for short range prediction purposes, but it may not be suitable when integrations are performed over long periods of time. Lorenz (1960) has introduced one procedure which through energy invariants attains inner consistency of the equations and avoids unbalanced sources of energy. Following the same idea we will perform a detailed analysis of the expressions describing energy conversions and will match the outcome of the equations of motion and the thermodynamics first principle. This procedure will also help to define quadrature schemes for integration in pressure which in a two-level model are heavily truncated. In the model we are using, the available potential and the kinetic energy enter the only invariant since the static stability is not allowed to vary (Lorenz, 1960). Kinetic energy is evaluated more precisely from the zonal velocities directly, instead of vorticities, because those are the outcome of an integration which in general will tend to compensate errors. The model provideSzonal velocities at two levels and these will enter in a quadratic: form in the kinetic energy evaluation. Available potential energy, on the other hand, can be evaluated from isobaric deviations of the thermal streamfunction which is given at one level only. Again the form involved is quadratic. Hence we have more and better information for the kinetic energy than for the available potential energy. On these grounds it was decided to choose the kinetic energy expression and adjust the potential energy formula so that conservation of energy prevails in the system in the adiabatic, frictionless case.

44 The kinetic energy of the model is expressed in terms of the solenoidal component of the isobaric flow only. This is a consequence of the hydrostatic and geostrophic assumptions mentioned in Section 2. 1. a. If we choose a trapezoidal rule for a quadrature scheme in pressure we can define the kinetic energy for the atmosphere by where the integration is performed over an isobaric surface. Considering zonal means and their deviations the total kinetic energy can be expressed as K = K+ KE where K.- c /j2 IVA + as (2. 3. 20) In terms of the streamfunction after an integration by parts we obtain Kz= + 2ttt3ZjS S then d: - I' t t S dK1 _ ~. (sl

45 This expression is suitable for use with the vorticity equations of the model. Remembering that the winds! and $ in equations (2. 1. 13) and (2. 1. 14) are solenoidal the zonal mean of the advection terms can be written:a A/57(E+ ) 4 a ceC6 alp ((UZ tX *P ) (2. 3. 21) Noting that by the solenoidal condition and after some arrangements we can obtain LT - ces an A3)ICY o9 (2.3.22) Substitution of (2. 3. 21) and (2. 3.22) in the zonal mean of (2. 1. 13) and (2. 1. 14) gives ata asOco _ __ where L-1,3 S I I A I F 1- = l 2 i^ )_^f] iSWAP~itl) -4A,, r

46 Forming the expression for the time variation of the zonal mean kinetic energy and after some manipulation we obtain... _ f......F r dt - 22 act 2URY |[; Uix A S + f ^ 2AA + ~ t32(3L y.- 4)IS S -(2. 3. 23) The interpretation of each term in the right hand side is evident. The conversion of zonal available potential to kinetic energy can be written as 2 -~C('.)=,/w2l8 dS,.(2. 3. 24) and the dissipation of kinetic energy by friction D(K.)= t/ 2 + i (3da5 (2.3. 25) In order to express the conversion of eddy to zonal kinetic energy in terms of zonal mean quantities we will make use of the relations

47 i cQ ( ) < cos: |')z + 9s ( 4 zX(2. 3. 26) which are deduced from the definitions of potential vorticities (2. 1. 18) and (2. 1. 19) making use of (2. 3. 22). Then the eddy flux parameterizations (2. 1. 28), (2. 1. 29) and (2. 1. 36) allows us to write down the first term on the right hand side of (2. 3. 23) as S C(K,.X)~= -X( E^,K; aA 5 + 4 — /| KL dS (2.3.27) The thermodynamic equation (2. 1. 34) can be used to obtain an expression for the rate of change of available potential energy. Since this is related to the isobaric variance of the thermal field (Lorenz, 1955) as a first step let us take the isobaric mean of (2. 1. 34) and subtract it from the original equation to obtain at(X+Vt e(~/'4r) 0I H1 (2.3. 28) where primes indicate isobaric deviations.

48 Taking zonal mean of (2. 3. 28) we get at& CIC ^ a^ (~V^V -4), co 2 c, Ht According to (2. 3. 24), which links the zonal kinetic to the zonal potential energy the conversion, C (Az, K ) must be generated from the vertical velocity term. To obtain the same form as in (2. 3. 18) we have to multiply by f and integrate. In this connection we must remember that )2Z-= -c0 by conservation of mass so that' 02z - = i)2Z. Introducing (2. 1. 36) we obtain ^ = - C(A.Ae)- C(AK+G(A) where fCP H dS (2.3.29) C( (A2.A A) 42^ K/ u d-CS (2. 3. 30) and the zonal available potential energy must be defined as A, = 2eS'2~YT SJ (2.3.31) S

CHAPTER III DIABATIC HEATING 3. 1 INTRODUCTION The purpose of this section is the design of a heating model suitable for a zonally averaged, two-level model with a time resolution adequate to simulate annual variations. This heating function must be as close as possible to reality in order to minimize its effect as a source of errors and permit a validation of results by comparison with observations. We will not hesitate to make use of climatological information to obtain this goal and therefore the thermal forcing to be formulated is of little value for climate simulation purposes. In the standard form, the two-level model possess a constant static stability which fixes the value of the lapse rate and limits the representation of the infrared vertical fluxes to some function of the mean atmospheric temperature. Another shortcoming of the dynamic model, from a radiative point of view, is the location of the level where thermodynamic variables are evaluated. If the atmosphere is to be considered a grey body a temperature representative of the lower troposphere is more relevant than the 50 cb value, since the most active absorber and thick clouds are found there. However, Smagorinsky (1963), Saltzman (1967 and 1968) and Saltzman and Vernekar (1971) have dealt successfully with this problem and provide us with useful basic information. The driving force of atmospheric motions is the radiative. transfer of energy. As a consequence, radiative processes play a central role in the forcing of the dynamical model. Radiative processes 49

50 depend strongly upon temperatures at the earth's surface and in the atmosphere. The latter is specified initially and subsequently forecasted by the two-level model. The surface temperature has to be determined. For this purpose use is made of the heat balance of a thin surface layer with negligible heat capacity. Admittedly this is an idealization which is not supported by reality over the oceans, where solar radiation penetrates several meters while long wave radiation, evaporation and turbulent transfer of heat are actual surface processes. Thus our "thin" layer is several meters thick and has an important heat capacity. This complication will be disregarded in this attempt. As a matter of fact temperatures at 100 cb will be used as surface values. When the surface temperature is known, the diabatic heating can be obtained from the energy balance of the air column. Simplicity and the nature of the dynamical model dictate a climatological water balance independent of vertical motions which the model determines poorly in tropical regions, where the release of latent heat is important. 3. 2 HEAT BALANCE FOR AN ATMOSPHERIC COLUMN Considering a unit vertical column of air as a whole, five processes contribute to its diabatic heating (Fig. 3. 2. 1): i) Shortwave solar incoming radiation ii) Atmospheric longwave emission iii) Terrestrial longwave emission iv) Transport of sensible heat across the earth surface v) Latent heat released by condensation Following Saltzman (1968) these five processes will be (i) labelled Ha, i = 1, 2, 3, 4, 5. The analytical expression for some of them will be those used by several authors, Smagorinsky

51 (1963), Saltzman (1968) and Saltzman and Vernekar (1971), namely H(= X(1-r )Ro a a 0 H(2) + 4 a ( 1 + 2 ) T2 (3. 2. 1) (3) 4 a 4 where X is the atmospheric opacity to solar radiation, ra is the atmospheric albedo, v1 is the downward emissivity of the atmosphere, V2 is the upward emissivity of the atmosphere, r is the atmospheric absorptivity to long wave radiation, R is the incoming solar radiation at the top of the atmosphere, T2 the mean tropospheric temperature and T4 the surface temperature. Specific analytical expressions for the transfer of sensible heat and latent heat release will be considered in the next section. To evaluate the diabatic heating consider the vertical divergence of the fluxes through the whole column. At the top of the atmosphere the net flux is Ft- [(1-a)- r S(1-r) (1-X)Ro] + v2 a T + (1 ) T4 the upward net flux across the surface is Fb - R (1-r) (1-r ) (1-X) -v c T2 =b + a 1Ta4 The only important heat source to be considered is the release of latent heat by consensatimo processes H. Other heating a

52 x(l-rI )R~ va2T (5) ATMOSPHERE T2 H 2 a EARTHS SURFACE /, /11,///' / / //, // 777 /77 7777 7T7 4 (3) (4) valT2 H H 12 a a Figure 3. 2. 1 Energy fluxes related to the diabatic heating of the atmosphere. H(1) H(2) T4 S OCEAN 3 H(6) Figure 3. 3.1 Energy fluxes related to the surface layer.

53 processes affect only small fractions of the total mass of the column, Hence the heating of the whole atmospheric column is given by Ha= (Ft-Fb) ^a - ^tor Ha = X(1 ra - ra v+) RT - ( + + H or (3. 2. 2) H = H(1) - H(2) + H(3) + H(4) + H(5) a a a a a a The evaluation of the surface temperature T4 will be done using the surface heat balance 3. 3 SURFACE HEAT BALANCE Consider now a thin layer around the surface of the earth such that its heat capacity can be ignored. Energy fluxes to be considered are i) Shortwave solar radiation ii) Atmospheric longwave emission iii) Terrestrial longwave emission iv) Vertical transfer of sensible heat v) Evaporation vi) Subsurface heat conduction The time integration of the equations of motion is done with time steps of 12 hours and for such a short time the surface layer cannot be considered in thermal equilibrium, therefore it is necessary to use a layer with negligible thermal inertia. The vertical divergence of energy fluxes and the sources satisfy then a time independent continuity equation for energy.

54 If the surface under consideration is the ocean surface the sixth process mentioned above becomes of paramount importance and cannot be neglected. The oceans receive heat in low latitudes throughout the year, while in middle latitudes they store heat from the warm to the cold season. The excess of heat received in low latitudes is transported poleward by the general circulation of the oceans. This is a fact that should be accounted for in the model. In the oceanic case the balance of the thin layer presents an additional problem, due to penetration of solar energy to appreciable depths. This difficulty will be overlooked and the entire incident solar radiation will be assumed absorbed at the surface. The horizontal advection of heat within the layer is negligible since by previous assumption its heat capacity can be disregarded. Ignoring, in addition, surface heat sources (such as dew) the heat balance for the surface becomes (Fig. 3. 3. 1): [(- H(1) H(2) + H(3) + H(4) + H(5) H(6) ) =0 S S S S S S (3.3. 1) with H(1) = ( -X)(1 - r )( 1 - r) R (2) 4 S V1 2T (3) Hs -CT4 and where rsis the albedo of the earth surface. Substituting in (3. 3. 1) we obtain ( 1 - X)(1 - r )(1 - r) R~ +( a ( H -H = 5 ~~ (3.3. 2)

55 3.4. VERTICAL TRANSFER PROCESSES AND THE WATER BALANCE Eddy processes which transfer water vapor and sensible heat to the atmosphere are due to turbulence generated by thermal and mechanical instabilities. Direct observation of this phenomena has so far been made only on a local scale and global figures given in the present literature are determined as small residues required for balance in a climatological time scale. Therefore, these are not well known quantities, especially over the oceans where observations are scanty. There are some indications that seasonal variations of sensible heat transport are relatively small when compared with similar variations in other items of the energy budget (Davies, 1963). Some authors have used parameterizations which depend linearly on the temperature difference between the atmosphere and the earth surface. It is difficult to ascertain the representation of such formulations since they certainly bear no relation with mechanically generated turbulence and as far as convective processes are concerned they are difficult to evaluate from observations due to the intermittent character of convection in time and space. Because of these uncertainties and for convenience, annual mean values of heat transfer and evaporation for each.latitude will be used throughout the year, i. e., (4) b(4) Hs b(P)= - H s Here b is the annual mean value of the sensible heat transfer as a function of latitude, E the evaporation, L the latent heat for water.

56 This simple formulation has the additional advantage of making evaporation time independent, facilitating a formulation for latent heat release, a most important process in generating zonal available potential energy. A third turbulent transfer mechanism to be considered is the subsurface conduction of heat. This phenomenon is important only over the seas, since, by comparison, the heat capacity of continents is negligible. An analytical expression (6) k( T TD) (3. 4. 1) was chosen. Here TD is an empirical temperature to which a loose physical meaning can be attached only when TD< T4i. e., when the stratification is stable. It is better to consider k and T as empirical factors to be determined from observations and the nature of the surface. E is a factor, time dependent, which relates the mean temperature at 100 cb, T4, with the surface ocean temperature at a given latitude. The factor k, a measure of thermal conductivity, was computed from 1 1 1 k k ( i 200 2 2 2 250 3 1000 (3 2) (3. 4. 2) after Saltzman (1968) where k* is a constant and jl, j, j3 and j4 are the fractions of the globe covered by oceans, sea-ice, continents and snow-ice, respectively. The relative importance of the different kinds of surfaces was fixed using values provided in Laikhtman et. al. (1970).

57 This extremely crude parameterization endeavors to represent the ocean behavior as a heat storage for the atmosphere. Although nobody can feel satisfied with this approximation it gives fairly good values for the oceanic heat transport towards the poles, Fig. (3. 6. 2), and a qualitatively correct ocean behavior, Fig. (3. 6. 3). (5) The release of latent heat H(5 can be parameterized a using the water balance and the observed rainfall distribution in time and space. From rainfall values at a given time a set of factors was computed in such a way that they are proportional to the amount of rain and its integral over the hemisphere, at any time, equals one. Using then 7T H(5) m(f, t) Ecosp d a 0 the release of latent heat has a distribution similar to rainfall and in every unit time condenses as much water as has been evaporated. Hence, there is no variation of water storage in the atmosphere. Representing the last integral by I (5) = m I H ml a All the processes present in the heat balances have been parameterized and equations (3. 2. 2) and (3. 3. 2) can be written in their final forms Ha = X (1-ra) R - (v1+ V2)aT + FaT4 + b + mI aa -o 1 2)VT-+UT2 T4 b (3.4. 3) (l-X)(1-ra)(lr) R + Y - Tb -E -k(eTT O v~-T uT4 -bE k(ET4-TD)(..4

58 As already stated seasonal variations are the main goal of the study. However data availability and simplicity dictates the inclusion of time dependence only where it seems to be of importance. Thus R, T2 T4 TD' m and E will be considered as functions of time and space, while X, ra, r, r, v, v2 b and E are functions of latitude only; I is a constant. The values of k will depend on latitude through the fractional cover factors ji 3.5 COMPUTATION OF THE HEATING FUNCTION The amount of heat available to an atmospheric column with unit cross section per unit time is given by (3. 4. 3). Given T2, the only obstacle to compute H is the unknown T4. To determine the a surface temperature we can use equation (3. 4. 4). Reordering this equation according to powers of T4 we obtain T4 + k T4 + C 0 where C= -(1-)(1-r)(l-r)R - T4+ b+E-kT a s 1 2 D This fourth order equation has the following properties: the function of T4 so defined has always positive curvature since co is always positive, the only real extremum is a minimum which is in the negative T half plane. Thus the function has at most two real roots, 4 one of which is always negative. Henceforth, if we use Newton's method, a real positive guess is sufficient to pick the physically meaningful root of the equation. The coefficient of the fourth power is an absolute constant while the other two coefficients are functions

59 of.latitude and time. Considering the right hand side of equation (3. 4. 3), the first and last two terms are functions of latitude and time, while 4 4 the coefficients of T and T4 are functions of latitude only. In order to compute the heating function H, we need to interpolate a three quantities in time and latitude and three others in latitude only. Latitude interpolation will be performed using a spline scheme (see Appendix II) with vanishing first derivative as boundary conditions. This choice seems sensible for the poles where symmetry is required. At the equator it was chosen for simplicity, since physical reasoning is of no help. Time interpolation was also done using the spline technique choosing ends where the linear estimation of the first derivatives is accurate, i. e., near the inflexion points of the seasonal fluctuation for middle latitudes. The quantities L1 = 1 L2 = I ( V1 + V2 ) L3= a F were interpolated in latitude only, so that values for every five degrees were available. The quantities LT1 k E LT = - (1-X)(1-ra)(l-r )R + b + E - kT a s o D LT3 X(1-r )R + b + m I 3 a o were interpolated in latitude first and then spline coefficients were stored for interpolation in time as close as necessary.

60 3. 6 TUNING OF THE HEATING MODEL The diabatic heating of an atmospheric column has been expressed in terms of a set of constants and parameters whose values has to be determined. Some of them are easily found in the literature but have to be adjusted to the specifications of the model. There are others which are not available explicitly. However, the whole set of parameters has to be consistent, so that the present state of the earthatmosphere system is compatible with their values. They also have to satisfy some global constraints. For instance, the atmospheric heating integrated over the globe and the year, must be very small. There must also exist an annual balance between the space and the earth atmosphere system, and the oceans must receive as much heat as they release over the year. In general terms we can divide the set of parameters and constants in four groups: (a) Those determined by geographical and astronomical conditions: R (b) Those determined directly by climatological observations: b, m, E, E (c) Those determined from global constraints: k, TD. (d) Those determined by fitting observational results: r, V1, v2, X, ra, rs Once determinations of groups (a), (b), and (d) are performed an adjustment is done so that mean monthly observed values of surface and mean tropospheric temperatures satisfy the global constraints with good approximation. Some words about the use of mean monthly temperatures for tuning purposes are necessary. These temperatures enter the

61 heating function in a nonlinear fashion so the use of mean values is not strictly correct in the determination of the parameters. At any time: T. T.' i = 2,4 1 1 1 where Ti is the mean value over say a month and Ti the corresponding deviation. Temperatures affect the heating function nonlinearly through the long wave emission terms where it appears raised to the fourth power, so 4 4 -3 22 T,23 4 T 4 T +~4T T + T' +T' taking monthly average T 4 4 1 6 T' 2 Since we are using zonal averages daily and local extrema Since we are using zonal averages daily and local extrema are smoothed out (except near the poles), the error made when substituting mean monthly values is very small and for a month in which the mean daily temperature changes as much as 200K around a 2500 mean amounts to less than 1%. 3. 6. 1 Numerical Procedure Computation for the tuning of the zonally averaged heating model were performed in latitude with an interval of 10 degrees, with the first point at the equator and the last at the pole: values at latitudes 0, 10, 20,..., 80,90 degrees were computed.

62 Global integrations were performed using the trapezoidal rule, and global balances were satisfied with an error less than 1 &y/day. To use the tuned heating model in the dynamic time integration, values were interpolated to every five degrees of latitude. The interpolation scheme, designed to give a smooth variation, did not preserve integrals, so a better accuracy in the global balance would not be justified on these grounds and also because of the approximate knowledge of the magnitudes of the energy fluxes. 3. 6. 2 Determination of Parameters Incoming solar radiation R: This quantity was obtained from List (1951), Table 132, interpolating with a spline scheme (see Appendix II) for twelve equally spaced centered points along the year. To ease the computational procedure in the integration of the initial value problem, a fictitious time unit was chosen such that 360 of them are equivalent to 365, 25 common days. Thus the year is divided into twelve equal fictitious months, each composed of 30 fictitious days. Each fictitious day equals to 1. 01458 days so that the correction is small. However, Ro values were multiplied by this factor. Before applying this correction, Ro values were adjusted to integrate, over the globe and a year, to one fourth of the solar constant, taken as 2 1. 95 cal/cm min approximately. Radiative fluxes were changed to the new time unit by a slight modification of Boltzman's constant c Turbulent transfers and latent heat released were not modified due to their approximate character. Evaporation, E This process was taken as time independent. Seasonal variations of the zonal mean evaporation are not

63 TABLE 3. 6.1 Boundary layer vertical transfers (ly day ) Latitude Evaporation Sensible Heat Transfer N. H. S. H. N.H. S. H. 0 210 197 28 28 10 213 224 36 28 20 215 237 55 38 30 188 217 68 37 40 140 180 57 27 50 96 121 43 29 60 64 58 31 32 70 32 21 11 16 80 9 6 -16 -26 90 0 0 -36 -38 Source: Newton (1972)

64 as small as the sensible heat, but its parameterization is particularly difficult in a simple model as the one used here. Attempts were made to use the temperature difference between ocean surface and 100 cb level without success. Fig. 9. 15 from Newton (1972), gives an idea of the seasonal variations involved. They are somehow important in the NH subtropics, but nowhere reach the variations shown by precipitation in low latitudes. This approximation will not affect directly the heating function but will modify the surface temperatures which in turn enter in the heating determination. The values shown in Table 3. 6. 1 were interpolated from Newton (1972). Transfer of sensible heat, b: This process is considered as time independent and annual mean values for each latitude were taken from Newton (1972) for both hemispheres. The validity of this approximation is substantiated by Fig. 3 in Davies (1967), which shows seasonal variations of this term together with latent heat release for some latitudes of the Northern Hemisphere. Distribution factor for latent heat release, m: For the computation of this factor the basic assumption is a constant water vapor storage for the hemispheric atmosphere. The amount of water evaporated per unit time is condensed in the same time interval with a geographical distribution m dependent on the season. An additional, less stringent assumption is used, namely that latent heat release can be estimated through the rainfall, i. e., water precipitates as soon as condensation occurs. With these two suppositions at any instant the release of latent heat can be expressed as TF (5) tT (5) ~ (3.6.1)

TABLE 3.8. 2 Latent heat release normalization factor' m Latitude Jan Feb Mar Apr May June July Aug Sep Oct Nov Dec 90 N 0.137 0 114 0. 117 0.119 0. 087 0 103 0. 151 0. 160 0.194 0, 178 0. 166 0. 158 80 O. 171 O. i53 O.143 O. 146 O.137 O. 137 O. 192 O. 197 O. 240 O. 225 O. 205 O. 201 70 O. 285 O.242 O.234 O. 225 O.224 O. 228 O.313 O. 310 O. 333 O. 337 O. 331 O.317 60 O. 763 O.700 O.651 O.597 O.511 O. 526 O.524 O.572 O. 648 O. 730 O.779 O.813 50 O. 945 1.013 1.041 1.061 1.034 O. 994 O.908 O.863 O. 860 O. 871 O.896 O.929 40 1. 014 1.107 1.145 1.180 1.121 O. 983 O.827 O.770 O. 777 O. 824 O.887 O.971 30 O. 820 O.955 1.002 1.008 O.897 O. 811 O.787 O.816 O. 758 O' 721 O.731 O.739 20 O. 569 O.446 O.390 O.464 O.623 O. 800 1.009 1.032 O. 953 O. 871 O.779 O.633 c.n 10 1. 822 1.527 1.301 1.193 1.370 1. 542 1.715 1.877 2.035 2. 060 2.026 2.006 0 1.594 2.036 2.473 2.519 2.279 2. 000 1.584 1.267 1.092 i.096 1.150 1.256 0 1.613 1.704 1.759 1.765 1.706 1. 594 1.452 1.319 1.238 1.254 1.364 1.497 10S 1.613 1.539 1.393 1.224 1.071 O. 947 O.865 O.848 O.933 1.128 1.378 1.558 20 1.251 1.178 1.033 O.880 O.769 O. 697 O.650 O.625 O.663 O.807 1.025 1.195 30 O.576 O.603 O.674 O.756 O.821 O. 870 O.907 O.933 O.920 O.845 O.722 0 620 40 O.732 O.793 O.960 1.146 1.300 1. 418 1.510 1.570 1.538 1.357 1.067 O.830 50 1.012 1.032 1.108 1232 1.398 1.575 1.709 1.736 1.636 1.439 1.224 1.074 60 O. 737 O. 749 0.-'769 O.790 O. 808 O. 826 O. 849 O. 879 O. 896 O. 872 O. 810 O.755 70 O. 214 O.228 O.248 O.258 O.242 O. 215 O.199 O.218 O.255 O. 273 O.252 O.225 80 O. 037 O.039 O.042 O.043 O.041 O. 038 O.~37 O.038 O.042 O. 044 O.041 O.038 90 O. 012 O.013 O.013 O.013 O.012 O. 011 O.010 O.010 O.010 O. Oll O.011 O.012

66 Using the constant storage condition JH5) (cos d^ H(5)cos d ) cost dtP = cosP d (3. 6. 2) o o fbr any time the integration of (3. 6. 1) gives fm cosP dp = 1 The second assumption mentioned before, permit us to write (5) H ( Lr (3.6.3) where r stands for rainfall and L is the water latent heat. Then (3. 6.1), (3. 6. 2) and (3. 6. 3) give (5) H r m = fHn5) cosV dip Jr cos( doi and so m can be computed from a knowledge of r(H) for any time. Taking then r((p) for different months one can get m(., t), Table 3. 6. 2. Monthly rainfall data for the Northern Hemisphere were taken from London (1957) and for the Southern Hemisphere from Mbller (1951). Oceans surface temperature factor, E: This factor reduces the 100) cb zonal mean temperature to the oceanic surface

67 temperature at the same latitude and month. Using values of ocean surface temperatures, T s, as provided by Neumann and Pierson (1966) Table 14. 3, for the North and South Atlantic Ocean and 100 cb temperatures given by Oort and Rasmusson (1971) for the Northern Hemisphere and by Jenne et al. (1968) for the Southern Hemisphere, values were computed for each month and every 10 degrees of latitude. Subsurface condition, TD and k: One of the most uncertain factors in the heating budget is the subsurface conduction. This process is very important over the oceans. In comparison it can be neglected over the continents and ice covered oceans. To determine the seasonal behavior of this process the following procedure was applied. Using values of T2, T4, Ts from the sources mentioned above, the surface ss balance equation (3. 4. 4) was used to determine the amount of heat conducted downward, H, as a function of latitude and month of the year. We chose to give a form like (3. 4. 1) to the process, where k is determined by (3. 4. 2) and k. was given the value 10 ly day deg. Then, suitable values of TD(l, t)were computed. It must be stressed that k and TD have no physical meaning attached to them; they only serve the purpose of proper evaluation of H) s In very high latitudes where there is no ocean coverage this term is absent, and the surface budget was approximately balanced by modifications of the radiative parameters. Radiation parameters: All six radiation parameters, rs, ra F, v, v2, X are considered latitude dependent only. This represents a fairly good approximation for some of them like F, X Y2 ra, but for v1 and rs the approximation is less reliable in high latitudes.

68 The short wave parameters X, ra, rs, were evaluated using studies of London (1957) and Smagorinsky (1963) for the Northern Hemisphere, and of Sasamori et. al. (1972) for the Southern Hemisphere Results presented in these articles allow an individual determination of X, r, rs for the year as a whole and for individual seasons. For the long wave part, results of Smagorinsky (1962) were used for the NH. The study of Sasamori et. al. (1972) allows a determination of vl and a relation between v2 and F the atmospheric flux divergence. Assuming values for', which in the literature do not vary much, values for v2 were fixed so as to satisfy the observed values for the long wave flux divergence. In this process again monthly mean values of surface and tropospheric temperature as observed were used. Once all the parameters were computed as described, small adjustments were made to satisfy the global annual equilibrium of the atmosphere and the earth-atmosphere system within a 1 ly / day tolerance. The outcome of this process for the Northern and Southern Hemispheres is shown in Table 3. 6. 3. Satellite observations have provided evidence for lower tropical albedos than those shown, (Vonder Haar and Suomi, 1971, and London and Sasamori, 1971), but no detailed radiation balance based on them was available to us. The values shown are to a large extent in agreement with the already mentioned sources. 3. 6. 3 Characteristics of the lhned Heating Function The behaviour of the heating function can be appreciated from the following results. Using observed mean monthly temperatures

69 TABLE 3. 6.3 Radiation Parameters ra: atmospheric albedo r: earth surface albedo s r s ocean surface albedo (used only for oceanixc heat transport estimate) X: atmospheric opacity to solar radiation r atmospheric absorptivity to longwave radiation vl: atmospheric downward emissivity v2: atmospheric upward emissivity Northern Hemisphere Lat. r v v rar rss x EQ 0.96 1.35 0 77 0. 27 0.07 0.06 0.30 10 0.96 1.35 0.78 0.29 0.06 0.06 0.30 20 0.96 1.35 0.79 0.26 0.07 0.06 0.29 30 0.96 1.33 0.80 0. 28 0. 08 0. 07 0. 26 40 0.96 1.30 0.81 0.31 0. 08 0. 08 0. 27 50 0. 96 1. 29 0. 83 0.36 0.08 0. 09 0. 29 60 0.96 1.26 0.84 0.38 0.10 0.10 0.30 70 0.96 1.20 0.84 0.41 0. 20 0.11 0.36 80 0.96 1.20 0.85 0.42 0.63 0.63 0. 38 NP 0.96 1. 09 0.87 0. 40 0. 63 0. 63 0.39 Southern Hemisphere Lat. r ra r r x ^ ^^1 ^2 ra ss SsX EQ 0.96 1.35 0.77 0.25 0.07 0.06 0.30 10 0.96 1.33 0.78 0.25 0.07 0.06 0. 29 20 0.96 1.32 0.80 0.25 0.08 0.06 0. 29 30 0.96 1.35 0.82 0.25 0.08 0.07 0.28 40 0.96 1 37 0.84 0.34 0.05 0.08 0.28 50 0.96 1.39 0.88 0.42 0.05 0.09 0.34 C60 0.96 1.37 0.90 0.51 0.10 0.10 0.37 70 0.95 1.16 0.86 0.41 0.45 0.15 0.27 80 04.95 0.82 0.83 0.30 0.80 0.65 0.23 SP 0.95 0.66 0.70 0.26 0.80 0.70 0.23

70 at 50 and 100 cb and the evaluated parameters the annual hemispheric energy budget was computed. The results, including observed values, are shown in Table 3. 6. 4. The agreement is quite good showing that, in this sense at least, the modifications performed to satisfy global constraints have not introduced substantial changes. The energy budget for the whole earth-atmosphere system is shown in Fig. 3. 6. 1 as a function of latitude. When one compares with observations of Vonder Haar and Suomi (1971) both curves are somewhat low in the intertropical region with a net excess of absorbed energy. Towards the north pole the agreement is very close. In the Southern Hemisphere the absorption is slightly underestimated. Discrepancies in the tropical region are due to the different albedos considered. The mean planetary albedos in the model are 34. 2% and 34. 6 % for the Northern and Southern Hemispheres, respectively, while a value of 30% reported by Vonder Haar and Suomi (1971). This difference is mainly due to the intertropical region. The poleward transport of heat by the oceans was estimated using the zonal mean values given in Table 3. 6. 1 for the turbulent transfer of heat and water vapor. Oceanic albedos, shown in Table 3. 6. 3, were taken from Manabe (1969), and the oceanic coverage factor was used to estimate the amount of heat that has to be conducted across the surface. An integration in latitude provides the transport polewards, if no heat is transferred across the equator. Fig. 3. 6. 2 shows the curves thus obtained. As expected, both hemispheres shows a net loss which is larger in the Northern Hemisphere where the mean turbulent transfers are less representative of oceanic conditions. The uncertainty in observations of these fluxes is very large, especially

71 TABLE 3. 6. 4 N.H. S. H. Globe London Model Sasamori Model London & 1957 et. al. Sasamori 1972 1972 I 100 100 100 100 1 00 AA 17.5 20. 1 20.7 20.5 20.7 SA 47. 5 45.5 44. 6 44.7 44. 6 AR SR 335. 34.2 34. 8 34.7 34. 8 SR AE 153. 161. 152.8 159. AEU 56.5 60. 8 60. 7 66. 5 66. 5 SES 5. 5 4.6 65 4. 66 SE 114.5 117.3 111.8 115.2 111.8 AED 96.5 100.1 96.1 98.4 96. 1 SEA 109. 112.7 110.4 SHT 11. 6. 4.1 E 18.5 22. 24. AEU SES I \ AR SR SPACE /SEA ATMOSPHERE AE AA \ SHT E ___V____!- ____________ SURFACE SA AED SE

72 0 o 500 0 - maL o o0, Y AZSOUTHERN HEMISPHERE I 0/ o/ 3o o Vonder Haar & Suomi (1971) 100 ) * Simulation 0 — Tuned Heating Function 600 _ - -- L-, |,zO; 7;NORTHERN HEM I SPHERE 300 100-I 0O90 60 40 30 20 10 0 LATITUDE Figure 3. 6. 1 Energy balance of the earth-atmosphere system. Solar energy absorbed in the tuning is shown by the continuous line. Longwave emission in the tuning is shown by the dashed line. Values obtained in the simulation and observations as indicated. Units: ly/day.

~- Tuned Model Oceans cft: Ccr-r~ Cr ~ -- Simulation c1 D / 12 I Tuned Model OcAtm 9l^ // \\' — Simulation' /cl~-~. 0 -3- 0^ i f ^.ICD 0 II Palmen & Newton Bryan (1962) I'I \\ / Zo \ ~ /r\67 0 00.~(D 0 I FU 0I PO r'~ ~ ~i___ 1 W' SP 60 40 30 20 0 20 30 40 60 NP D f, c LAT ITUDE

74.with the low values for the tropical albedos (Vonder Haar and Oort, 1973). The profiles obtained should be compared with those evaluated with high albedos, for instance as shown by Bryan (1962). The order of magnituae and general shape of the curves are similar. There exists a maximum between 20 to 30 degrees and improvements in the right direction are to be expected from a proper account of the turbulent processes. Another source of uncertainty in the observed values is the equatorial flux, which seems to be from South to North. The time behaviour of the subsurface conduction (3. 4. 1) is presented in Fig. 3. 6. 3. Roughly speaking the oceans store heat from March to September in the Nbrthern Hemisphere, to release it in the rest of the year. In very low latitudes there is no cycle, the sea gaining heat throughout the year. The annual mean heating function of the tuned model for the Northern Hemisphere is shown in Fig. 3. 6. 4, together with those obtained from observations by several authors. The main discrepancy occurs at latitudes 20 and 30 where the model overestimates the cooling. Otherwise the curve falls among the observations, in particular at high latitudes. Fig. 3. 6. 5 shows how this mean heating is formed by the processes represented in equation (3. 4. 3). It is remarkable to what extent the release of latent heat by condensation determines the meridional variations of diabatic heating. In this regard a thermal forcing that does not include latent heat effect will be in error. Newell et. al. (1969) have estimated the heating for the four seasons. Fig. 3. 6. 6 compares the tuned heating function with their results. The modelling for the Northern Hemisphere is in good agreement with the observations except at the equator. In the Southern Hemisphere the fitting is poorer

75 Figure 3. 6. 3 Heat released by the oceans. Upper part corresponds to the Northern Hemisphere while the lower part belongs to the. Southern Hemisphere. Units: ly/day. Hemisphere. Units ~ ly/day.

76 LATITUDE 0 10 20 40 60 90 ------ Simulation \<'4 - Model Tuned ///// Palmen & Newton (1969 )\ ~} Sellers (1967) \ -6.- ~ Brown ( 1964) \. [ -~ Smagorinsky(1963) \.-2 10 Figure 3. 6. 4 Mean annual heating function. Results of tuning simulation and observations are shown as indicated. Units: 103 kj/(t sec).

77 I=J~~~~ C \ ^ >^C CD CD C: B -~~~~~~~~o oo lr=~~~~l ^ 1 1-1 I S\/ a' ~ O ~1 LI\ C) cs v- c Sa ) CO latitude. Units: ly/day. (f) v -- f~ L -- I I I I Figure 3.6.5 Atmospheric heat budget. The mean annual diabatic heating of an atmospheric column and its components according to latitude. Units: ly/day.

78 at all latitudes, with some evidence of the model marching somewhat ahead. Finally Fig. 3. 6. 7 shows the time and space variation of the diabatic heating for the actual atmospheric temperatures. 3. 6. 4 Diabatic Heating in the Simulation The dynamical model was unable to reproduce the observed temperatures and reached a different quasi-steady state. In the Southern Hemisphere the annual mean heating remained about the same at all latitudes with exception of the equator, see Fig. 3. 6. 8. This is not the case in the seasonal apportionment, Fig. 3. 6. 6. The hemispheric balance remained small. In the Northern Hemisphere the mean annual values remain fairly close to the tuned profile with discrepancies in low latitudes and some improvements near the tropic, Fig. 3. 6. 4. The seasonal fluctuations'- were retained with a good accuracy northward of latitude 1 0N. The hemispheric imbalance increased to 1. 7 ly day, a quantity small enough as to generate little variation between the second and third year of simulation. For instance, the maximum zonal wind at 25 cb in January increased from 34. 98 to 35. 03 m sec and the temperature at level 2, latitude 40, in January, changed only by one 1/100 of degree Kelvin. The available potential energy went from - 2 4837. 9 to 4838. 9 kj m. These amounts are in all cases negligible allowing us to state that the model reached a steady state for almost all purposes. The differences between the tuned heating function and the simulation function are most prominent near the equator. This is a consequence of the two main shortcomings of the model. First the diagnostic values of the eddy exchange coefficients were evaluated (Wiin-Nielsen and Sela, 1971) under the quasi-geostrophic assumption,

~'BP /I': slunfl'OauI snonui.uoo -e Aq umoqs a~i (6961) "'I I [TIAmaN iq sanlEA IeuoIT.AIasqo'uoITTnI uTfs aoqI IUOJJ sanl^A oJa 1 S0a10 pu-e soaI'!uiun% aOq uT soanlsA sa89 1pu T auII paqsBC'utLunoo3 DTJxaqd -soUS uae jo guuisaq zoi qaeTp aq1 Jo uoTinq.Jas.p IsuoseaS 9 9 onE 9'nIlTaJ V 03 0I 03 2 o~OI~~ I C) 00 c0 ( 0'Ji CC>~~~~ I~~~~I> Z) 6 L 00 CI. __

ap/ [. ~stUr'IaqdsL'oiuH uj'1l+no8 aGl ol ixed jaAOt aoq puB aaoqds.uauH UJaL44ON aql ol spuodsaaaoa lped uaddfl muunloc xi.aaqdsomLue u- jo ~u.-eaq ~.ijeqeT.C L'9' an. d 00.. c % n lz::w P- cj w 8 %.n Cy, s - co CD C:J O ) CJ1 C:.J ) i rN ~ ~'J N 40'-P n O O ~ O % C:> C) ~ ~ ~ ~ ~ ~~~~C "C9 ca c~~~ o 6~~ C~~ ~' bOs kaO O', OO Z,'C I~ I I IL I I~ I~~~ I Il Q~~~~~~~~~ ~I \I C) [i~8, \ I __ Ir o~'., t~~ % 08

81 100- oQI 0 O! I ~.. I NORTHERN HEM IS PHERE 100: 8s c' I,,~: SOUTHERN HEM IS PHERE * Tuned Model ro o 0Simulation 1I1 I I, I I I 1 90 60 50 40 30 20 10 0 LATITUDE Figure 3. 6. 8 Comparison of mean annual heating obtained in the simulation and the tuning. Upper part Northern Hemisphere. Lower part Southern Hemisphere. Units: ly/day.

82 so that they are unable to account for mean meridional circulations, which are the most effective agents of transport in low latitudes. Besides, the model itself is quasi-geostrophic and ignores all mean meridional circulations, as far as transport effects are concerned (see for instance relation 2. 1. 36). The second reason for the misrepresentation in low latitudes is due to the artificial boundary at the equator, which does not allow heat fluxes across. The energy balance of Table 3. 6. 4 is maintained to a remarkable extent. The largest modification amounts to. 5%. Modifications in the transports of heat are shown in Fig. 3. 6. 2, next to observational results of Palmen and Newton (1969) and Sellers as quoted by Lorentz (1967). As a final comment, Figures 3. 6. 2 and 3. 6. 4 show that observational values cover a wide range, specially if one considers the new values for tropical albedos. The heating function model we are going to use falls within these variations everywhere but near the equator, fitting observation quite well in middle and high latitudes of the Northern Hemisphere.

CHAPTER IV LINEAR SIMULATIONS Equations (2. 1. 30) and (2. 1. 31) were integrated over three years using diagnostic values for the eddy exchange coefficients for heat and potential vorticities. These depend only on latitude and appear in Fig. 2 of Sela and Wiin-Nielsen (1971). This integration was performed for each hemisphere separately. The same values of the eddy coefficients were used in both cases due to a lack of information about the southern half of the globe. A schematic flow chart of the computation procedure is shown in Fig. 4. 0. 1 for the time step n. The computation was performed with a grid interval of 5 degrees of latitude (550 km. approximately) and a time step of twelve hours. The initial condition, of little importance for our purposes, was a motionless and isothermal ( T2z = 273 K ) state. The spin up of the model took about 200 time steps ( 3 a months ) after which the seasonal fluctuation repeated almost exactly in the second and third simulated year. This periodicity also serves as a control of truncation errors. The computational stability of the time integration was insured by an implicit scheme (Crank-Nicholson, see Appendix I). In exchange for this advantage the thermal forcing and the dissipation had to be used with a half time step lag with respect to the other terms in Equations (2. 1. 30) and (2. 1. 31). This lag of 6 hours in the forcing is assumed to be of little consequence on the results which have a much longer typical time scale. 83

START IN XT IAL CONDITION Equations Used Hn n V 4,4-3 2ZT'4Z 1 4,4.4 an an 2. 1.18 i Z' 3Z 2.1.19 n+Q l Qn+l 2.1.30 z' 3Z 2.1.31! -nal | 1 {2.1.26 OUTPUT n+l-.n 2U n —n-l Figure 4. 0. I Flow chart for computation of linear simulations.

85 The basic cycle in the computation was the following: i) With Qlz' Q3,' qn given by initial conditions or a previous time step, equation (3. 4. 4) provides the surface temperature, T4. Then equation (3. 4. 3) gives the diabatic heating of an atmospheric column, H2, and the potential vorticities definitions allow the calculation of the relative vorticities 7 n and. lz 13i ii) Equations (2. 1. 30) and (2. 1. 31) are integrated, using boundary contitions described in Section 2. 2. The numerical procedure is developed in Appendix I. Thus we obtain n+1 and Q3z iii) Equation (2. 1. 26) is integrated with boundary conditions (2. 1. 10) n+1 and (2. 1. 22) and by the technique described in Appendix I, to give T. Tz iv) If desired at this point one can evaluate the three dimensional flow, temperature at 50 cb, heat and angular momentum transport and energy conversions by the relations deduced in Chapter 2. This was done once every six time steps (once every three days). v) Return to i) and repeat as long as necessary. A three-year simulation took 112 sec of CPU time in an IBM 360/67 computer. 4. 1 NORTHERN HEMISPHERE SIMULATION Results for the Northern Hemisphere are represented in Figs. 4. 1. 1 to 4. 1. 13. 4. 1. 1 Temperature Field The model provides temperature values for 50 and 100 cb, which are shown in Fig. 4. 1 1 and 4. 1. 2 respectively. The field at 50 cb represents a drastic improvement over that obtained by Sela and Wiin-Nielsen (1971) hereafter referred to as S-WN. The results show

86 NP 230 \ 240 60 250 26 30 ^ ^ 270~O EQ J J F M A M J J A S O N D Figure 4. 1.1 Thermal field at 50 cb for the Northern Hemisphere in the linear simulation: Units: K. NP// / 269240 30 200 Fiur_ 4..' aea iue411bta 0. J F M A M J J A S O N D Figure 4. 1.2 Same as Figure 4. 1. 1 but at 100 cb.

87 OK.~1 O. 270 0 o 00 250~- LU — 20 -20 Jan - 40 230 -._ — -50 0 30 60 NP LATITUDE Figure 4. 1. 3 Meridional temperature profile for July and January at 50 cb, Northern Hemisphere. Simulation values shown by dots and circles. Observed values after Oort & Rasmusson (1971) shown by the continuous and dashed lines. ~C 300 -~.,~-%o_ 30 300^U^-^ ^ 3 "'%'..July - 20 10 280 \ XJanJ - 260 1- _ Observed Values * - 20 from Oort & Rasmussort 1971) 3 240 -- O Simulation 0 30 60 NP LATITUDE Figure 4. 1.4 Same as Figure 4. 1. 3 but for 100 cb.

88 that most of their discrepancies with observations were due to the oversimplified thermal forcing. However the weak meridional circulation in low latitudes manifests itself in excessive thermal gradient, especially in the summer. See Fig. 4. 1. 3 for a comparison with observations. At this level the model produces a larger temperature gradient in middle and high latitudes in winter as if the eddy exchange were not as strong as observations suggest. Surface temperatures, shown in Figs. 4. 1. 2 and 4. 1. 4, seem to be better simulated than 50 cb values. Again low latitudes show disagreement to be blamed on a weak Hadley cell in summer. In Figs. 4. 1. 1 and 4. 1. 2 it is interesting to note that the date of maximum temperature lags the summer solstice by a larger amount the lower the latitude is. This is a consequence of the oceans behaviour shown in Fig. (3. 6. 5). 4. 1. 2 Zonal Motion The model specifies the zonal wind at 25 and 75 cb (the 100 cb wind is obtained through a linear extrapolation in pressure). These are shown in Figs. 4. 1. 5 and 4. 1. 6. The flow at the upper level simulates the observed variations in some aspects such as the trade winds appearance in summer and the intensity of the jet stream. On the other hand polar easterlies are present during a short period of time and reach higher latitudes than shown by observations. The maximum wind in summer is too low by about 10 m/sec and does not move to lower latitudes in winter. At 75 cb, Fig. 4. 1. 6, the zonal flow shows a summer anomaly due to the inconsistency introduced by the time independence assumption of the exchange coefficients (see 2. 3. a). The winter

89 LAT 1 ^ X 10-1 20 60 30 2 300 J F M A M J J A S O N D Figure 4.1.5 Zonal flow at 25 cb for the Northern Hemisphere in the linear simulation. Velocity correction for annual angular momentum balance: 3. 4 m/sec. Units: m /sec. LAT 6 1 1 J F M A M J J A S O N D Figure 4. 1. 6 Same as Figure 4. 1.5 but at 75 cb.

90 magnitudes are correct, but in summer, easterlies appear at all latitudes. The extrapolation at 100 cb shows the same feature. The figures improve if we correct the zonal velocities by a constant amount (3. 4 m/sec) such that the torque on the earth surface vanishes, when integrated over the hemisphere and the year. At 75 cb we obtain some westerlies between 25 and 50 degrees of latitude. However, such correction alters the fields of vorticities and therefore is illegitimate. In comparing with S-WN there is improvement in the magnitudes of the zonal flow. Most of the defects just mentioned are also present in their work and to this extent can be attributed to the values assigned to the exchange coefficients. 4. 1. 3 Meridional Circulation The meridional circulation can be shown through the meridional flow or the vertical motion since they are related by the continuity equation (2. 1. 38). Fig. 4. 1. 7 presents the vertical motion in the isobaric system of reference, so that negative values indicate ascent. There exists a three cell configuration all year long with some variations in intensity but little in geographical location. When these results are compared with the observations of Oort and Rasmusson (1971), hereafter referred to as O-R, two discrepancies are apparent. First the tropical ascent is too weak and in middle latitudes the upward motion is too strong. Second, seasonal variations of the extension of easterlies and westerlies are not well simulated. S-WN did not obtain a three cell system which again must be due to the simplified thermal forcing.

91 NP E I — 2(^J I I I I I ____I I 60 t -- 20 o # -20t 60 20 30 m _20 EQ I l u i 1 1 1 1 I I J F M A M J J A S O N D Figure 4. 1. 7 Vertical velocity at 50 cb for the Northern Hemisphere in the linear simulation. Units: 10-6 cb/sec. NP 60 0 0 30 J F M A M J J A S O N D Figure 4.1. 8 Heat transport by the atmosphere for the Northern Hemisphere in the linear simulation. Units: 1011 kj/sec.

92 4. 1. 4 Transport of Heat and Momentum The transport of heat shown in Fig. 4. 2. 8 compares quite well with observed estimates of the annual mean as provided by Lorenz (1967). This is a new check on the suitability of the heat transport parameterization. The seasonal and geographical variation is also in good agreement with results reported by O-R. Angular momentum transport is not well modelled, particularly in summer. In winter months this transport shows the deficiencies of the Hadley cell. Fig. 4. 1. 9 shows the actually computed curve for January (dashed line) and the variation adjusted so as to get zero transport across the equator (continuous curve). This last one is a very good representation of the annual mean both in magnitude and latitudinal variation, see Lorenz (loc. cit). In summer, the transport becomes negative at all latitudes and the adjustment decreases their magnitude but does not alter the sign. This behaviour is in agreement with the zonal flow features already commented. It seems that the diagnostic values of the eddy exchange coefficients may be adequate for winter conditions but not for modelling summer features. 4. 1. 5 Energetics Energy amounts and conversions between them were computed once every third day. Fig. 4. 1. 10 shows the tine variation of zonal available potential and kinetic energies computed by equations (2. 3.31 ) and (2. 2.20) respectively. To judge the quality of the modelling we can compare with results of observational studies by WiinNielsen (1967), Krueger, Winston and Haines (1965), Newell et. al. (1969) and Oort (1964), hereafter referred to as WN, KWH, NVDFK

93 30 20 Jan 10 i 30 60 NP LATITUDE / Figure 4. 1. 9 Angular momentum transport for the Northern Hemisphere in January. Dashed line shows the values obtained in the simulation. Continuous line are values corrected as to obtain zero cross equator flux. Units: 1015 ton m2/sec2.

94 and O respectively. WN based his calculations on a one-year data sample for the Northern Hemisphere including the atmosphere below 10 cb. KWH used a five-year data sample containing information of the 85 and 5C cb levels only and for latitudes north of 200N. NVDFK collected data from a number of sources covering the globe and up to 10 cb. Oort's article is a review paper and includes results of several authors for the Northern Hemisphere. Available potential energy showed a wave with the maximum in the first week of February. The springtime decrease occurs with a larger rate than the fall build up. The amplitude of the variation 2 2 amounts to about 3500 kj/m around a mean level of almost 5000 kj/m In a qualitative sense this is in excellent agreement with results of WN. The mean level is somewhat high compared with Oort's estimate of 4000 ~ 1000 kj/m2, but the discrepancy is larger with results of WN and NVDFK which are around 3000-3500 kj/m. However, the annual mean is very sensitive to the value assigned to the static stability, and the streamfunction constant, f0. In particular the latter has been given arbitrarily the value 104 sec which is used in models with cartesian geometry. Had we considered all latitudes weighted accord-1 ing to area a more proper value would have been 0. 729 x 10 sec which reduces the value of Az to about one half of its original magnitude. Therefore the value obtained is within the range that can be expected from the model. Besides, the poor modelling of the Hadley cell tend to increase the available potential energy. The amplitude seems too large by a factor of three. In this respect the results of KWHI are helpful because they have been evaluated for latitudes north of 20 N and exclude most of the tropical regions. The seasonal fluctuation of available potential energy for the 85 to 50 layer is reported

95 6 - / 2 0 J F M- A M J Jt A S 0 N D Figure 4. 1. 10 Simulated variations of zonal available potential energy (A ) and zonal kineti% energ (KZ) along the year for the Northern Hemiphere. Units: 10 kj/m J F M A M J J A S O N D Figure 4. 1. 11 Simulated variation of the generation of zonal available potential energy (G) and dissipation of zonal kinetic energy (D) along the year for Northern Hemisphere. Units: watts/m2.

96 2 with an amplitude of 1200 kj/m, which compares well with the value obtained from the model for the whole atmospheric column. With respect to S-WN's results, the available potential energy variation shows a smaller amplitude and higher values in the summer months but a somewhat higher mean value. In other aspects they are similar. The kinetic energy, see Fig. 4. 1.10, also represents a simple wave variation. The maximum lags that of the available potential energy by about 12 days. The amplitude reaches about 850 kj/m, about the same value as the annual mean. Observed values provided by WN are in very good agreement with the simulation in all respects. This is also valid for values reported by O and KWH. This is an indication that the zonal kinetic energy of the tropical regions is correctly represented. Results of NVDFK do not agree with those cited above, they are about one half smaller. The results obtained represent a considerable improvement over those reported by S-WN both in amplitude and mean value. The generation of available potential energy is represented in Fig. 4. 1. 11 beside the dissipation of kinetic energy. Their maximums occur approximately by January 15 and February 30 so there exists a 1. 5 months lag as in S-WN's results. It is to be noted that WN, KWH, and S-WN, computed generations and dissipations as residual terms, while in our case equations (2. 3. 23) and (2. 3. 19) were used. The mean values of both quantities over the year agree with estimates by O. The seasonal fluctuation of the generation has an amplitude similar to those reported by NVDFK and somewhat smaller values than those of KWH. Results of WN are not directly

97 comparable because they are in terms of generation of total potential energy, but do not contradict the outcome of the simulation. The general features of the seasonal fluctuation of both quantities reported by S-WN agree with our simulation, with the exception that neither the generation nor the dissipation attain negative values in our case. The three energy conversions that are present in the model are shown in Fig. 4. 1. 12. That from zonal to eddy available potential energy, C ( A, AE), is related to the meridional transport of heat and its maximum is around January according to KWH and WN. The simulation peaks in mid February. The annual mean agrees with values reported by 0, KWH and WN; those of NVDFK are again lower. With respect to S-WN's results, there is a slight improvement, mainly in the summer months. The conversion between eddy and zonal kinetic energies C (KE, KZ ), presents fluctuations which do not agree with those reported by WN and NVDFK. Both authors find that a maximum is attained in the fall while this simulation and that of S-WN show a maximum in winter. The location of the minimum is less clear. For NVDFK it is in the summer, being positive, while in WN's results it appears in winter and is negative. This transfer was evaluated through equation (2. 3. 21) and therefore presents one basic difficulty because it is a small residue obtained by subtraction of two relatively large quantities. Furthermore these quantities depend on parameterizations which are not valid in tropical regions (which cover about one half of the earth). At this point there is one indication that tends to support our parameterization. KWH have estimated energy diagrams for each season from 20~N to the pole. Taking the annual mean diagram and

98 C (A,A E) A C(KKEKz) 2 C(A zK) J F M A M J J A S O N D Figure 4.1.12 Simulated annual variations of the conversions between zonal and eddy available potential energies C(A A ); eddy and zonal kinetic energies C(KE, K \ and between zonal vailbble potential to zonal kinetic energies C(, KZ) for the Northern Hemisphere. Units: watts /m2. AZ KZ 3.~1.0* 4000~ 1000' 0. ~ 0.2* 800 ~ 300' 0.5 ~ 0.2 G D 2.1 4948 -0.7 818 0. 3.0 ~ 1.0* 10.4~ 0.2 2.8 1.2 AE KE Figure 4. 1.13 Simulated energy box diagram in a typical year for the Northern Hemisphere, zonal available potential energy (A ) and zonal kinetic energy (K ) in kj/m2. Generation (G), dissipation-D) and conversions in watts/m2. Observed values after Oort (1964) with asterisk.

99 introducing the value for the zonal kinetic energy dissipation of say 0. 5 watts/m2, a generally accepted value which must remain unchanged or increase when the tropical regions are omitted, it turns out from the principle of conservation of energy that a value of 1. 16 watts/m for C ( KE, Kz ) is obtained. Applying the same procedure to each season we obtain D(KZ) + C(KZ, AZ) = C(KE' KZ) Winter 0.5.83 = 1.33 Spring 0.5.77 = 1.27 Summer 0.5.40 = 0.90 Autumn 0.5.74 = 1.24 Year 0.5.66 = 1.16 This shows that the modelling is in agreement with the observation in extratropical latitudes. S-WN obtained negative values for this conversion during the spring and summer months. This is contrary to results of NVDFK and KWH, but presents some agreement with the analysis of WN. The same piece of information provided by KWH is useful in appreciating the results shown for the conversion of zonal available potential to kinetic energy, C( Az, KZ ). When compared with hemispheric averages the values shown have the wrong sign and are about one order of magnitude too large, but estimates of KWH show that when the tropical regions are not included, the values obtained in the simulation agree, in magnitude, sign and seasonal fluctuation, with observations. S-WN provided an annual mean value which agreed quite close with observations. Nevertheless, they were unable to obtain a

100 realistic field of vertical motion, a fact that detracts from their result; Finally, Fig. 4.1. 13 shows the annual mean energy diagram obtained from the study together with observed values from Oort (1964). The agreement is good except for the two conversions just discussed. In the light of the results of KWH one carn conclude that the model simulates the middle and high latitudes quite well as far as the energetics is concerned. 4. 2 SOUTHERN HEMISPHERE SIMULATION When dicussing the modelling of the Southern Hemisphere it is necessary to keep in mind that the eddy exchange coefficients were estimated from Northern Hemisphere data so that to some extent we are assuming the same cyclonic activity in both places. The simulation was completely analogous to the northern counterpart and it differs only in one aspect related to the presence of the Antarctic continent. The tuning of the thermal forcing function was performed with surface temperatures as given by Jenne et. al. (1968), instead of 100 cb values. 4. 2.1 Temperature Field At 50 cb, Fig. 4. 2. 1, the thermal fields behaves very much like the Northern Hemisphere's. Again the maximum meridional gradient fluctuates anomalously, from about 30~ in the winter to 40~ in the summer. Similarly the maximum temperature date is between the end of July to mid August in low latitudes. The field at 100 cb, Fig. 4. 2. 2, shows some substantial differences with respect to its northern counterpart. The oceanic coverage produces a more homogeneous regime. This is apparent

101 SP 230 60 240 250 260 27 EQ Llll l l l l 1/l J F M A M J J A S O N D Figure 4.2. 1 Same as Figure 4. 1. 1 but for Southern Hemisphere. SPI. 220 60 260 280 30: 25 I60 30 — EQ\I I I X I I I 1 I400 EQ J F M A M J J A S O N D Figure 4.2.2 Same as Figure 4. 1.2 but for Southern Hemisphere.

102 OK- ~250 July -20 - 270 ~"-50 ~ oo.,4.."5Q 2350- "o^^ 23I0I - o.. 0 0 30 60 NP LATITUDE Figure 4. 2. 3 Same as Figure 4. 1. 3 but for Southern Hemisphere. Observed values after Jenne et al. (1968). OK I C 26 0J M' 30 " OV -30 300 o — -20 240 Jenneetal. (1968) Simulaion Jan -610! 10 260 JuLlyAT E 6 \ %\ \ - -20 e 4 Observed Values from Soutn 30 240 Jenne et al. (1968) \ — 40 o Simulation 220 - \ -60 0 30 60 NP LATITUDE Figure 4. 2. 4 Same as Figure 4. 1.4 but for Southern Hemisphere. Observed values after Jenne et al. (1 968).

103 from a comparison of isotherm 280 K, for instance. In general the whole field in high latitudes follows a smaller amplitude variation; the larger latitudinal gradient is due to the orography effect included in the tuning of the heating function. The model reproduces the flat minimum winter typical of the Antartic stations. Fig. 4. 2. 3 and 4. 2. 4 compare latitudinal variations with observation for the months of January and July. The same characteristics discussed in Section 4. 1. 1 appear, namely, misrepresentation in the low latitudes and better adjustment in surface values. 4. 2. 2 Zonal Motion Figs. 4. 2. 5 and 4. 2. 6 show the zonal flow at 25 and 75 cb respectively. The seasonal fluctuation is somewhat less than in the Northern Hemisphere in agreement with the smaller seasonal change in the thermal field. This agrees with observations. However, the wind maximum is smaller and its seasonal variation in latitude is 6 months out of phase with observation. The trades behave properly, and there are no polar easterlies. At 75 cb the simulation is in fair agreement with observation. Polar easterlies begin at 650S in summer and 72~S in the winter according to Van Loon (1972). In this respect the model attains better agreement with observations in the Southern than in the Northern Hemisphere. The extrapolated flow at 100 cb does not show substantial differences between the two hemispheres except that the area covered by westerlies is smaller in the austral case.

104 SP 10 ~ 60 30 20 301 10 EQ t-~,,, J F M A M J J A S 0 N D Figure 4. 2.5 Same as Figure 4. 1.5 but for Southern Hemisphere. Velocity correction 3. 6 m/sec. SP 0 6a~~~~~~~0 60 -! 0 - 30 EQ... l,. III! J F M A M J J A S O N D Figure 4. 2. 6 Same as Figure 4. 1. 6 but for Southern Hemisphere. Velocity correction 3.6 m / sec.

105 4. 2. 3 Meridional Circulations The only observational study that contains information about the zonal mean meridional circulation of the Southern Hemisphere is Newell et. al. (1969). This source reveals that the Hadley cell is much stronger than the other two meridional cells in the fall, winter and spring and almost vanishes in the summer. Fig. 4. 2. 7 shows that the seasonal fluctuation is not reproduced and that the tropical circulation is too weak. The three cell structure is very steady in location and intens ity. 4. 2. 4 Transport of Heat and Momentum Fig. 4. 2. 8 shows the transport of heat by the atmosphere through the year. A comparison with the Northern Hemsiphere, Fig. 4.1. 8, shows that the major difference is a more moderate seasonal fluctuation, a consequence of the oceanic coverage. The mean level of the transport is in agreement with results of Sellers (1965), but is somewhat larger than the results by Newton (1972). This last author also presents seasonal fluctuations which are about one third of those in the Northern Hemisphere. The angular momentum transport, Fig. 4. 2. 9, shows the same defects that appear in the Northern Hemisphere. The main difference is a smaller maximum, a fact that does not seem supported by observation, Newton (1972). 4. 2. 5 Energetics Seasonal variations of energy amounts and conversions are shown in Fig. (4. 2. 10) to (4. 2. 12). It is apparent that the whole energetic picture is characterized by smaller seasonal variations than the

106 SP 0 -2o~ 60 20 3 20 30 EQ I I I I J F M A M J J A S O N D Figure 4. 2. 7 Same as Figure 4. 1. 7 but for the Southern Hemisphere.

107 10 3_3 0 EQ..... I J F M A M J J A S O N D Figure 4. 2. 8 Same a-s Figure 4. 1.8 but for Southern Hemisphere. 20 July 10 \ 30 60 90 LATITUDE -10io Jan Figure 4.2. 9 Angular momentum transport for the Southern Hemisphere in January and July. Both corrected for zero cross-equator flux. Units: 105 ton m2/sec2.

108 Northern Hemisphere. A peculiar feature is the early maximum in generation of zonal potential energy with respect to the summer solstice. The maxima of available potential and kinetic energy follow with time lags very much as those of the other hemisphere. The mean levels of generation and the two energies are in very good agreement with observed values, Fig. (4. 2. 13). However, the aforementioned time lag is not supported by observations, Newell et. al. (1969). This discrepancy can be traced back to the thermal forcing shown in Fig. (3. 6. 6). It is seen there that the heating function in the simulation (small circles) tend to underestimate the correlation between temperature and the diabatic heating in the fall and winter months. In low latitudes this is probably due to an incorrect assessment of the release of latent heat and its seasonal fluctuations. The only source for this data was Mbller (1951), who presents the information on a seasonal rather than monthly basis.

109 8 6 A Z 4 2 Kz F M A M J J A S O N D Figure 4.2. 10 Same as Figure 4. 1. 10 but for Southern Hemisphere. 32 J F M A M J J A S O N D Figure 4.2. 11 Same as Figure 4. 1.11 but for Southern Hemisphere.

110 i LC(AZI, ie, 3 1 I I I I I I I C( KI I J F M A M J J A S O N D Figure 4. 2. 12 Same as Figure 4. 1. 12 but for Southern Hemisphere. AZ KZ 1.9' A 0.16' 4300~ 620~ 2.0 4 0.1 0.3 4498 -0. 663 n. ** 1.2* 0.353 2.8 1.1 Figure 4.2. 13 Same as Figure 4. 1.13 but for Southern Hemisphere. ()bserved values (with asterisk ) after Newell et al (1969).

CHAPTER V SENSITIVITY STUDY The governing equations of the model contain some parameters introduced by the quasi-geostrophic approximation, others are needed for the representation of sub-grid processes, and finally, exchange coefficients are used for heat and potential vorticity parameterization. These parameters were assigned magnitudes according to average values suggested by observations. The results of the simulation must be considered according to their sensitivity to variations of these mean values. In this chapter we will report results of simulations performed with parameter values changed by reasonable amounts, (Table 5.0.1). We will be concerned with the quasi-geostrophic constants, fo and 5, experiments #1 and #2, the sub-grid constants E and A, experiments # 3 to #6, and the three eddy exchange coefficients, experiments #7 to #15. 5.1 STREAMFUNCTION CONSTANT -4 -1 The value adopted for f was 10 sec, corresponding o to Coriolis parameter at 43. 3 degrees of latitude. An additional -4 -I simulation was performed with f = 0. 6 x 10 sec1 a value of f at 23. 6 degrees, other factors unchanged. The variation of f introduces modifications in the values 2 2 of X2 and q. In the last case the dependence is quadratic and 2 q reduces to about one third of its original value. This implies that the baroclinic effects represented by q2YVT in the potential vorticities 111

112 TABLE 5.0.1 Exp. # Parameter Standard Value in experiment in study value 1 f 10-4 0.6 x 10 4sec o 4 -2-2 ~2 C 2 3 m sect -6 -1 3 | 3x10 6 1.Ox10 sec^ ~~~4 E ~ ~9. 0 x 10- sec 4 ^. 9,0x106 sec~1 5 A 0.6x106 0. x 106 03 10sec-6 -1 6 A 1. 0x10 sec 7 K2 (*) 100 shift south 6 2 -1 8 K2 j 2. 0x10 m sec 9 K2 (**) 10 K1 & K3 (') K1=K3 from (*) 6 2 -1 6 2 -1 12 K1 | (*) K1 2. 0 x0lm sec K3 from (*) 13 K1 100 shift south, K3 from () 14 K (:) 100 shift south, 3 K1 from (*) 15 K3 K3 from (44), K1 from (*) (*): values taken from Table 5. 5. 1 (44* ): values taken from Wiin-Nielsen & Sela (1971)

113 become smaller, although still dominant. As a direct consequence the available potential energy decreases proportionally to f. This strong 0 dependence indicates that the relatively high value of this energy (by about 20%) must not be of great concern, for a change of only 1% in f would give the observed magnitude of the mean zonal available potential energy. Of greater importance is the modification in the meridional gradient of potential vorticity since with the new value for 2 q 2contributions from the beta and thermal terms become comparable. As far as the vertical motion is concerned equation (2. 1. 37) indicates that a decrease of f will decrease the importance of the 0 contributions from horizontal advection and local changes of temperature. Consequently the tropical convection increases in magnitude. In January the maximum in trades doubles its value. On the other hand, the high latitude ascent becomes weaker and covers a narrower latitude belt. This improved correlation between the vertical motion and the diabatic heating changes drastically the conversion of energy between zonal available potential and zonal kinetic from -0. 8 to + 1. 1. This corroborates our former statement about a weak Hadley cell being the reason for the negative value of this conversion. A rough guess indicates that a decrease of f by 20% would give a value of C(AZ, K ) 0 Z Z about equal to the observed one. However, such a change brings in other undesirable effects. The use of a smaller f produced a decrease of the 0 zonal flow. For instance, the maximum wind in the westerlies diminished by 10 m sec at 25 cb and by 13 m sec at 75 cb in January. The fraction of the globe covered by easterly winds increased and consequently the transport of momentum became worse. This is an

114 2 effect of the weaker thermal forcing (which is proportional to 2 ) and the relation between wind and streamfunction. The kinetic energy decreased to about one half of the standard case. However the dissipation increased, because in the summer months the zonal flow at 75 cb reached larger absolute values (they are mostly negative ) than in the winter. The temperature difference between pole and equator decreased by about 20 degrees in winter and 10 in summer at 50 cb. About the same change occurred at 100 cb. As a consequence the heat transport diminished to about one half. In conclusion, f is an important parameter, and the 0 results are quite sensitive to its value. However, with the present formulation the value- assigned seems adequate. Changes in f will tend to deteriorate some results while improving others. 5.2 STATIC STABILITY The second quasi-geostrophic constant is the static stability, 4 2t-2 r, for which a value of 2 m sec t was chosen. Typical values provided for the months of January and July by Gates (1960) are Layer (cb) Jan July 30-50 3.73 3.66 50-70 2. 26 2. 20 The value chosen seems thus somewhat low. A simulation was done -_ 4 2-2 with c = 3 m sec t, experiment #2. The value of cr enters 2 2 linearly in the denominator of the constants X2 and q. For this reason an increase of 50% in static stability has many of the effects

115 of a decrease of f by 40%. For instance, both increase the geo0 graphical domain of easterlies and worsen the simulation of the angular momentum balance. The meridional temperature gradient decreases, but only by a few degrees. The influence of r on the different processes that generate vertical motion is the same. Hence, variations in r are not so selective as those of f. In January the tropical maximum ascent went down by about 30% while the high latitudes ascent diminished by about 40%. Because the temperature changes were restricted mainly to high latitudes there was a higher correlation with the vertical motion and the conversion of energy from zonal available potential to zonal 2 kinetic attained a value of -. 15 watt/m. Kinetic and available potential energies decreased to 2700 2 and 350 kj/m, respectively, and the same happened to the generation 2 that reached 1. 5 watts/m. In general, the sensitivity of the results to this variation was smaller than in the f case, although still of considerable magnitude. 5.3 BOUNDARY LAYER FRICTION The transfer of momentum between the atmosphere and the underlying surface was parameterized in terms of a constant ~. In Section 2. 1. d we showed that this factor was proportional to the wind speed near the surface, V4, to a drag coefficient, c d which depends on the roughness of the surface, and other factors which are less variable. ~-6 - Using typical values for V4 and cd a value E= 3 x 10 sec was chosen. To study the influence of the uncertain value of E two additional -6 1 -6 -1 simulations were done with ~= 1 x 10 sec and = 9 x10 sec, experiments #3 and #4.

116 Variations in ~ have only slight consequences for the energetics and the thermal field. The most important changes occur in the zonal flow. As an example, changes in the kinetic energy are 6% between simulations using the extreme values of &. Analogous figures for available potential energy, C(A, AE) and dissipation are 1%, 12% and 40%. The latter is the most significant change in the energy cycle. The total pole-equator temperature difference varies by 10 K at most. -6 -1 Changes in the zonal flow are such that with ~= 9 x 10 sec a remarkable improvement has been obtained in the angular momentum balance. The flow at 25 cb shows a maximum of 33. 4 m/sec in -1 0 January and 11 m/ sec in July, both at 40 N. There exist westerlies at 75 cb all year around. In July they occupy the belt between 20 and 50 degrees. At 100 cb there are westerlies between 30 and 70~N in January but they disappear in July. The zonal velocity correction necessary for a long term hemispheric balance of angular momentum decreases from 3. 4 m/sec in the prototype simulation to 1. 9 m/sec1 -6 -1 in the case with =- 9 x 10 sec The effect of an increase in ~ on the zonal flow is to decrease the amplitude of the seasonal variations keeping the maxima almost invariable, i. e., it intensifies the summer westerly flow without much change of the winter circulation. From this result, one can infer that in the summer the zonal circulation is frictionally driven while in winter the eddy exchange of momentum overshadows the interaction between the air and the earth surface. The meridional circulation is increased by a very small amount so that even a difference of one order of magnitude in ~ gives

117 an increase of 4%in the maximum ascent in the equatorial region and about 10%o in the middle latitude upward motion. In conclusion the value of E hardly effects the results of the simulation with the exception of the zonal flow which improves when is increased. The optimum value of & seems to correspond to a drag coefficient near the upper bound of those recommended in the literature, Cressman (1960). 5.4 INTERNAL FRICTION The internal friction coefficient A was shown in (2. 1. d) to depend on a kinematic viscosity v, representing the effects of subsynoptic processes, and the temperature T2z squared. The latter fluctuates by about 30% with most of the variance coming from the high latitudes. Estimates of the dynamic viscosity,,, have been made by Palmen (1955) and Riehl (1951). This information gives A -6 -6 -1 values of 0. 6 x 106 and 1. 3 x 106 sec, respectively. In the standard simulation we have used the value derived by Palmen. -6 -1 Two experiments were run, the first with A = 0.3 x 10 sec -6 -1 and the second with A = 1. 0 x 10 sec. In general, changes generated were small and the last case produced results in better agreement with observation. The zonal flow experiences only slight modifications: a decrease in west winds. Temperatures were also changed by small amounts in high latitudes. The temperature difference between pole and equator at 50 cb in January decreased by 6%. The most important variations were in the energy balance. The available potential energy decreased from 4950 to 4400 kj/m. The kinetic energy diminished by 160 units and the conversion C(A, K ) was reduced to -0. 4 watts/m.

118 The value of this parameter does not play a significant role in the outcome of the experiments although it helps to improve some of them to a limited extent. 5.5 EDDY COEFFICIENTS The model includes three eddy exchange coefficients. One for transport of heat, K2, and two for transport of potential vorticities, K1 and K3. All of them were taken as time independent and assigned values used by Sela and Wiin-Nielsen (1971) which were determined on the basis of a diagnostic study by the same authors (Wiin-Nielsen and Sela, 1971). Since the study was done under quasi-geostrophic conditions values of the coefficients were determined for latitudes between 25~ and 800N. The basic data extended over the whole year 1963. Monthly mean values for each latitude were computed. They represent therefore a rather smoothed distribution. However, the fact that the coefficients represent zonally averaged effects decreases the importance of using monthly mean values. Due to the small size of data sample monthly mean values were averaged and the annual mean used in the simulation. This additional smoothing filters out seasonal changes, and the extent to which they are important is shown in the results of the simulations discussed in Chapter 4. Between the Equator and 25~N of latitude the magnitudes were fixed by subjective extrapolation with no observational evidence to support them. The values used in the simulations are presented in Table 5. 5. 1. Since each exchange coefficient is a function of latitude the numerical study of their influence is not straightforward. Based on the results of Wiin-Nielsen and Sela (1971) we will assume that the meridional variation has one maximum in mid-latitudes, and we will

119 TABLE 5.5. 1 EDDY EXCHANGE COEFFICIENTS Latitude K3 K2 K1 6 2 -1 deg. 106 m sec 0 1.3 0.1 0.5 5 1.3 0.1 0.5 10 1.3 0.1 0.5 15 1.3 0.2 0.5 20 1.3 0.4 0.5 25 1.3 0.5 0.5 30 1.3 0.6 0.5 35 1.3 1.3 0.5 40 1.5 2.2 0.4 45 3.5 3.3 0.6 50 6.6 5.0 1.1 55 7.4 6.0 1.8 60 5.8 5.4 2.4 65 4.5 4.3 2.6 70 3.7 3.5 2.7 75 2.9 2.5 2.2 80 2.0 1.3 1.5 85 1.0 0.6 0.7 90 0.0 0.0 0.0 Source: Figure 2 from Sela and Wiin-Nielsen (1971).

120 try to assessthe importance of its magnitude and location. It is important to keep in mind that the transfer coefficient for quasi-conservative quantities is unique (Bretherton, 1966 a) and that the difference we make here between K, as a heat transfer coefficient, and those of potential vorticity, K1 and K3, is due to the different role they play in the numerical model and to the diagnostic evaluation used to determine their values. We will start with the heat exchange coefficient and continue with the potential vorticities. 5. 5. 1 Heat Eddy Exchange Coefficient The K2 coefficient is used to evaluate the vertical motion from equation (2. 1. 37) and in the energy conversions between eddies and the zonal component. It is irrelevant as far as the basic dependent variables are concerned. Because the integral constraints that relate the eddy exchanges of momentum and heat are not satisfied, changes in K2 do not alter K and K3. According to equation (2, 1. 37) the vertical motion is determined by four contributions: the local change of temperature, the product of the meridional variations of K2 and temperature, the product of K2 and the Laplacian of the temperature, and the diabatic heating. The local change of temperature makes, normally, a very -6 -1 small contribution (less than 0. 5 x 10 cb sec ). When the diabatic term reaches a level of 50 ly/day, it contributes about 12 x 10 cb/sec. Thus the thermal forcing guarantees ascent in the tropical regions and polar descent, unless K2 or its gradient happen to be very high. Typical contributions in mid-latitudes for all the terms are:

121 s 2 _ R 4WI - VK2 VWYt2 z 2 Tz 2f c 2z 4f t o p 0 0.5 8. 1. 12. (5.5.1) 6 1 where the figures represent contributions in 10 cb sec. Three experiments were run changing the values of K2 as shown in Fig. 5. 5. 1. In experiment #7, the standard profile of the coefficient was shifted 10 degrees equatorward. In experiment #8 a constant value equal to 2 x 10 6m2 sec was chosen. A profile with smoother variations was chosen for experiment #9 so that the gradient of K2 was different from the standard case. Actually, profile 9 is taken from Wiin-Nielsen and Sela (1971), Fig. 13. From relation (5. 5. 1) and the profile of K2, as determined by these authors, it is apparent that a constant eddy exchange coefficient is not a good representation. The outcome of the experiments can only be considered in terms of vertical motion and energetics since the coefficients for potential vorticity were not modified. The most important consequence of using a constant K2, experiment #8, were the effects in tropical regions where descent replaced upward motion. The cause of this drastic change is the large increase of the third term in (5. 5. 1) which overcame the diabatic heating contribution. The region covered by the polar front ascent shifted about 10~ to the south, narrowed and the maximum ascent reduced to one half. The energetics associated with the vertical motion remained almost invariable indicating a compensation between low and high latitudes effects: the Hadley cell becomes indirect, weak and

122 6,. i-: _, A% 5 or 30 60l90 4 7, \STD I2 I 3 I> I,~ ~ 2 em...s_ " mes e a....IseeOmes I'/ \ L / \.. 0 30 60 90 LATITUDE Figure 5. 5. 1 Profiles of exchange coefficient K2 in experiments #7, #8 and #9 beside the profile in the standard simulations. Units: 106 m2/sec. OK 270 250..- T12 STD *** * IO 230 0 30 60 NP LATITUDE Figure 5. 5.2 Temperature profile for January at 50 cb in experiments #10, #11 ~nd #12 bes-ide the profile in the standard simulation. Units: K.

123 coversless area, while the polar cell gains intensity and becomes more extense. Conversions expressed in terms of K2 remain very much the same probably because the value chosen is not far from the global mean. Experiment #7, with K2 values shifted towards the equator, showed that the vertical velocity profile has the middle latitude zero determined almost completely by the peak of K2. This is due to the dominant role of the gradient of K2 in equation (5. 5. 1) and suggests, because the observations do not support a large gradient of vertical motion around 500 of latitude, that the maximum of K2 must be less sharp than used in the simulations. The shifting of the profile towards lower latitudes produces an increase in its mean value. This increase is felt in the energetics and the three conversions showed values notably larger, C(A, K ) =-2. 2, C(A, A)= 4. 1 and C(KE, K)= Z Z Z' z E' Z 2. 6 watts/mr. Experiment #9 used values of K2 with a lower maximum and gentler variations. It presents also larger values in tropical and polar regions. The results show that this is a more realistic profile for middle latitudes where smoother transition and smaller values for the ascent are obtained. In low latitudes the larger K2 tend to decrease the ascent by small contributions of the third term in (5. 5. 1). In high latitudes the diabatic term dominates. Since the mean value of K2 is smaller in this experiment, the conversions associated with it decrease: -2 -2 C(A, AE) =. 8 watts m and C(KE, K) = 0. 3 watts m. Because -2 of the weaker Ferrel cell we also obtain C(A, K ) = 0. 3 watts m. z' z In conclusion the K2 shape has strong influences on the vertical motion in middle latitudes and through its mean value

124 determines the energy conversions from the zonal to the eddy components. 5. 5. 2 Potential Vorticity Eddy Exchange Coefficients The fact that potential vorticity is specified at two levels makes the discussion more cumbersome so we shall restrict it to the more outstanding features. K1 and K3 will be modified keeping K2 invariable. A first step in the study of the diagnostic eddy coefficients for potential vorticity is the evaluation of the importance of their functional dependence on the independent variables pressure and latitude. In this respect our efforts will be guided by some considerations of Green (1970) and the results obtained by Wiin-Nielsen and Sela (1971). As far as the pressure dependence is concerned equation (2. 3. 5) is of some help. For the two-level model it originates relations (2. 3. 10) and (2. 3. 11) which can be written as Y/2 coS K (f+f1) + (K2-K) q2 z d 0 1 a30 2 aS (5.5. 2) 2 (/2 Cs2- 2 cos2 K3 aa(f+) + (K3-K2) q a dW = (5 ^^r~/2~~~~~~~~~ g (5.5.3) Since, except near the equator, f >i for any level, they indicate that the eddy coefficients must, on the mean, decrease with height: K3 > K2 > K1. In particular this must happen in the baroclinic

125 regions where the temperature gradient is large. By the same argument equation (2. 3. 15) indicates that the annual mean behavior associates westerlies with regions where K3 > K1. Equation (2.3.3) also suggests an exponential variation with the square value of the pressure wherever the meridional temperature variation is important and fairly constant with height. In terms of exchange efficiency, large values of the coefficients in the lower troposphere have been explained as a consequence of the presence of the critical level (Green, 1970). The values calculated by Wiin-Nielsen and Sela (1971) are in agreement with these considerations as can be seen from Table 5. 5. 1. The difference (K3-K1) is large between 40 and 70 degrees of latitude. Experiment #10 shows the consequences of using K1 and K3 equal to the diagnostic values for level 3 shown in Table 5. 5. 1. This modification represents a substantial increase of K1. Under these conditions equation (5. 5. 2) cannot be satisfied. The model force the dependent variables trying to decrease the magnitude of the inconsistency. For instance, the meridional temperature gradient is decreased by 14~K in January and 6 K in July bringing the annual mean of avail2 able potential energy down to 2000 kj/m. The small value of the temperature gradient almost destroys the ascent in mid-latitudes forced by the contribution of the second term in equation (5. 5. 1). As a consequence, very large values of C(Az, Kz) are obtained (1. 4 watts/m2). The zonal flow is dominated by easterlies, the trades reaching large magnitudes. This that can be interpreted as an effort.to satisfy equations (5. 5. 2) and (5. 5.3) which if added give c s 123) a1s (Q1 +Q ) +(K3 -K1 (a''s)' d 0 ^~~~~~~13 (5.4)

126 In this experiment the second term in the brackets vanishes so that Q1 + Q3 = 2(f + 7 2) has to be negative, particularly in low latitudes, to satisfy the relation. The very large easterlies at 100 cb generate 2 a high boundary layer dissipation (4. 2 watts/m ). In conclusion, the eddy coefficients must decrease with height especially in mid-latitudes where westerlies are the observed prevailing winds near the surface. The next experiment was related to the latitudinal dependence of the coefficients. K1 and K3 were taken independent of latitude, equal to 106 m /sec and 3 x 106 m /sec, respectively. In winter, the absence of the gradient of the exchange coefficient shifted the zonal wind maximum to 65 N. Likewise, the westerlies maximum moved to 650N and 700N at 75-and 100 cb respectively. The trade winds covered a broader latitude belt. In summer, easterlies covered the hemisphere at 100 and 75 cb. Temperatures changed comparatively little. The total pole-equator difference in January was reduced by 40K at 50 cb. Polar easterlies seem to be related to small values of K3 since they are not present at any time of the year. The larger average value of K3, and to a lesser extent of K1, increased the conversion C(KE, Kz). The latitudinal variation of K1 and K3 seems thus to be of importance mainly in the zonal flow structure. The next experiment was done with K = 2. 0 x 106m2sec and K3 as in Table 5. 5. The results corroborate that a large value of K1 generates easterlies, shifting the maximum west wind polewards. Theeas theral field was modified decreasing the pole-equator contrast by changes in the very low and high latitudes (Fig. 5. 5. 2).

127 The last experiment with K1 was performed with the values shown in Table 5. 5 shifted by 10 degrees of latitude equatorwards. Changes in the zonal flow included a larger area covered by easterlies and a weaker maximum (24 m sec at 25 cb in January), probably caused by the global increase of K1. The most important effect, nevertheless, was the displacement of the maximum wind from 45 N to 25~N at 25 cb and the tendency to generate a secondary maximum around 80 N. The cause of the second feature is related to the decrease of the gradient of K1. Due to the decrease of zonal speed, kinetic energy is reduced -2 by 50% and C(KE, Kz) diminished to 0. 37 watts m. We conclude that the location of the maximum zonal flow is determined to a large extent by the profile of K1. To determine the influence of the location of the maximum of K3 another experiment was run with K3 shifted by 10 degrees equatorwards. Its results showed that the influence of the K3 profile is less drastic: it displaced the jet by 5 degrees only. It also increased the zonal flow maxima by about 14% at 25 cb, 25% at 75 cb and 45% at 100 cb. This differentiated addition gave a substantial improvement in the momentum balance. The thermal field remained almost unchanged, and the energetics showed only small changes. Finally a last experiment was done with K3 taken from the results of Wiin-Nielsen and Sela (1971), Fig. 12. This represents a sharper maximum and a lower global value. The only substantial change was in the zonal flow which diminished in intensity and wider areas were covered with easterlies. 5. 6 AN IMPROVED CASE To substantiate the influence of some of the modifications

128 just discussed an improved case was evaluated. For this case the following modifications were done: i) The K3 profile was shifted 15 degrees of latitude towards the equator so that its maximum was located at 40 N. ii) The K2 profile was modified to adjust to that reported by WiinNielsen and Sela (1971) and shown in Fig. 5. 5.1 as experiment #9. iii) The boundary layer friction was increased by making ~= 10 sec-1 iv) The three eddy coefficients were modified in the neighborhood of the poles so as to make them approach unity at the pole itself. The outcome of such an experiment is shown in Fig. 5. 6. 1 to 5. 6. 5 for those variables significantly affected by the aforementioned changes. The thermal field was modified only in the proximities of the pole by the larger coefficients. The temperature at the pole in January increased by 4. 5 K at 50 cb and 100 cb; in July it increased by about 2 K at both levels. The zonal flow improved because of the shift of the K3 profile and the stronger friction. The wind maximum moved to 40~N in better agreement with observations, Fig. 5. 6. 1, but still kept its latitude throughout the year. At the same time, the area covered by the easterlies decreased and some westerlies are present in summer in midlatitudes at 75 cb, Fig. 5. 6. 2. The velocity correction necessary to attain angular momentum balance, globally and over the year, decreased from 3. 4 to 0. 9 m/sec. The winter intensity of the zonal flow remained very much the same but it became stronger in summer. The vertical velocity, shown in Fig. 5. 6. 3, presents modifications due to the change in K2. The subtropical descent became

129 NP 60 20 330 ~~~~~~~30 / 2o o0 EQ L J F M A M J J A S O N D Figure 5. 6.1 Same'as Figure 4.1.5 for the improved case. NP 60o-\ ro/ j5 0 Q I I I I I J F M A M J J A S O N D Figure 5.6.2 Same as Figure 4.1.6 for the improved case.

130 NP ~ 20 60-a 30 30 _ EQ I 0 J F M A M J J A S 0 N D Figure 5.6. 3 Same as Figure 4.1. 7 for the improved case. NP 60 5 <10 30 EQ J F M A M J J A S 0 N D Figure 5. 6. 4 Same as Figure 4. 1. 8 for the improved case.

131 slower, a change not supported by observation. The ascent associated with the polar front narrowed and decreased in intensity improving agreement with observation. The strong meridional gradient around 50 N disappeared. However, the seasonal fluctuation is still too weak. The change in K2 introduced some modifications in the magnitude of the heat transport but not in its geographical configuration, Fig. 5. 6. 4. The momentum transport kept its negative values at all latitudes during summer while it attained very large positive values in winter. After correcting for zero flux at the equator there was a maximum at 30 N of about 83 x 1015 ton m2/sec2, for January about 2.7 times the observed annual mean (Lorenz, 1967). The mean annual energy box diagram, Fig. 5. 6. 5 shows important differences with that of Fig. 4. 1. 13 but only some of them can be considered improvements. The low value of C(Az, AE) can be explained by the lack of a tropical mean circulation but the large value of C(Az, Kz) indicates that the subtropical descent became too small.

132 A K 2.2 0.3 0.6 4700 920 1.9 0.5 Figure 5. 6. 5 Sare as Figure 4. 1. 13 for the improved case.

CHAPTER VI EDDY EXCHANGE PROCESSES We have seen in Chapters 4 and 5 that the diagnostic coefficients are able to model most of the important features of the axisymmetric circulation. However, the summer conditions in the zonal flow are not well reproduced and the wind maxima do not appear at the right latitude nor oscillate with the season. Some of these features are difficult to simulate in a quasi-geostrophic model because they depend heavily on the tropical circulation. But some improvement can be expected from a more realistic parameterization. The macroscale transport processes are mean meridional circulations and eddy fluxes. The first process is essentially ageostrophic, and it is excluded from the governing equations of the model, but not from the diagnostic values of the exchange coefficients. The eddy fluxes can be separated into those due to the very long waves and those due to the synoptic scale perturbations. The importance of this division is related to the fact that a parameterization by means of eddy coefficients is not very suitable for the ultra-long waves so that their effects are misrepresented in the modelling. Nevertheless, the diagnostic eddy coefficients include their effects. Had the coefficients been determined as functions of the internal variables of the model the ultralong waves would be in error to a larger extent. From observational studies it is apparent that all kinds of fluxes have important seasonal variations. Therefore, if the eddy exchange coefficients are representations of the eddy activity it seems natural to expect them to have a strong time dependence. When these 133

134 variations are determined by the current state of the simulation a feedback mechanism is included which increases the utility of the model for climate change studies. In the past, the problem of eddy flux parameterization has been closely linked to results of baroclinic instability studies (Charney, 1959; Green, 1970; Kurihara, 1970; Saltzman and Vernekar, 1971; Stone, 1972 a and b). An exception i;s Smagorinsky (1964) who estimated the meridional transport of momentum through the lower boundary stress. This technique is useful under steady conditions only. Green (1970), using a mixing length approach, also elaborates on a finite amplitude justification which rests on the importance of the steering level and the validity of Eady's optimum slope for finite displacements. Both features are related to small amplitude studies, but stand a better chance to be correctly extrapolated into the finite amplitude domain than the phase speed and growth rate provided by a linear baroclinic theory. We will suggest an extension of some results by Green (1970) and Stone (1972 b) which relate the eddy exchange coefficient for quasi-conservative quantities to the most relevant features of the mean zonal field. With this goal in mind we will digress to study the nature of the eddy exchange process and its applicability to the case in which we are interested. 6. 1 NATURE OF THE EDDY EXCHANGE Eddy exchanges are the macroconsequence of the excursions of a population of microelements. Every one of these possesses characteristics related to the environment where its migration started. They are modified, relatively slowly,to those of the new environment.

135 The idea of exchange coefficients rests on the assumption that a sufficiently stable mean field can be defined. Its macrovariations are due to the action of a large number of small eddies having a short lifetime with respect to the time scale of the changes in the macrofield. In other words, it is assumed that the frequency spectra of the transfer presents a gap. Being more specific, the transfer of mean zonal properties, which vary with a typical time scale of three months, must be performed, in the bulk by the transient waves and to a lesser extent by the stationary ones. It is important to note that in contrast to the interpretation in the kinetic theory of gases our case is related to quantities that change slowly (quasi-conservative). There is not an exact conservation followed by a catastrophic exchange. For a quasi-conservative quantity one can distinguish between three time scales. First, the very slow variation of the mean zonal field, Tr. At the other end, the fast displacement takes the microelement from its origin to a substantially different environment in a typical time Tj. In between is the relaxation time, T, required to smear out the differences between the element and its new surroundings. This is a measure of how conservative the process is. It is crucially important that the first time scale is very long compared with the other two. The displacement time must be shorter than the relaxation time. Hence we require > > j (6. 1. 1) If the quantity in question is not conservative the second inequality is not satisfied and we cannot be certain that we are transferring anything at all.

136 The first requirement is satisfied when dealing with the seasonal variations, and the daily cycle has been smoothed out (then TF - 3 months). Even for the very long waves the criterion is fairly well fulfilled. The relaxation time is determined by diabatic and subgrid phenomena which, with the exception of latent heat release, have a time scale of several days. This is somewhat larger than the time required by a typical particle to travel an appreciable fraction of the earth radius. This criterion then is marginally satisfied especially in places where there exists important release of heat by condensation. 6. 1. 1 Exchange Properties of Eddies If G represents a quasi-conservative property for isobaric motion, like potential vorticity, over periods of time r such that T < < T it is possible to define a field G(X,P, p) with a zonal mean Gz(f,p) z Bretherton (1966 a) has shown that the flux of property G is given by (G'v) a2 d Z aP [tb9 2 (6.1.2) defining the eddy coefficients as K= 2 d 1 2, K ~[ a K -.z (6.1. 3) where j denotes the south-north displacement of a parcel starting from latitude: po. Relation (6. 1. 2) was first found by Taylor (1921) and is valid under the requirement

137 62G _ z 2.. (6.1. 4) If G presents typical variations across distances of the order of the earth radius this inequality gives <<( 13000 km. It includes the cyclonic waves but is marginal with respect the ultra long waves. Equation (6. 1. 3) allows a direct computation of the eddy exchange coefficient when some tracer provides knowledge of the meridional displacements Yl (using results from the GHOST Project, Wooldridge and Reiter (1970) have computed values for the Southern Hemisphere at 200 mb). This equation also points out the universal character of the coefficient for quasi-conservative properties that comply with condition (6. 1. 4). This relation for K does not depend on the random character of the eddy motion nor its small amplitude and is therefore valid for stationary or finite amplitude disturbances. The small perturbation theory permits some additional interpretation. Identifying G with the potential vorticity Q, the transport in terms of zonal mean and its deviations is BQ (vQ') ( Y )z'' aTo Consider the mean zonal flow at an isobaric level dx dt = U(QP)

138 where dx = a cosYd?, and a small perturbation superposed on it dy - v(Xp, t) = v(t) cos a(x-c t) dt o where dy= a dp. Let us compute the trajectory of a particle that starts at y(O) = yo x(o) 0. To the lowest order x= U(Qo) and dt PO Cos [ U( ct] hence V y o + - sin (U -c )t o a o(U -c o where Uo= U(Po) and v v(^). Then 0 0 a "- y - Yo v Co) sin [-(U o-Co)t] The refore Uo o sin 2a - - x 2 o0 a(vi )z o 2a(Uo o z a ) = L 2a(Uo-C0) z (6. 1. 5) In general co is a complex eigenvalue that varies with latitude independently of Uo. If c is real the zonal mean of the bracket vanishes, since the inner quantity is harmonic in x ( or 0). Therefore a neutral wave is unable to support exchange unless UO = c0. Due to the fact that in most cases this condition will be satisfied locally this mechanism will be operative in narrow latitude zones where

139 Z"> d1 2 2 z so that a parcel is caught in the wave and taken away from its original latitude with the meridional displacement v increasing monotonically with time. If c0 is complex the bracket is not harmonic anymore and does not vanish when averaged. An unstable wave always support transfer, but in a more effective way when c U. The importance of the critical layer (where the basic flow of velocity equals the phase velocity of the wave) has been stressed by Dickinson (1969) who for an ensemble of statistically stationary, small amplitude, geostrophic waves, in adiabatic reversible flow with no sources of potential vorticity pointed to the critical layer mechanism as the only one capable of forcing a zonal flow by making the eddies transport potential vorticity. In this case the mean value of the rate of particle dispersion can be interpreted in terms of the power spectrum V (v) of the meridional eddy motion at the critical frequency, i. e. d 1 2 1 V dt (2 )z 2a2 rV(=UO) where the wedge brackets indicate a time mean and (v2)z) - 2 Diabatic effects and other sources of potential vorticity represent additional forcing mechanisms of a zonal flow perturbed by small amplitude, statistically stationary, geostrophic waves.

1.40 Finally we still have the finite amplitude waves that in the neutral and unstable cases will support transport processes. Relation (6. 1. 3) indicates that in all likelihood these waves are very effective agents. However, very little is known about them. If we ignore the ultra-long stationary waves and accept that finite amplitude eddies originate from small unstable perturbations, we have at our disposal the results from baroclinic instability theory that at least shed some light on the processes that generates them. Summarizing we have identified three possible causes for the eddy transfer: a critical layer mechanism, small amplitude unstable waves and finite amplitude waves. It is not easy to asses their relative importance. As a working hypothes is we will assume that baroclinic unstable waves are the responsible agents of transport. Because unstable waves in mature stages tend to present a smaller growth rate than small amplitude waves and because nonlinear effects eventually stop the development, information gained through linear theories has to be used with care. We will only try to determine from baroclinic instability studies the most relevant features of the mean zonal state that favor the generation of waves and its growth. 6. 1. 2 Generation of Eddies: Baroclinic Instability The fundamental cause for the existence of the large scale eddies is the differential solar heating and its immediate consequence: the meridional temperature gradient. The birth and initial growth of the disturbances embedded in such an environment has been studied extensively. However, in most of the cases the problem has been simplified to such an extent that although qualitative conclusions may still be valid, the magnitude of the results is subject to considerable

141 error. Very frequently the assumptions include: hydrostatic, quasigeostrophic, adiabatic and frictionless motion, the beta plane approximation, and a basic flow independent of the meridional coordinate with only small deviations from it permitted. Some of these simplifying hypotheses rest on observations and therefore do not detract from the applicability of the conclusions. Others have been introduced to overcome mathematical difficulties and their implications have to be considered with care. Such is the case with the smallness of deviations, introduced to gain linearity, beta plane approximation and basic flow independence of meridional coordinate, which makes the differential equations separable. The change from the beta plane approximation to spherical geometry has several consequences. First the stabilizing action of the Rossby parameter will change with latitude, but the effect of such a change is probably small because the order of magnitude of,8 remains invariable and our main concern are the most unstable waves and not the marginal cases. The spherical geometry has little effect on the equations of motion adding only small metric terms, but the differential operators have to be modified. In general we can expect little modification by removing this approximation. The smallness of deviations from the basic state is the key assumption that relieves the nonlinearity of the mathematical problem. At the same time, it limits the validity of results to those cases where waves have amplitudes so small that they do not modify the basic flow substantially nor do they interact among themselves. The most one can obtain from linear baroclinic theory is an initial growth rate and the identification of those features of the basic flow which are most important for eddy motion generation.

142 A basic flow independent of the meridional coordinate gives results quantitatively different from those where both vertical and meridional shear are included. Pedlosky (1964 b) has proved that omission of horizontal shear overestimates the growth rate and phase speed of disturbances, but does not change the minimum vertical wind shear necessary for instability. It changes also the spectrum of unstable waves by a shift towards longer waves. The most important difference for our purposes is that when only vertical shear is considered the parameters determining instability are the thermal wind, the 8 effect and static stability, (see for instance Miles, 1964); when only horizontal shear is included the important feature is the meridional variation of the absolute vorticity (Kuo 1949, Fj~rtoft 1951); but if both effects are considered, the meridional gradient of potential vorticity becomes the determining factor (Charney and Stern (1962), Pedlosky 1964 a, Bretherton 1966 a). The most relevant results for our particular case are those of Pedlosky (1964 a and b) who considers an incompressible two-level model in a beta channel of width 2f. His results have been adapted to our model in Appendix III. He establishes a necessary condition for instability, equation (A. III. 6), and sets an upper bound to the growth rate of the small eddies, equation (A. III. 11), both in terms of the meridional gradients of potential vorticity of the basic flow. He also concludes that eddies that behave like the atmospheric ones, i. e., eddies that converts zonal potential energy into eddy kinetic and then into zonal kinetic energy, occur in a basic flow with no extrema of potential vorticity in the upper level. This result, in the light of previous conclusions, implies a lower level with potential vorticity

143 gradient of opposite sign to that in the upper level if the necessary condition for instability is to be satisfied. This, in turn, can be interpreted in terms of the thermal wind, which has to be large enough to make such gradient in the lower layer negative. The upper bound on the growth rate is proportional to the aQ. maximum value of the product ui a ~. Since potential vorticity i a o% is very uniform for the kind of geostrophic flow under considerations (Phillips, 1963) this maximum is determined by the sign of the derivative and the wind profile, so that, with some accuracy, we can state t QI that u1 a~ z evaluated near the jet stream will determine the growth rate, 6. 2 RECENT MODELS OF THE HEAT TRANSPORT Lately two parameterizations of the heat transport have been published by Green (1970) and Stone (1972 b). The first used perturbation solutions to obtain a normalized vertical profile of the exchange coefficient. The magnitude of the coefficient was obtained for the vertically integrated value by assuming that one half of the typical zonal available potential energy is converted into eddy kinetic energy, that the meridional and zonal eddy kinetic energies are equal,and that the correlation between the eddy meridional velocity and the eddy entropy when normalized to the product of the root mean square value of these two quantities is constant. He obtains a relation (K) ( - (1 - (62 (A 0)2 -"(~~~~~

144 where the wedge brackets represents a zonal-time-pressure mean, the subscript s a standard value and A 9 the typical tropospheric difference of potential temperature 0 across the baroclinic region. Stone (1972 b) used Eady's (1949) model in his analysis. This model assumes a Boussinesq, adiabatic, inviscid and hydrostatic fluid located between two horizontal planes rotating with an angular speed f/2. For the perturbation analysis a basic state consisting of a zonal flow with constant vertical shear and no horizontal variations is considered. Stone (1972 b ) expands the dependent variables in powers of a common amplitude and then considers the first and second order solutions which are necessary to evaluate the eddy fluxes. The amplitude of these fluxes cannot be determined from the linearized theory where it is assumed to be small but increasing exponentially with time. To assign a magnitude it is assumed that the basic state can be taken as the mean state of the real atmosphere, thus making identical the constant Richardson number of the Eady's model with that obtained from mean values of the temperature field in the meridional plane. The next assumption is that the perturbations attain the same vertically integrated amplitude as the basic flow. This condition determines the amplitudes and provides, for a large Richardson number, the relation: r-2;5 _0 (6.2.2) where the wedge brackets represent now the mean value on the meridional plane and the bar a time-zonal mean.

145 Equations (6. 2. 1) and (6. 2. 2) agree in the linear dependence of the eddy coefficient on the meridional gradient of temperature but differ with respect to the static stability. Because the static stability shows very small seasonal changes it is difficult to decide which dependence is better. Using a relation similar to (6. 2. 2) for the vertical transfer of heat Stone (1973 b) has concluded that assuming a constant value for the static stability is not a bad approximation when studying climatic changes. Therefore, the difference between (6. 2. 1) and (6. 2. 2) is of little consequence. 6.3 A PARAMETERIZATION FOR MOMENTUM TRANSPORT In previous sections we have presented some arguments and results which suggest a way to parameterize the meridional eddy flux of momentum in an indirect fashion. The first conclusion is that it is easier to model the eddy transfer of a quantity the more conservative it is. In spite of the fact that this is a direct consequence of the nature of the mechanism of transfer it is only in Green's (1970) paper that the use of potential vorticity is indicated as a way to parameterize the momentum transport. A second idea derived from Bretherton's (1966 a) formulae for the eddy coefficient, is that any conservative quantity provides the same information if it satisfies the relation (6. 1. 4). This restriction indicates that some difficulty can be found in the representation of the ultra-long waves. A third and very suggestive piece of information comes from the baroclinic instability theory: the decisive parameter, when both vertical and horizontal shear are included in the basic state, is the meridional gradient of potential vorticity. Thus the necessary condition

146 for instability is a generalization of the criteria for barotropic and baroclinic instabilities, see Appendix III. Finally we have the results of Green (1970) and Stone (1972 b) for the heat flux. A conclusion can be drawn from them due to the universal character of the eddy coefficient. Because the transfer coefficient is the same for any quasi-conservative quantity the heat transport formula obtained by the aforementioned authors provide an estimate for the transfer coefficient of potential vorticity by cyclonic waves, and can be used in the parameterizations of momentum transport. Their results also suggest a more general parameterization. The eddycoefficient expressions (6. 2. 1) and (6. 2. 2) show that the main dependence is on the thermal wind in terms of which the stability criterion for baroclinic waves is expressed. Following the generalizations obtained by Pedlosky (1964 a) for the instability criterion it seems natural to think that the result of Stone (1972 b), who ignored the effect of variations of the Coriolis parameter, can be extended to make the eddy coefficient proportional to the meridional gradient of potential vorticity. Such a refinement may not be important in winter when most of the contribution to the gradient come from baroclinic effects but should become important in summer as the barotropic effects get at least as important as the baroclinic. Let us then postulate as a more suitable parameterization for the potential vorticity exchange coefficients an expression like K(P, t) = f(p,t) t ~ or |,-Qz( d (6 31)

147 where the interval ( 41, P2 ) encompasses the bulk of the baroclinic region, f ( ) is a shape factor and a a proportionality constant to be fixed by numerical experimentation or from observational studies. The shape factor f( t) was introduced because the evidence provided by Sela and Wiin-Nielsen (1971) indicates that the maximum of the exchange coefficients is much sharper than those of the mean zonal variables. This selectivity has been attributed to the critical layer mechanism (Green, 1970). Within the baroclinic region f ( ) is normalized to integrate to one, but corrections to the coefficients are applied in the tropical regions, where potential vorticities are not conservative quantities, and the very high latitudes, where the zonal mean losses statistical meaning. Considering the two-level model we are going to define one exchange coefficient for each level: Ki(P, t) f( Q, t) I i a i=, J_ 3 i= 1,3 (6.3.2) The exchange coefficient K2 will be expressed as a weighted mean of K1 and K3: K2(l, t) - f(P, t)O 0. 3 1 - + 0. 7 3 a dip The bias towards level 3 is justified by the evaluations of Wiin-Nielsen and Sela (loc. cit. ) and by the vertical asymmetry due to the location of the steering level (Green, 1970).

148 The shape factor f( ) is taken at every time step from a Beta distribution, 3 (xi/n, r), where the parameter n was fixed according to the observed profiles of K1 and K3, and r is evaluated from the first moment of the distribution of l I aQlz/ / with respect to the pole and considering contributions from 30 to 90 degrees of latitude. This last feature allows for changes in the location of the maximum. The independent variable in the Beta distribution is defined by x. = 3 1 3 so that 0 < xi 1. In the tropical regions, where we have no quasi-conservative quantity evaluated from the model, we used constant values for the exchange coefficients very much like those of Table 5. 5. 1. In high latitudes a linearly decreasing profile was used. A detailed description of these modification as well as numerical values for all constants is included in Appendix IV. In summary, the time dependent coefficients change in magnitude according to the mean gradient of potential vorticity in the baroclinic region and the geographical location of their maxima varies to a lesser extent because of the determination of the parameter r at each time step. This parameterization was introduced in the model described in Chapter 2 and a two year simulation was run for the Northern Hemisphere. 6.4 RESULTS OF A NONLINEAR SIMULATION The annual variation of the exchange coefficients is depicted in Fig. 6. 4. 1. The profiles shown are instantaneous values corresponding to the 15th day of each month. It is seen that the seasonal variation of the three coefficients is quite important reaching very low values in summertime. In particular K3 becomes so small that

Oas/ LU 901 S Iufn 9J idqdstluaq uxJq -aoN aq Jo uoT eInuTls aeau.Tuoi aq ui X'qU L puE'qD og I)'qc g Le suatOIJJao o appa jo SUOljleXeA TlUOsaS' I' 9 aanTzi Hanii ivi dN 09 OE 0 01~0 ~0 woa / O I A, i I I^ -wI9 3anll Vi dN O 09 O 0 qel i/ oj sdN 09 o, WoJsefloeA o....'. /;, -I,lnr -- 1! A1 9 Aienur ~ ^..-/ 9 9 z ]anll lvi dN 09 0o 0 w;.?,lo:, =o 0.. 0.... 6 T GM 1

150 a substantial discontinuity is produced in the transition from the tropical fixed values. Even more important, the normal size sequence is altered and K < K2 < K1 for a short period of time (the month of July). On the other hand the peak of the distribution shows a seasonal change of location: a ten-degrees migration poleward in summer. It is evident from the graphs that the very low values have a short lifetime and that during more than six months the peak is quite high. Comparing with the diagnostic values shown in Table (5. 5. 1) these are very similar to winter magnitudes for K2 and K3 but for K the maximum is shifted polewards by about 10 degrees and the magnitudes are about an annual mean value. The thermal field, Figs. (6. 4. 2) and (6. 4. 3), shows little change with respect, to that obtained using time independent exchange coefficients at 50 cb. The seasonal amplitude in high latitudes is smaller and there is aconcentration of baroclinicity in winter between 20 and 45 degrees of latitude. The latitudinal gradients for the extreme months, Figs. (6. 4. 4) and (6. 4. 5) show some improvement at both levels with a better simulation in summer and near the surface level. The zonal flow at 25 cb is shown in Fig. (6. 4. 6) which represents a substantial improvement over that obtained with the diagnostic coefficients. In January the maximum zonal wind is now correctly located (at 30 N) and its magnitude is slightly overestimated. There is also a five degree migration of the maximum towards higher latitudes in summer which is reasonable for a model in which the tropical circulation is underestimated. In winter a secondary maximum appeared near 75 N which is not present in tropospheric observations of the Northern Hemisphere. In summer maxima is correctly estimated in magnitude.

151 NP r NP 220 / / 230 40 60 24 \250 25 30 270-o EQ J F M A M J J A S O N D Figure 6. 4. 2 Same as Figure 4. 1. 1 for nonlinear simulation., NP -2402 60-/2 \60 30 9 270 300 EQ 1 I I I I I J F M A M J J A S O N D Figure 6. 4. 3 Same as Figure 4. 1.2 for nonlinear simulation.

152 OK ^o0 0 230 00 ~ ~ 0 ~ O 30 60 NP LAT ITUDE Figure 6.4.4 Same as Figure 4. 1. 3 for nonlinear simulation. 2 8 0 240 o~.~^ July 250 ~"~""< —,,..o o o 230.-..I I 0 30 60 NP LATITUDE Figure 6. 4. 5 Same as Figure 4. 1. 4 for nonlinear simulation. LATITUDE Figure 6.4.5 Same as Figure 4.1.4 ~~~~Vfornniersmtion

153 NP 10 ( - Q (ftft 60 I 3 ~0~0 Figure 6.4. 6 Same as Figure 4. 1. 5 for nonlinear simulation. Velocity correction 1.2 m/sec. 30 j EQ J F M A M J J A S 0 N D Figure 6.4. 7 Same as Figure 4. 1. 6 for nonlinear simulation. Velocity correction 1.2 m/sec.

154 The flow at 75 cb, Fig. (6. 4. 7), appears too high in winter but in summer is correctly simulated. There has been improvement also in the momentum balance since there are some westerlies in July between 20~N and 50~N. The field of vertical motion, Fig. (6. 4. 8), is now subject to stronger seasonal changes of intensity and structure. The three cell scheme is well defined in the winter time. In summer there are also three cells but the polar front ascent has been replaced by a narrow region about 450 with weak vertical motion. In the transition months there exist three regions of ascent although the important ones are still those associated with the polar front and the tropical convection. This irregular behaviour of the vertical motion is present in the observations reported by Qort and Rasmusson (1971). The heat transport shown in Fig. (6. 4. 9) presents strong seasonal fluctuations mainly due to a lower value in the summer months, although the winter maximum is also smaller than in the linear simulations. Considering the transient eddies transport only, the winter and summer maximum are about twice the observed values, indicating that the value of K2 is too large, but that seasonal fluctuations in magnitude are simulated fairly well. A comparison with results obtained using diagnostic values of the exchange coefficients is not straight forward because they include the effects of the ultra-long quasi-stationary waves, while the time dependent formulation might not include their effects. The angular momentum transport shows the effects of the global imbalance, although there is improvement. For instance the zonal velocity correction necessary to attain a zero transfer of angular

155 NP 360 is- > 30 J F M A M J J A S N D Figure 6. 4. 8 Same as Figure 4. 1. 7 for nonlinear simulation.

156 NP 60 30 EQ J F M A M J J A S O N D Figure 6. 4. 9 Same as Figure 4. 1. 8 for nonlinear simulation. 20 Jan 10 60 ~ r 30 \ NP \/July \ LATITUDE Figure 6. 4. 10 Same as Figure 4. 2. 9 for non linear simulation.

157 momentum between the earth and the atmosphere is 1. 2 m/sec. Fig. (6. 4. 10) shows schematically the variations of the poleward angular momentum transport, corrected for zero crossequator flux, for the months of January and July. The energetic evolution of the model is presented in Figs. (6. 4. 11) to (6. 4. 13). The seasonal fluctuations of zonal available potential and kinetic energies is only slightly different from their counterparts of the linear model. The time lag between generation and dissipation is about 24 days and their time evolution is again very similar to former results. The conversions show some change. First C( A A ) which is proportional to the heat transport presents values which are about one half of the corresponding results of the linear simulation. 2 The difference between the observed value of 3 watts/m2 and the simulation is due to the absence of a mean meridional circulation and the exclusion of very long waves to some extent. The conversion between the eddy and zonal kinetic energies shows a small positive value except in the fall when a barotropic behaviour of the eddies dominates. The small value of C(KE, Kz) can be blamed on the misrepresentation of the ultra-long waves which make up for a large fraction of the momentum transport in winter (Wiin-Nielsen et. al., 1963). The conversion between the two zonal energies is small and almost invariable throughout the year. The large positive value obtained is due to the vertical velocity field which in summer gives a large correlation with temperature.

158 J F M A M J J A S O N D Figure 6.4. 11 Same as Figure 4. 1.10 for nonlinear simulation. 4 J F M A M J J A S O N D Figure 6. 4. 12 Same as Figure 4. 1,11 for nonlinear simulation.

159 4 2 C (AZPAE/ 0 C(KEKZ) I I I I.. I..1. I..I I.. I,. J F M A M J J A S O N D Figure 6.4. 13 Same as Figure 4.1.12 for nonlinear simulation. 1.88 AZ 0,32 z 0.44 ~G 35396 1116 D G D 1.33 0.16 AE KE Figure 6.4. 14 Same as Figure 4. 1.13 for nonlinear simulation.

CHAPTER VII CONCLUSIONS This attempt to simulate the axisymmetric part of the general circulation of the atmosphere is based on a hemispheric twolevel quasi-geostrophic model. The eddy fluxes have been parameterized in terms of quasiconservative quantities only,through a diagnostic and a time dependent formulation. The model itself, very simple and inexpensive to operate, has several shortcomings which hamper the validation by comparison with observations. Among these the most important is due to the quasi-geostrophic character of the governing equations. Because of this feature the tropical regime is poorly represented. By the same reason the insertion of an artificial boundary at the equator is not important. When fluxes are parameterized in terms of diagnostic coefficients the transports by the very long waves is present in the eddy coefficients. This is not a complete representation of their effects, but a partial one. When the time dependent formulation of the coefficients is used, the influence of the ultra long waves is excluded to a larger extent. The parameterizations of the eddy fluxes of heat and potential vorticities, both considered as quasiconservative quantities in terms of unique exchange coefficients,determine their vertical variation through the global constraints of the angular momentum transport. Therefore, the exchange coefficients at the three active levels of the model are not independent of each other. Because of the 160

161 misrepresentation of the tropical mechanisms of transport, which cover about one half of the earth, it is safer to ignore the integral constraints for the sake of those latitudes properly represented in the model. An additional fact that detracts from the use of these integral constraints is the nonconservative character of the thermal streamfunction. This feature of the thermal streamfunction led us to avoid a formulation of K2 in terms of its meridional gradient in the time dependent coefficients case. The model is very sensitive with respect to the streamfunction constant fo and the boundary layer friction coefficient 8. Some of the results, like the zonal available potential energy, has to be considered within somewhat large uncertainties compatible with the model. Increasingthe empirical coefficient of skin friction ~ allows an improvement of the zonal flow in summertime. Hence, some defects, like polar easterlies extending too far away from the pole, are features that can be corrected by an increase of E. This is an inference which is substantiated by the improvement of the zonal flow for the Southern Hemisphere case. The results of the study have to be considered within the frame that the restrictions and flexibilities just mentioned are able to define. The first conclusion that can be drawn is that many of the unrealistic features obtained in former studies can be relieved by the use of a more complete heating function. In particular, we must include the release of latent heat which determines to a large extent the meridional variations of the thermal forcing. This modification produced a much more realistic temperature field and quantities

162 associated with it are improved correspondingly. As a consequence, the approach of casting the governing equations in terms of quasi-conservative quantities and through them parameterizing the eddy flux of nonconservative quantities, such as momentum, seems to be quite successful. Through the use of diagnostic time independent eddy coefficients it is possible to reproduce the observed thermal and flow fields and most of their seasonal evolutions. In some aspects the energetics of the hemispheric atmosphere is simulated quite well in the annual mean as well as in the seasonal variations. A second conclusion relates to the variation of the eddy exchange coefficients in space. A characteristic feature of the results by Wiin-Nielsen and Sela (1971) is the sharp maximum of the exchange coefficients between latitudes 50 and 70. Such a peak appears also in the results of Clapp (1970) for the heat transport by cyclonic eddies. The omission of the latitude dependence in the mixing coefficients results in a zonal flow where the maximum wind moves to high latitudes (70~ N). This omission is unreasonable. Baroclinicity and cyclonic activity are characteristically concentrated in midlatitudes. The pressure dependence of the eddy exchange coefficients is also of primary importance as theory and observation show. Therefore, we must conclude that the use of austauch coefficients constant in time and space is not suitable for our purposes. The time dependent formulation of the eddy coefficients introduced produces so mewhat stronger seasonal fluctuations which in most of the cases are supported by observations. However the implicit exclusion of the ultra long waves makes it difficult to verify the results. This approach restores the nonlinearity of the

163 model and for that reason itself it is an improvement for climate simulation purposes. We feel that although the postulated parameterization is far from being an ideal one, it represents a promising improvement. The results presented here cannot be taken as a conclusive proof, but show that the mechanism introduced does work in the right sense in the seasonal fluctuations and is able to reproduce observed features otherwise excluded. In the experiments reported here, the thermal field behaved in a comparatively stable way with respect to the field of zonal motion. The energetic cycle is very stable. This is probably a result of their being integrated quantities, but not only the annual mean values share this stability. The phase relations of the annual cycle also appeared to change little from one experiment to the other. Finally, we would like to suggest some avenues that, in our judgement, are worth further exploration. In relation to the thermal forcing, there has been an appreciable modification of the derived values of the surface albedo in low latitudes during the last few years. It seems desirable to modify the thermal forcing so as to include this new information. However, this is not a straightforward change because other fluxes may have to be altered at the same and/or other latitudes. The Southern Hemisphere seems to be worth additional efforts directed toward an improvement of the heating function. As far as the model is concerned, the most important limitation is in its quasi-geostrophic nature. Hence, an obvious suggestion is the use of a model based on the primitive equations.

164 In such a case, the model can be made global and the eddy fluxes expressed in terms of quasi-conservative quantities such as potential vorticity. Retaining a quasi-geostrophic formulation there is one modification that can produce improvement. Because the vertical divergence of radiative fluxes reaches higher values in the lower troposphere and condensation processes are also concentrated there,and because the heat transport is biased towards the same levels, it seems natural to modify the vertical structure of the model as to specify temperature at levels 1 and 3 as is made with vorticities. This kind of model suggested by Lorenz (1957) allows in addition changes in static stability and avoids the necessity of using a constant fo. Therefore, the equator can be crossed without a discontinuity in the zonal flow. However, no significant improvement in the low latitudes can be expected as long as quasi-geostrophy is present.

APPENDIX I FINITE DIFFERENCES EQUATIONS A.. I1 Prognostic Equations Potential vorticities are obtained from two equations of the form aQ 1 Qa aQ at a cosa ap(K os-p- + S i=1,3 (A I. 1 1) a2 CO O These relations can be considered as two separate parabolic partial differential equations with source terms that couple them: S= -X H2 - 2A~t I -2Az Tz (A. I.1.2) 3 H2z + 2A3 <3 1i ) S3 2 zTz (2- 3z, 2 1 z We will deal then with the general form (A. 1. I. 1) omitting the subscript until it becomes necessary. Boundary conditions for the solution of (A. I. 1. 1) have already been determined in Section 2. 2. The initial condition for our experiments will be a simple state, say Qiz = f + 6 q2 Tz Qiz Tz where YTz is independent of latitude. This initial condition and non vanishing values for Ki imply a state such that there is no mean relative motion, but well developed cyclonic random action still exists. 165

166 It is certainly contradictory in itself because cyclonic activity derives its energy from meridional temperature differences and this in term implies a vertical wind shear on a rotating earth. However, we will not be concerned with the spin up problem and will concentrate the attention on the phenomena of the steady regime that will show little trace of the initial condition from which they evolved. From a mathematical point of view the initial condition can be arbitrary. Given boundary conditions at (= 0 and = -- (for the NH) and an initial condition the problem stated by (A. I. 1. 1) has a unique solution, Richtmyer and Morton (1957). The stability of the computation is insured by the use of the Crank-Nicholson scheme. The source terms (A.I. 1.2) however, contains a part that cannot be expressed in terms of the corresponding potential vorticity Q and the scheme cannot be maintained for these terms. They will be used as of the last time step computed, hence they lag the diffusion term in (A. 1. 1.1) by half time step (6 hours). This must be of minor importance since the typical time scale we are concerned with is of the order of three months and source mechanisms have a characteristic period considerably larger than the time step used. With a time step At equal to 12 h. and a latitude grid length eA of 5 degrees, using a superscript n for the time level and a subscript m for the grid point, with m=l for the equator and m=19 for the pole, we can write the difference form corresponding to (A.. I. 1): n+n n - 2K (Q n + Q) + Km-/ ( Q Q + S r m -m /m -C5 p "m-l -m n

167 for 1 < rn ( 19 and n > 0. Rearranging -. - mtt cosg p + Qm [1 + kp 2 K n+l [K c SP ] Qn+l K Qm~+l 2 (a A?). 2 M+Y.1+l CO+ (a A P2 m + -Qnt c^cos m At COS2m+ Q+ K m- 1/Z Q n K M m-1 2(a a2)2 Km-cOSm Qm+ 2(a0 m )2 m+ ]cOSpm -at m t cos r n n m + Sn +- Q [a K +Q K m (aAY )2 m M-2 2 mrCos m (A.I. 1..4) Define + At cost C K ___ m 2(a? )2 m+ cosp 2C = C +C (A I. 1.5) mm rn m Equation (A. I 1. 4) can now be written as + n+1 n+1 n+l n - C Q + (1 + 2C ) Q -C Qn+ Bn m m+l m m -m m m- m (A.I. 1.6) where Sn n n n B =S +C Q + (1 -2C) Q + C Qmm m m m+1 m m m m-1 and m= 2, 3,..., 18 and n 0. The limit form of equation (A. 1. 1 1) at the pole were (2. 2. 12) and (2. 2. 13) from where it is apparent that we have to provide a onesided finite difference form for the first and second derivatives of

168 potential vorticity. It is in this process that the boundary condition (2. 2. 23) imposed at the end of Section 2. 2 is useful. At this point we note that all numerical schemes adopted thus far, have truncation errors of second order. The first and second derivatives have been approximated by central differences and for quadrature the trapezoidal rule has been used. To be consistent we must require the same order forms for the boundary equations. Note also that equation (A.I. 1.6) has a tridiagonal matrix which can be easily decomposed into an upper and a strictly lower triangular matrix. We could satisfy the truncation order requirement without imposing zero values of first and third derivative of Q, but the matrix of the system (A.I.. 6) would not be tridiagonal anymore. We chose to satisfy the second order truncation error, keep a tridiagonal matrix and impose (2.2.23). With aQ/ a = 3 3Q/33 = 0 we can approximate 2 2Q 2 ~2 (p)2 = )2 (Q18 -Q19 )+ (^t) and the finite difference equation for the pole becomes n+l n Q19 - 19 2 n n+ n n+l n+ n l "(a= K19 (18 + Q18 Q19 Q19 ) + S19 at (aul) or 2 t n+1 21t n+1 ( 2 018l K9 + (1 + 2 K19 Q1 9 (a^ (a 19 (A.I. 1.7) 2 t 2 t 2Ast n 2t n (a 2 Q18 K19.+(1Q r + S.

169 So if we define At C19 ( 2 K19 (aMASO) n n + n 19 = S9 + (1 - 2C19) Q19 + 2C19 Q18 We can write (A. I. 17) as n+1 n+1 n (1+2c19)Q 19 2c19 18 B19 (A.I. 1.8) Likewise for the equator (2. 2. 20) and (2. 2. 21) give for n = 1, with ^t C = K -1 2 (aA%)2 1 n n n n B S1 + (1 - 2C1 )Q1 + 2C1 Q2 the equation (1 + 2C1 ) Q1 - 2C1 Qn+2 = B 1 2 1 ^ (A.I. 1.9) We can write a matrix relation for the equation (A. I. 1. 1.) Ax=b where the transposes xT and bTof x and b respectively are.-2T n+l n+l n+l B n n Q1 Q9' Q19 a b B 2 B... B1I

170 anc the matrix A is 1+2C1 -2C1 -C2 1+2C2 -C -C3 1+2C3 -C ^3 ^^3 3. A C18 1+2C18 -C -2C19 1+2C19 Considering (A.I. 1.5) we can see that this matrix is strictly diagonal dominant and therefore nonsingular. A.I.2 Helmholtz Equation for Thermal Streamfunction T. The equation for fT is (2. 2. 26), valid for all inner grid points, which in finite difference form becomes c ms m- /2 n 2 2 n + (al)2cos T. m-l L (a^)2 + q n.m cos Y,m+ 1 n) +~ TYT, m+l 2- Q1, m 3, m (aA) cos' (A. I. 2. 1) where n 0, m=2,3, 3..., 18.

171 The limit form of the equation at the pole is (A.. I. 2. 10) and its finite difference form is 4 4 [(a)2 q T, 19 (a+ )2 T, 18 2 - 1, 19 3, 19 I. 2.2) For the equator we had equation (2. 2. 19), Since we know that the boundary condition /3 aP = 0 chooses an even VT with respect to the equator we can write 2 = )2 (W 2T, 2 T,1) [() and the finite difference form of (2. 2. 19) is 2 2 n2 + n 1 Qn _ n " ~ ^2' q + 2 2 ~" J, (a (A... 2.3) In matrix form Ax=b where,T I fn yn Vn x~' _ T, 1 T, 2' T, 19, T /n 1n 1 n I n \] - 2 (Q1, 1 i) Q 2Q, 2 3,2 2.2(Q1,19 Q193,1) and if we now define + I cosmP,/ Dm (a p )2 cospm D =(a b)2 D =-D+ +D m m m

172 We can write 2 -(2D+q2) 2D D -(2D2+q) D 2 + D3 -(2D3+q2) D3 D18 -(2D18-q ) D18 -(4 D+q2 The strict diagonal dominance of the tridiagonal matrix is evident, hence its inverse exists.

APPENDIX II SPLINE INTERPOLATION The divided differences of a function f(x), k times differentiable in [a b] given at a set of discrete points xi, distinct or not, can be defined as f(xi) k 0 f[x x fi+] E x ] ~f~~i xi'., ]|Xi+l' * * i+k iv * xi+k=-1 -^ i 1 X i+k x f(k)! k>O, Xi=xi+l=.= Xi+k k1 ^ i i -+1 i+k (A. II. 1) Then k ipk(X) = l f x.. xi+j (x is the unique polynomial of degree less or equal to k which interpolates f(x) at xi, Xi+.,., xi+k. The spline scheme interpolates a function f(x), given at N+1 points xi, i=l, 2,..., (N+1), by a cubic piecewise polynomial gr(x)=Pi(x), where Pi(x) is a cubic polynomial such that Pi(x.) f(x) Pi(Xl ) - f(x ) 1,, N i i+l i+l The other two conditions needed to fully determine Pi(x) will be chosen as to make Pi(x) agree with f(x) not only in value, but also in slope at Xi and xi+1. Therefore, g.r(x) is continuous and differentiable. 1 73

174 Using the divided differences (A. II. 1) we can write P(x) = f(x) + f [xi xi] (x-x) + f, xi xVi~ (x-x)2 + + f x, xi? X. X. 1 1 i+' i+l (x-x)2 (x-x which using the identity x-xi+ = ( i) + (i+1 can be transformed into P.() c, i + c2, i(x-x) + c3, i(x-xi + c4, i(x-x (A.)32) where c, i= f(xi) f 1 1^^i 2, i= f'(xi) = Si 1, 1 i 1'3 i= (AYi)( {3f[xiix+,1] -2i - i+l} 4, i= (Yi) ( {2f[i Xi+. + + Si + l Yi =Yi+l - Yi. 3) These c.. are such that g (x) defined by g (x) Pi(x) if XE [xi, Xi+1] 1 i + 4 =- c;i (x-x)'i^ -= 1, 2,..., N j=l satisifies g (x.)- f. 7 ij = 1,2..., (N+1) g(xj) = s. n 3 3

175 So if we know xi, fi and s. for i=l, 2,.., (N+1) we can find g (x), which is called the cubic Hermite interpolant. However, normally we have no information on the derivatives s., and we must choose some method to estimate them. In the spline scheme we determine them so that g (x) is twice differentiable, i. e., Pi+ (x) = P (x) i= 2 3,..., N i+1 1 Using (A. II. 2) this condition gives OYi Si-1 + 2(hyi -BYi-l) Si + 4yi-lsi+l Afi_ fif 3 b'Yi -L Yi + hYi i t i=2, 3,.., N These are N-1 relations among N+1 unknowns si, i=l, 2,...,(N+1). If we pick the end points derivatives, sl and SN+1, we can in principle determine s2, s3,..., N by (A. II. 4). With xi, fi., and si known the relations (A. II. 3) serve to evaluate cj i and (A. II. 2) to obtain interpolated values. Note that (A. II. 2) is a picewise cubic polynomial continuous and twice differentiable at the data points.

APPENDIX III INSTABILITY OF A TWO LEVEL MODEL A. III. 1 Basic Relations In this appendix we will consider the stability of a zonal current which varies with latitude and pressure, so that a perturbation can grow at expense of the zonal available potential or the zonal kinetic energy. This will be done specifically for the two-level standard model, which we are using, and in terms of potential vorticities whenever possible. This analysis follows Pedolsky's (1964a) presentation. The governing equations of the standard two level quasigeostrophic model for an adiabatic and frictionless flow can be written in terms of the potential vorticities (equations (2. 1. 16) and (2. 1. 17) ). t )y 3x +ax Dy O22 2 2 3 (A. III. 1.1) where d ={1 when i= 3 For simplicity we will work in a beta channel of width 2. Suppose we have a steady basic state which satisfies equation (A. III. 1. 1), is independent of x and varies with y and the level of the model, say, U1 (y) i-1 U. - U3(y) i=3 Let us introduce a perturbation to this basic flow p3 i== 3 176

177 The perturbation equation is 2 2 t + Ui x) 2 + 2 8 (43 (I)1).. 2 ax f + U3 U1) dy, Since the coefficients of perturbation quantities are independent of x and t we can assume that the perturbation is harmonic in x and t with wavenumber a =i Re i (y) ei(t)} (A, III. 1.2) Substitution of (A. III1.. 1.2) into the perturbation equation gives 2 (U. ) ( - (2i)+ -i ( - U + q (U - + 2 + - (Ui -c)68(3-1) 0 (A. III 1. 3) df where = dyNoting that the basic flow potential vorticities are 2 Qi = - Ui +f + q 8( 3 - 1 ) 1 2 3 1 we can write (A III.1. 3) 2 (U -c) ( 2 cri)+i (U -c)-((-) 1' ~ 1'~ ~ 1'~ 2'~ 1^~'~3)~ ^(A. III.1.4)

178 At the boundaries of the channel, y=~l, the perturbation meridional velocity must vanish for all x, hence i= ). Since both the governing equation (A. III. 1. 4) and the boundary conditions are homogeneous, solutions will exist for certain values of c=c ici. A. III. 2 Necessary Conditions for Instability Assume a growing disturbance (ci # O), so that division by (Ui-c) is permissible. Then multiply the first equation of (A. III. 1. 4) by 2 4'k/ q (U1-c), where asterisk denotes a complex conjugate quantity, * 2 and the second by 2 /I q (U3-c), and the two equations are integrate from y= 1 to y= -1. Using integration by parts and the boundary conditions one gets dy' -— ] y J91 i= e, 3L fY ~E (I|i 12 +a2iI2)+ - 2 2 (A. III. 2. 1) The imaginary part of this relation provides IdY U IPil c b) (A. III. 2.2) showing that a necessary condition for instability is that the potential vorticity gradients cannot have the same sign everywhere. Because 2 Q.1 2 ay'y 1 (A. III. 2.3)

179 this statement indicates that 1 is a stabilizing factor, the thermal wind UT= I(U1-U3) is a destabilizing influence, and the static stability 2 included in the denominator of q tend to make the flow stable. If we identify the basic flow with the zonal mean wind, observations show that its relative vorticity is unimportant compared with the other terms in the right hand side of (A. III. 2. 3), we have approximately Q. 2 1 -= + 8q2 U y +q UT (A.III. 2.4) and a necessary condition for instability is UT > 2 q a well known relation deduced by Phillips (1954), for a basic flow independent of y, by the solution of the eigenvalue problem. In the barotropic case (UT= 0) the necessary condition for instability is the same, but in terms of absolute vorticity. For a continuously stratified fluid Pedolsky (loc. cit. ) has shown that the condition reduces to the well known result of Kuo (1949) and Fj rtoft (1951). The real part of (A. III. 2. 1) gives, using (A. III. 2. 2) ~i2 Qi dy' EU c|2 i ~y fdY)Z (J2 +i a Jdy (I(2 + aI i12)+ 02 |1i 12.5 i= 1,3 (A. III. 2. 5) Hence the product of the zonal flow and the meridional gradient of potential vorticity must be somewhere positive in unstable flows,

180 A. III. 3 Growth Rate and Upper Bound Relation (A. III. 2. 5) allows the determination of an upper bound to the growth rate of the perturbation. Considering that IU -c2>, c2 j=1, 3 and that if g(y) is a well behaved function that vanishes at y=~ 1 2 Jg' dy > 2 dy J1' -2|41 (A. III. 3.1) we can deduce from (A, III. 2. 5) 9.1 B Qi c- 2 z max ( U. 41 i (A. III. 3.2) A. III. 4 A Semicircle Theorem Starting from Equation (A. III. 1. 3) and following Pedlosky (loc. cit. ) one can find that the results of this author are unchanged for our model, since the terms which could make a difference are lost in the simplification of the inequalities. The results are f <U Umin 2 2 r max 2 + 2(a +412 (A. III. 4. 1) and U +U 2 U -U U -Um U maxmin 2+ 2 max min 2 max Uin Cr i ) + 2 - ~ r 2 i2 ^2 ~ 2 2 Ir ^ 41 (A. III. 4. 2)

181 Ci max3~ Umax-Umi 2.~'/ 2 41z ~/2~2 or not we see that instability can be due to the vertical or horizontal shear. To make the flow stable it is necessary to have both zero. The beta influence results in unstable waves moving slower max than the minimum basic wind, but Uax is still an upper bound to theax phase speed. A. III. 5 Short Wave Cut-off The short wavelength cut-off of the unstable waves can also be expressed in terms of the potential vorticity gradient (Pedlosky, 1964a1 Bretherton, 1966h).

APPENDIX IV TIME DEPENDENT COEFFICIENTS Expressing the latitude? in degrees the expressions adopted for the three eddy exchange coefficients are K1 = f(~,t) / a dzI 12J d K3 = f( t) | 3 Q3z I 330 3 130 K2 f(,t) (0. 3 a. + 0. 7 a3 ) d The constants al and a3 were fixed by a combination of numerical experiments and the information provided by the diagnostic study of Wiin-Nielsen and Sela (1 971). The values used in the nonlinear simulation reported in 17 2 17 2 Chapter 6 were a1= 0. 25x10 m and a3=2. 5x10 m. The shape factor f(e) at a given time, was taken from a Beta distribution /(x/r, n). The parameter n was taken as time independent and its value was determined by adjusting Beta distributions to the diagnostic profiles determined by the aforementioned authors. Since the variations of n for different levels differed slightly a mean value of n=12 was taken. The adjustment was made considering only the values from 30 degrees to the pole and using the mean and variance of the distribution with respect to the pole. Then for the distribution r-1 n-r-1 f() (x/rn)1 (1x) for O-xxl 0 3xr ) 1(r, n-r) where x = (90-~ )/60 and 182

183 3(r, n-r)= / x (lx)nldx the mean and the variance are respectively r and n-r A1- n and A2 n(n+l) The magnitude of r was obtained at every time step from an evaluation of Al for the integrand of K2: 6 (0.,Z 3 + 0. 7 a31 o ) (li 3 18Y 13 r n _ -Q t 0. 31 lz a l zl 2. 3I c + 0.7 73I ) (9i i=6 Therefore 90 f(O) dP = 1 0 Inthe tropical regions the eddy exchange coefficients were modified in the following way: for 0. T 0 265} K= j=1, 2, 3 11 i,6^ K, i,i 3- 2 for 25 <i 83 Kli K1,i Kj i - K j m=2 K 1 In the high latitudes the coefficients were modified as follows: for 70 < 90. = K'j 1 23 15 < i 19 Kj,i i j 1, 2, 3 where -1 6 2 -1 K'= 0.2x106m2sec1- j=2 Kj' i 0. 2 x106 m2sec-1 j= 2 1.5xl0 m sec j=3 "' _= KT (90-) j,i, 15 20 = 2, 3 and j, i= Kj( pi)

REFERENCES Bretherton, F. P., 1966a: Critical Layer Instability in Baroclinic Flows, Quart. J. Roy. Meteor. Soc. Vol. 92, 325-334. Bretherton, F. P., 1966b: Baroclinic Instability and the Short Wavelength Cut-off in Terms of Potential Vorticity, Quart. J. Roy. Meteor. Soc., Vol. 92, 393-403. Bryan, K., 1969: Climate and the Ocean Circulation: the Ocean Model, Mon. Wea. Rev., Vol. 97, 806-827. Clapp, P., 1970: Parameterization of Macroscale Transient Heat Transport for Use in a Mean-Motion Model of the General Circulation, J. Appl. Meteor., Vol. 9, 554-563. Cressman, G. P., 1960: Improved Terrain Effects on Barotropic Forecasts, Mon: Wea. Rev., Vol. 88, 2-12 and 327-342. Charney, J., 1959: On the General Circulation of the Atmosphere, Rossby Memorial Volume: The Atmosphere and the Sea in Motion, Rockefeller and Oxford University Presses, New York, 178-193. Charney, J., 1960: Integration of the Primitive and Balance Equations, Proc. of the Int'l. Symp. on Numerical Weather Prediction, 7-13 November 1960, Tokyo, 131-152, (published March 1962). Charney, J. and M. E. Stern, 1962: On the Stability of Internal Baroclinic Jets in a Rotating Atmosphere, J. Atmos. Sci., Vol. 19, 159-172. Charney, J., 1971: Geostropic Turbulence, J. Atmos. Sci., Vol. 28, 1087 -1095. Davies, P., 1963: An Analysis of the Atmospheric Heat Budget, J. Atmos. Sci., Vol. 20, 5-22. 184

185 REFERENCES (continued) Dickinson, R. E., 1969: Theory of Planetary Wave Zonal Flow Interaction, J. Atmos. Sci., Vol. 26, 73-81. Eady, E. T., 1949: Long Waves and Cyclone Waves, Tellus, Vol. 1, 33 -52. Fjortoft, R., 1951: Stability Properties of Large-Scale Atmospheric Disturbances, Compendium of Meteorology, Amer. Meteor. Soc., Boston, 454-463. Green, J. S. A., 1970: Transfer Properties of the Large-Scale Eddies and the General Circulation of the Atmosphere, Quart. J. Roy. Meteor. Soc., Vol. 96, 157-185. Jenne, R., H. Van Loon, J. J. Taljaard, andH. L. Crutcher, 1968: Zonal Means of Climatological Analysis of the Southern Hemisphere, Nots, Vol. 17, 35-52. Krueger, A. F., J. S. Winston and D. A. Haines, 1965: Computations of Atmospheric Energy and its Transformation for the Northern Hemisphere for a Recent Five Year Period, Mon. Wea. Rev., Vol. 93, 227-238. Kuo, H., 1949: Dynamic Instability of Two-Dimensional Non-Divergent Flow in a Barotropic Atmosphere, J. Meteor., Vol. 6, 105 -122. Kurihara, Y., 1970: A Statistical-Dynamical Model of the General Circulation of the Atmosphere, J. Atmos. Sci., Vol. 27, 847-870. Laikhtman, D. L., et al., 1970: Problems in Dynamic Meteorology, World Meteorological Organization, WMO- N 261. TP. 146, Geneva, Switzerland, Edited by A. Wiin-Nielsen, 245, Appendix 15.

186 REFERENCES (continued) List, R., 1951: Smithsonian Meteorological Tables, Smithsonian Misc. Tables, Vol. 114, Washington. London, J., 1957: A Study of the Atmospheric Heat Balance, Report: New York University, Dept. of Meteor. and Ocean., College of Engineering. London, J. and T. Sasamori, 1971: Radiative Energy Budget of the Atmosphere, Space Res., XI, 639-649. Lorenz, E. N., 1955: Available Potential Energy and the Maintenance of the General Criculation, Tellus, VII, 157-167. Lorenz, E. N., 1960: Energy and Numerical Weather Prediction, Tellus, XII, 364-373. Lorenz, E. N., 1967: The Nature and Theory of the General Circulation of the Atmosphere, WMO, Geneva, 161 pp. Manabe, S., 1969: The Atmospheric Circulation and the Hydrology of the Earth's Surface, Mon. Wea. Rev., Vol. 97, 739-774. Miles, J. W., 1964: Baroclinic Instability of the Zonal Wind, Rev. of Geophy., Vol. 2, 155-176. Mbller, F., 1951: Vierteljahrsten des Niederschlags fur die ganze Erde, Petersmanns Geograph. Mitt., Vol. 95, 1-7. Neumann, G. and W. J. Pierson, Jr., 1966: Principles of Physical Oceanography, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 545 pp., Table 14.3, pg. 437. Newell, R. E., et al., 1969: The Energy Balance of the Global Atmosphere, The Global Circulation of the Atmosphere, Edited by G. A. Corby, Pub. by Roy. Meteor. Soc., London, Joint Conf. held 25-29 Aug. 1969 at the Roy. Soc. London.

187 REFERENCES (continued) Newton, C. W., 1972: Southern Hemisphere General Circulation in Relation to Global Energy and Momentum Balance Requirements, Met. Monographs, Vol. 13, Amer. Meteor. Soc., Boston, 215240. Oort, A., 1964: On Estimates of the Atmospheric Energy Cycle, Mon. Wea. Rev., Vol. 92, 483-493. Oort, A. and E. M. Rasmusson, 1971: Atmospheric Circulation Statistics, NOAA Professional Paper 5, U. S. Dept. of Commerce, Nat. Ocean. and Atmos. Admin., Rockville, Md., 323 pp. Palmen, E., 1955: On the Meridional Circulation in Low Latitudes of the Northern Hemisphere in Winter and the Associated Meridional and Vertical Flux of Angular Momentum, Soc. Scient. Fennica., Comm. Phys. Math., Vol. 17, 1-33. Palmen, E. and W. Newton, 1969: Atmospheric Circulation Systems, Academic Press, New York, 603 pp. Pedlosky, J., 1964a: The Stability of Currents in the Atmosphere and the Ocean: Part I, J. Atmos. Sci., Vol. 21, 201-219. Pedlosky, J., 1964b: The Stability of Currents in the Atmosphere and the Ocean: Part II, J. Atmos. Sci., Vol. 21, 342-353. Phillips, N. A., 1956: The General Circulation of the Atmosphere: A Numerical Experiment, Quart. J. Roy. Meteor. Soc., Vol. 82, 123-164. Phillips, N. A., 1963: Geostrophic Motion, Rev. of Geophy., Vol. 1, 123 -176.

188 REFERENCES (continued) Richtmyer, R. D. and K. W. Morton, 1957: Difference Methods for Initial Value Problems, 2nd Edition, Interscience Publishers, New York, 405 pp. Riehl, H., 1951: The North-East Trade of the Pacific Ocean, Quart. J. Roy. Meteor. Soc., Vol. 77, 598-626. Saltzman, B., 1967: On the Theory of the Mean Temperature of Earth Surface, Tellus, Vol. XIX, 219-229. Saltzman, B., 1968: Steady State Solutions for the Axially Symmetric Climate Variables, Pure Appl. Geophys., Vol. 69, 237-259. Saltzman, B. and A. Vernekar, 1971: An Equilibrium Solution for the Axially Symmetric Component of the Earth's Macroclimate, J. Geophys. Res., Vol. 76, 1498-1524. Sasamori, T., J. London and D. V. Hoyt, 1972: Radiation Budget of the Southern Hemisphere, Meteor. Mon., Vol. 13, N 35, Meteorology of the Southern Hemisphere, Edited by C. W. Newton, Published by the American Meteor. Soc., Boston. Sela, J. and A. Wiin-Nielsen, 1971: Simulation of the Atmospheric Annual Energy Cycle, Mon. Wea. Rev., Vol. 99, 460-468. Sellers, W. S., 1965: Physical Climatology, The University of Chicago Press, Chicago, 272 pp. Smagorinsky, J., 1963: General Circulation Experiments with the Primitive Equations: I the Basic Experiment, Mon. Wea. Rev., Vol. 91, 99-164. Smagorinsky, J., 1964: Some Aspects of the General Circulation, Quart. J. Roy. Meteor. Soc., Vol. 90, 1-14.

189 REFERENCES (continued) SMIC, 1971: Inadvertent Climate Modification, Report of the Study of Man's Impact on Climate, The MIT Press, Cambridge, Mass., 308 pp. Stone, P. H., 1972a: A Simplified Radiative-Dynamical Model for the Stability of Rotating Atmospheres, J. Atmos. Sci., Vol. 29, 405418. Stone, P. H., 1972b: On Non-Geostrophic Baroclinic Stability: Part III, The Momentum and Heat Transports, J. Atmos. Sci., Vol. 29, 419-426. Stone, P. H., 1973: The Effect of Large Scale Eddies on Climatic Change, J. Atmos. Sci., Vol. 30, 521-529. Taylor, G. I., 1921: Diffusion by Continuous Movements, Proc. Lond. Math. Soc., Vol. 20, 196-212. Van Loon, H., 1972: Wind in the Southern Hemisphere, Meteor. Mon., Vol. 13, Ed. by C. L. Newton, American Meteor. Soc., Boston. Vonder Haar, T. H. and V. E. Suomi, 1971: Measurements of the Earth's Radiation Budget from Satellites During a Five-Year Period, Part I: Extended Time and Space Means, J. Atmos. Sci., Vol. 28, 305-314. Vonder-Haar, T. and A. Oort, 1973: New Estimate of Annual Poleward Energy Transport by Northern Hemisphere Oceans, J. Phys. Ocean, Vol. 3, 169-172. Wiin-Nielsen, A., J. A. Brown and M. Drake, 1963: On Atmospheric Energy Conversion Between the Zonal Flow and the Eddies, Tellus, Vol. 15, 261-279.

190 REFERENCES (concluded) Wiin-Nielsen, A., 1967: On the Annual Variation and Spectral Distribution of Atmospheric Energy, Tellus, Vol. XIX, 540-559. Wiin-Nielsen, A., 1969: On Atmospheric Response to Large Scale Seasonal Forcing, Proc. of the WMO/IUGG Symp. on Numerical Wea. Prediction, Tokyo, Japan, Nov. 20-Dec. 4, 1968, Japan Meteor. Agency, Tokyo, lV-21-1V-27. Wiin-Nielsen, A., 1970: A Theoretical Study of the Annual Variation of the Atmospheric Energy, Tellus, Vol. 22, 1-16. Wiin-Nielsen, A. and J. Sela, 1971: On the Transport of QuasiGeostrophic Potential Vorticity, Mon. Wea. Rev., June, Vol. 99, 447 -459. Wooldridge, G. and Reiter, E., 1970: Large Scale Atmospheric Circulation Characteristics as Evident from GHOST Balloon Data, J. Atmos. Sci., Vol. 27, 183 194.

UNIVERSITY OF MICHIGAN 9015 032III 1946III IIIIIIIIIIII 3 9015 03127 1946