Three-dimensional parameterizations of the synoptic scale kinetic energy and momentum flux in the Earth ’ s atmosphere

We present a new set of statistical-dynamical equations (SDEs) which can accurately reproduce the threedimensional atmospheric fields of synoptic scale kinetic energy and momentum flux. The set of equations is closed by finding proper parameterizations for the vertical macroturbulent diffusion coefficient and ageostrophic terms. The equations have been implemented in a new SD atmosphere model, namedAeolus. We show that the synoptic scale kinetic energy and momentum fluxes generated by the model are in good agreement with empirical data, which were derived from bandpass-filtered ERA-40 data. In addition to present-day climate, the model is tested for substantially colder (last glacial maximum) and warmer (2 ×CO2) climates, and shown to be in agreement with general circulation model (GCM) results. With the derived equations, one can efficiently study the position and strength of storm tracks under different climate scenarios with calculation time a fraction of those of GCMs. This work prepares ground for the development of a new generation of fast Earth System Models of Intermediate Complexity which are able to perform multimillennia simulations in a reasonable time frame while appropriately accounting for the climatic effect of storm tracks.


Introduction
Transient atmospheric eddies play a crucial role in forming the Earth's climatic state and hence require accurate representation in climate models.The high and low pressure systems associated with individual eddies comprise much of the day-to-day weather variability and their ensembleaverage form the so-called atmospheric storm-tracks in the mid-latitudes (Blackmon, 1976;Chang et al., 2002b).In the Correspondence to: D. Coumou (coumou@pik-potsdam.de) mid-latitudes, synoptic eddies carry out the bulk of the poleward atmospheric heat-transport and contain roughly half of the total kinetic energy in the atmosphere (Lau and Oort, 1982;Barry et al., 2002).The interaction between ensembles of synoptic eddies and the basic flow plays a crucial role in the stability of the latter and is thus of paramount importance for the overall structure of atmospheric circulation (Charney, 1947;Eady, 1949;Held and Hoskins, 1985;Held, 1999).Furthermore, in regions of high synoptic activity, exchange of momentum between atmospheric eddies and the oceanic mixed layer notably contributes to the intensity of oceanic flow, thereby directly affecting the thermohaline circulation (White et al., 1980;Alexander, 1992a,b;Weaver et al., 1999).Storm tracks thus have a critical role in both the atmosphere and the climate system as a whole, and changes in them are hence likely to play an important role in shaping future climate change (Swanson, 2007).
The current availability of high resolution atmospheric data-sets over several decades, like ERA-40 (Uppala et al., 2006) and NCEP-NCAR (Kistler et al., 2001), allows quantification of the synoptic scale field in terms of its kinetic energy and heat, moisture and momentum fluxes (Petoukhov et al., 2008).Blackmon (1976) and Blackmon et al. (1977) showed that 2-7 days time filtering of meteorological datasets isolates the transient eddies with spatial scales of the order ∼1000 km, allowing to quantify storm tracks.This and other early classical studies (Lau et al., 1978;Oort, 1983) identified the two northern hemispheric regions of major synoptic activity, in the vicinity of the east coasts of the North American and Eurasian continents at about 40-50 • N. In these regions, sharp horizontal gradients of large scale temperature and wind-speed are observed, and hence both strong baroclinic and barotropic transient eddy forcing are generally thought to fuel mid-latitude storm tracks (Pedlosky, 1979).Several decades of observations since then have provided a consistent picture of synoptic scale activity, quantified the seasonal variability and revealed variability on Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.D. Coumou et al.: Synoptic scale kinetic energy and momentum flux interannual (e.g.due to El Niño; Alexander, 1992a,b) and decadal timescales (Chang et al., 2002b,a) State-of-the-art atmospheric general circulation models (AGCMs) explicitly solve for synoptic scale eddies and their monthly averages can reasonably well reproduce the climate statistics of the real atmosphere (Peixoto and Oort, 1984;Roeckner et al., 1996;Laine et al., 2009).To explicitly resolve the eddies, AGCMs require a high resolution in space and time and hence need ample computational time.The fact that AGCMs are able to reproduce key atmospheric features like the mid-latitude storm tracks by solving the principle equations makes them essential tools in studying the Earth's climate.Still, substantial structural differences may exist between observed and AGCM-simulated patterns of the synoptic scale activity (e.g.their kinetic energy) for extended atmospheric domains notably in the North Pacific and in the upper troposphere (Robinson and Black, 2005).Also, due to their computational expense, both the number of simulations and the simulated time period are limited.
In contrast to the AGCMs, statistical-dynamical atmospheric models (SDAMs) are based on time-averaged equations in which the effects of transient eddies are parameterized in terms of the large scale field (Saltzman, 1978;Peixoto and Oort, 1984).These models directly resolve the ensemble-averaged synoptic eddy characteristics (second and higher order moments) rather than calculating the contribution of each individual storm.The essential difference with AGCMs is thus the point of truncation in the frequency spectrum of atmospheric motion (Saltzman, 1978).This different approach allows much coarser spatial and temporal discretizations, making SDAMs computationally efficient and allowing climate simulations up to multi-millennia timescales (Kubatzki et al., 2000;Bauer et al., 2004).For this reason, several of the Earth System Models of Intermediate Complexity (EMICs) (Claussen et al., 2002;Petoukhov et al., 2005), including Climber-2 (Petoukhov et al., 2000;Ganopolski et al., 2001), IAP RAS (Petoukhov et al., 1998;Mokhov et al., 2002;Eliseev et al., 2008), MIT (Prinn et al., 1999;Kamenkovich et al., 2002) and UVic (Weaver et al., 2001), use statistical-dynamical atmosphere modules.
A major challenge in developing SDAMs is to derive useful parameterizations for the synoptic scale moments appearing in the equations (Saltzman, 1978).Traditionally, the most common way to achieve this has been by means of a diffusive approximation for the second order moments (Green, 1970;Branscome, 1983;Pavan and Held, 1996;Handorf et al., 1999;Barry et al., 2002).Early diffusive models described synoptic scale moments only in zonally averaged form, as a function of the gradients of zonal wind (Starr, 1968), temperature (Williams and Davis, 1965), or the combination of them (Green, 1970;Stone, 1972a,b;Branscome, 1983) and thus lacked geographical variation.However, based on observational data and primitive-equation model results (Blackmon et al., 1977;Gall, 1976a,b;Schneider, 1981), the importance of longitudinal variations in the climatological zonal flow was soon acknowledged (Niehaus, 1980;Fredericksen, 1983;Robertson and Metz, 1990).
An alternative to the diffusive models is provided by baroclinic adjustment models, which postulate an equilibrium between baroclinic wave forcing and the large scale atmospheric temperature field (Stone, 1978;Cehelsky and Tung, 1991).Farrell and Ioannou (1996a,b) were the first to develop linear stochastic models, which parameterize the energy and momentum transfer, associated with nonlinearity of synoptic eddies, as a combination of stochastic excitation and enhanced damping.Using such models, Farrell and Ioannou (1996a,b) showed that synoptic eddies can develop on a baroclinically stable mean flow.Their approach, which used idealized zonally symmetric flows, was extended by Whitaker and Sardeshmukh (1998) to zonally varying background flows with prescribed Gaussian white noise excitation and linear damping.Applying the same method, Zhang and Held (1999) performed a case study of the so-called "midwinter suppression" of the Pacific storm track (Nakamura, 1992), and DelSole ( 2001) made an attempt to specify the statistics of the stochastic forcing as a function of the background flow and low-level eddy heat flux.
The goal of this study is to derive a set of statisticaldynamical equations which properly capture the synoptic field for climates not too different from present day, i.e. typically from the last glacial maximum (LGM) to a doubled CO 2 -concentration climate (2×CO 2 ).This is part of our effort to develop an Earth System Model of Intermediate Complexity (EMIC) which can be used to do long-term simulations within this climatic range.Starting from early work by Lau (1979) and Lau et al. (1978), we derive a set of equations for the kinetic energy and momentum flux due to transient synoptic eddies as well as parameterizations for the vertical synoptic and ageostrophic terms.This set of equations is implemented in a newly developed SDAM named Aeolus.We compare synoptic fields generated by Aeolus with those of reanalysis and GCMs for present-day, LGM and 2×CO 2 climates, and show that they are in good agreement.

Governing equations
We start from the principle equation governing atmospheric velocity (V = {u,v,w}) in the hydrostatic approximation (Lorenz, 1967), as given by Eq. (A3) in Appendix A. A key assumption of our approach is that we can split V into a large-scale long-term ( V ) and a synoptic scale (V ) component such that V = V +V = { u , v , w }+{u ,v ,w }.As our focus is on parameterizing the synoptic scale field, for the remainder of this article we assume the large scale field to be known.Though in a 1-D atmospheric power spectrum (i.e. as function of period or wavelength only) maxima from long-term (∼30 days period) and synoptic-scale (2-6 days period) fields are not well separated (Van der Hoven, 1957;Oort and Taylor, 1969;Vinnichenko, 1970;Mitchell, 1976), in 2-D spectra (versus both period and wavelength) the gap is more pronounced (Fraedrich and Böttger, 1978).Even more, in a 3-D power spectrum, with phase velocity as a third dimension, the long-term field and the field of transient synoptic scale eddies form distinct and well-separated maxima (Petoukhov, 1991;Petoukhov et al., 1998Petoukhov et al., , 2003)), indicating that they are governed by different physical processes.Therefore, splitting V into a large-scale and synoptic component and parameterizing the ensemble characteristics of the synoptic scale eddies in terms of the long-term field is principally allowed.By applying the appropriate averaging procedures, we can derive statistical-dynamical equations for the averaged kinetic energy of transient eddies E k defined as and the momentum flux by transient eddies u v (see also Lau, 1978): Here, K fh and K fz are internal atmospheric small/meso-scale friction coefficients in the horizontal and vertical direction respectively, K fs is the surface friction coefficient, f the Coriolis parameter and subscript "ag" denotes ageostrophic terms.In Eqs.
(1) and (2) (equivalent to Eqs.A43 and A44 in Appendix A), third order moments as well as the metric terms have been neglected.The fourth till sixth terms on the right hand side of Eqs. ( 1) and ( 2) are friction terms whereas the first describes advection by large-scale, longterm wind.In Eq. ( 1), the second and third term on the righthand side capture fluxes coming from horizontal and vertical shears in the large scale flow.These thus capture the exchange of kinetic energy and momentum flux between eddies and the "basic" flow, due to its horizontal and vertical shear.Lau et al. (1978) neglect the vertical components of these flux terms based on the assumption that u w and v w are of the same order.This assumption, however was never supported by empirical observations.In fact, the continuity equation in combination with the quasi-geostrophic character of the transient synoptic eddies implies that the correlation between v and w is generally much higher than that between u and w .This results in particular in much higher correlations between v w and ∂ v ∂z as compared to those between u w and ∂ u ∂z .Figure 1 plots fluxes arising from zonal, meridional and vertical shear in the large scale wind field, which enter the right hand side of Eq. (1), calculated from band-pass filtered ERA-40 data (Petoukhov et al., 2008).This figure shows that zonal and meridional flux terms tend to have an opposite sign reducing their combined magnitude, most notably in the North Atlantic region.Further, it shows that the vertical flux terms are in fact of the same order as the horizontal ones.We thus emphasize, based upon theoretical and observational considerations, that the vertical flux terms should not be neglected.Motivated by Williams and Davis (1965), we describe their sum as (see Eq. A7): where K syn is the vertical macro-turbulent diffusion coefficient (Williams and Davis, 1965).We assume equipartitioning of this vertical (baroclinic) flux term between the zonal and meridional kinetic energy components, allowing us to split prognostic Eq. (1) into two separate ones for u 2 and Observational data indicates that cross-terms between these two prognostic equations are generally small and only become important when the eddy kinetic energy is small, i.e. in near-equatorial regions.Here we aim at parameterizations which hold in the mid-latitudes, where the cross-terms can be neglected.In doing so we follow earlier work by Eady (1949) and Kurihara (1970).The same considerations imply that the vertical fluxes in Eq. ( 2) should be approximately equal ( u w ∂ v ∂z ≈ v w ∂ u ∂z ), and hence the flux maintenance equation can be reduced to: We can express v w in terms of v 2 using the continuity equation and assuming a quasi-geostrophic character of synoptic-scale motion: with z st as the "steering" level (Chan and Gray, 1982) for synoptic scale eddies, H 0 the air density scale height and γ (1 + ζ,z/H 0 ) the lower incomplete gamma-function (Abramowitz and Stegun, 1964).
To close the set of Eqs. ( 4)-( 7), we derive suitable parameterizations for the vertical macro-turbulent diffusion coefficient (K syn ) and ageostrophic terms ( u v ag − v u ag and v v ag − u u ag ), resulting in 4 Eqs.( 4)-( 7) and 4 unknowns ( u 2 , v 2 , u v and v w ).We assume internal friction coefficients (K fh , K fz ) to be constant.Appendix A describes in detail the derivation of parameterizations for ageostrophic terms (Eqs.A45-A34) and K syn (Eq.A11).

Present-day climate
Equations ( 4)-( 7), and all parameterizations needed (Appendix A), were implemented in the SDAM Aeolus and numerically solved on a regular 3.75 • × 3.75 • -grid consisting of 5 pressure levels (850 mb, 700 mb, 500 mb, 400 mb and 300 mb), for each calendar month.Monthly averages of large-scale long-term wind field and surface temperature, which enters the equation for the Rossby length scale (Eq.A12), were taken from the ERA-40 reanalysis dataset for 1976-2002 (Uppala et al., 2006).We compared the calculated numerical solutions for the synoptic scale parameters against empirical data, which were obtained by 2-6 day bandpass filtering of 6-hourly ERA-40 windfield data (Petoukhov et al., 2008).Two different numerical experiments were performed.First, an "uncoupled" test experiment, in which the synoptic scale terms appearing in the right-hand side of Eqs. ( 4)-( 7) were taken from empirical data.Second, a "coupled" experiment, in which the righthand-side synoptic moments were calculated using Eqs.( 4)-(7) as well.Thus, in the uncoupled experiment, each of the synoptic moments ( u 2 , v 2 , u v and v w ) is calculated for prescribed large-scale wind field ( u and v ) and prescribed other synoptic moments.In the coupled experiment only large-scale wind is prescribed and the synoptic moments are calculated interactively.Time-integration of Eqs. ( 4)-( 6) is achieved by means of finite volume and finite difference methods, whereby initially the prognostic variable is zero everywhere.Table 2 provides values for all constants appearing in the equations as used in both experiments.We present "coupled" solutions for synoptic scale kinetic energy and "coupled" and "uncoupled" solution for momentum transport u v and v w .Note that the terms involving u v in the right hand side of Eqs. ( 4) and ( 5) are an order of magnitude smaller than other terms and hence "coupled" and "uncoupled" solutions for E k are essentially the same.
Figure 2 compares empirical data versus Aeolus results for the synoptic scale kinetic energy (E k ) at three different pressure levels for Northern-Hemisphere winters.In general, the areas of high synoptic activity are accurately captured, both geographically as well as their vertical variation.In the Northern Hemisphere, the location of the key storm track regions in the Atlantic and Pacific oceans at about 40 • N are well resolved.The amplitude of the Pacific one is somewhat overestimated, whereas the amplitude is Kinetic energy of synoptic transients Eq. ( 1) Zonal kinetic energy of synoptic transients Eq. ( 4) v 2 m 2 s −2 λ,φ,z Meridional kinetic energy of synoptic transients Eq. ( 5) u v m 2 s −2 λ,φ,z Momentum flux by synoptic transients Eq. ( 6) u w m 2 s −2 λ,φ,z correlation between u and w Eq. ( 3) v w m 2 s −2 λ,φ,z correlation between v and w Eq. ( 7)  (d-f) show the results from Aeolus (Eq. 1) solutions at same pressure levels.
underestimated above southern Europe.In the model, synoptic activity is more confined to the mid-latitudinal storm track regions, whereas in observations it is more diffused with a broader latitudinal extend.In the Southern Hemisphere, the elongated storm track region between 40  4b).However, synoptic activity in the model is much more confined to a region stretching from southern Africa to southern Australia, whereas in the observations activity is circum-global.
In Fig. 4, the zonally and vertically averaged meridional distribution of E k are shown for Northern-Hemisphere   v w (c, d) and u v (e, f), all in m 2 s −2 , for December-January-February (a, c, e) and June-July-August (b, d, f) averages, plotting ERA-40 (solid), "coupled" (dotted) and "uncoupled" (dashed) numerical solutions.(d-f) show results of "uncoupled" and (g-i) "coupled" experiments using Aeolus (Eq.2) at same pressure levels.
winters (Fig. 4a) and summers (Fig. 4b).The location and amplitude of the mid-latitude storm regions are well resolved, as well as the low-energetic near-equatorial regions.Only in the near-polar region (>70 • ), the model tends to underpredict the synoptic scale kinetic energy.Figure 4e and f plots zonal-mean u v in boreal winter-and summer respectively, showing that both the "uncoupled" and "coupled" numerical solutions of Eq. ( 6) are in good agreement with observations.Also geographically, the fit between observed and model-generated u v is reasonable, for both "uncoupled" and "coupled" experiments (Figs. 5 and 6).The main high-amplitude regions, associated with storm track regions, are well captured geographically and also their vertical increase is in agreement with empirical data.Again the model generally produces u v patterns which have a less diffusive character than those observed.The amplitude and meridional extend of the Southern Hemisphere storm track region are correctly resolved in both the "uncoupled" and the "coupled" model.In the Northern Hemisphere, results are generally similarly good, though the "coupled" results predict too low u v over central Europe in wintertime (Fig. 5f).In the "uncoupled" experiment, over some regions at high altitude, a sign inversion is observed: i.e. south-east of Australia in JJA (Fig. 6b) and Central Asia in DJF (Fig. 5b).This is not observed in the "coupled" experiment.Shown in Fig. 4e and  f, the zonally and vertically averaged solutions of u v are in good agreement with those from observations.Finally, Figs. 7 and 8 compare "coupled" and "uncoupled" model results of v w with bandpass filtered results, respectively for Northern Hemisphere winter-and summertime.For both periods, the model accurately captures location, magnitude and vertical variation of v w ."Coupled" results show a narrower band of activity in the Northern Hemisphere compared to observations, which is a direct consequence of the same model behavior of E k (Fig. 2).In the "uncoupled" experiments this effect is less pronounced.Nevertheless, zonal averages (Fig. 4c and d), geographical spreading and seasonality of v w are in good agreement with ERA-40 data for both "coupled" and "uncoupled" experiments.

Other climates
In the previous section, we showed that the presented statistical-dynamical equations for kinetic energy and momentum flux perform well under present-day climate conditions.In addition, we tested their behavior in warmer and colder climates and compared results against state-of-the-art atmospheric general circulation models (GCMs).(d-f) show results of "uncoupled" and (g-i) "coupled" experiments using Aeolus (Eq.2) at same pressure levels.
As a case-study for warmer climate, we used the IPCC (Intergovernmental Panel on Climate Change) 2×CO 2 equilibrium experiment (IPCC-4AR, 2007), which represents a nearly 4 • C warmer world (yearly-averaged global mean near-surface air temperature) in the GCM reference simulations used here.On the one hand, we bandpass filtered daily windfield data from the GFDL-CM2.0(Delworth et al., 2006) and ECHAM5/MPI-OM (Roeckner et al., 2003), 2×CO 2 equilibrium experiments as well as their 1×CO 2 control runs (all available at http://www-pcmdi.llnl.gov) to extract their synoptic scale kinetic energy.On the other hand, we used the large-scale, long-term wind field (monthly averages) from the GCM simulations to force Aeolus, just as we used ERA-40 large-scale windfield for the present-day climate experiment.Figure 9a-f shows the results of this experiment, plotting the synoptic scale kinetic energy anomaly, i.e. the energy in the 2×CO 2 climate minus the 1×CO 2 climate.
Generally storm activity shifts polewards with an abatement around 40 • and an intensification around 60 • in both hemispheres.(Fig. 9c and f).This poleward shift has been discussed in detail by Yin (2005) and is more pronounced in the GFDL simulation than in the ECHAM5 one.Geographically, the two GCMs can differ substantially from each other over extended regions, most notably in the Northern Hemisphere.The zonally averaged solution of Aeolus fits reasonably well with the GCM solutions (Fig. 9c and f).The locations of abatement (near 40 • ) and intensification (near 60 • ) are predicted correctly.The sensitivity, i.e. the magnitude of abatement or intensification in the 2×CO 2 climate, of Aeolus seems slightly higher than GFDLs and slightly lower than ECHAM5's sensitivity.These small differences might be caused by differences in the description of the static stability parameter in the three atmosphere models (see Eq. A13).Geographically we capture the most prominent features, though some differences exist.The Southern Hemisphere, where zonal variation is limited, is rather well resolved.In the Northern Hemisphere, which has much more zonal variation and the GCMs differ themselves substantially, the most important features are captured.For example, the ECHAM5 simulation shows strong intensification in three mayor regions (above the UK, Canada and the northern Pacific) which are all present in our results.A notable difference occurs in the eastern Mediterranean region where the anomaly from Aeolus forced with GFDL larg-scale wind has an opposite sign compared to the GFDL anomaly.Note here however, that also ECHAM5 predicts an opposite anomaly in this region as compared to GFDL.We performed an analogous  (d-f) show results of "uncoupled" and (g-i) "coupled" experiments using Aeolus (Eq.7) at same pressure levels.(d-f) show results of "uncoupled" and (g-i) "coupled" experiments using Aeolus (Eq.7) at same pressure levels.
experiment for a last glacial maximum (LGM) climate, for which we used HadCM3M2 (Pope et al., 2000), LGM and Pre-Industrial simulations from the Paleoclimate Model Intercomparison Project (PMIP2) database (Braconnot et al., 2007).The results of this experiment (Fig. 9g-i) show that storm activity is reduced substantially in an LGM climate as compared to pre-industrial, which is in agreement with previous studies (Laine et al., 2009).Aeolus rather accurately reproduces the zonally averaged HadCM3M2 results though the small intensification (∼ 2 m 2 s −2 ) at ∼40 • S and ∼40 • N are not resolved.We associate the observed differences in the Northern Hemisphere to dynamics associated with the thick icesheets not captured in out model.Geographically Aeolus manages to capture the most important features, though some differences are observed too, resulting in a fair comparison.
Finally, we compared the seasonal cycle produced by Aeolus for the three climates, i.e. 1×CO 2 , 2×CO 2 and LGM, with those from GCMs and ERA-40.Here, we define the synoptic scale kinetic energy anomaly as the energy in winter time (December-January-February averaged) minus the energy in summer time (June-July-August averaged).Figure 9j, k and l plots this kinetic energy anomaly in zonally averaged format for 1×CO 2 , 2×CO 2 and LGM, respectively.These show that the seasonality produced by Aeolus is comparable to that of the GCMs and ERA-40.

Sensitivity experiments
We performed sensitivity experiments to determine the effect of changes in the value of parameters on the synoptic scale kinetic energy.The last 8 parameters listed in Table 2 affect kinetic energy.The characteristic drag coefficient CD has the equivalent, though opposing, effect as the characteristic eddy diffusion coefficient in the PBL, KPBL , as can be seen from the equation of the cross-isobar angle (A31).Also, friction coefficients K H and K 2,z together capture the magnitude of internal small-scale friction in the atmosphere.This thus reduces the total number of free parameters to 6.The sensitivity of the eddy kinetic energy E k to each of these parameters is shown in Fig. 10.This figure plots the change in E k when the value of a parameter is changed from minus 10 % to plus 10 % of the value given in Table 2, thus covering a total parameter range of 20 %.
E k is most sensitive to changes in static stability parameters m and U 0 .If m, which captures the larger static stability high-up in the troposphere, increases, so will the static  and (d, e) and (g, h), respectively with GCM (solid) and Aeolus (dashed).Plots (j-l): seasonal E k -anomaly (DJF-JJA) for (j) 1×CO 2 , (k) 2×CO 2 and (l) LGM climate.In (j) and (k): GFDL (solid), ECHAM5 (dashed), Aeolus (dotted), in addition in (j) ERA-40 (dashed-dotted).In (l) HadCM3 (solid) and Aeolus (dashed).stability, and hence E k will decrease as seen in Fig. 10.A 20 % increase in m reduces E k by about 30 m 2 s −2 , or by roughly 50 % of its value.Vice versa, when U 0 increases, static stability will decrease and thus E k will increase.Gates (1961) argues that static stability increases by roughly a factor 4 between sea-level and the tropopause.In our parameterizations this is captured by m = 0.66.E k is much less sensitive to the choice of other parameters.Surface friction, via τ E,min , has a stronger effect than internal meso-scale friction, via K H and K z,2 .An increase in the last tends to broaden the storm tracks with a decrease in E k in their centers and an increase at their periphery.The surface friction parameter τ E,min has a near-linear effect on E k with a change in τ E,min causing a similar magnitude change in E k .Finally, the vertical eddy diffusion coefficient in the PBL, K PBL , and the meso-scale eddy diffusion coefficients have a limited effect on E k ., (d) m, (e) U 0 and (f) K meso .Plotted is the change in E k when the value of the parameter is changed from minus 10 % to plus 10 % the value given in Table 2, thus covering a range of 20 %.

Discussion
We presented and tested a set of statistical-dynamical equations which can accurately reproduce the geographical and vertical distribution of storm tracks, including their seasonality.We close the set of equations by finding physicallyrealistic parameterizations for the vertical macro-turbulent diffusion coefficient (K syn ) and the ageostrophic terms.These parameterizations are used for solving all synoptic quantities ( u 2 , v 2 , u v and v w ) and do not require any individual tuning to obtain proper results.Also, those parameterizations, as well as all constants (Table 2), are valid globally without any local adjustments.This way we obtain proper results for the four synoptic scale quantities for present-day, 2×CO 2 and LGM climates.
Nevertheless, direct comparisons of the parameterized intermediate terms (i.e.K syn , u w , v v ag − u u ag ) versus observations or independent estimates would further validate this model.This is complicated as these quantities cannot straight-forwardly be calculated from meteorological datasets like ERA-40.To our knowledge, only few studies have published reliable estimates of them so far, notably Tucker (1960) for K syn , Lau et al. (1978) for v v ag − u u ag and Stone and Yao (1987) for u w .Tucker (1960) and Stone and Yao (1987), unfortunately, used unfiltered data for their estimates hence including the effect of standing eddies as well.A comparison therefore only makes sense for the mid-latitudes as here kinetic energy of transient eddies dominates over that of standing eddies.Here Tucker (1960) estimates the zonal average of K syn to be roughly 70 m 2 s −1 at 400 mbar (see Table 6 of Tucker (1960)) which is in fair agreement with our estimate of ∼ 80 m 2 s −1 .Similarly, the values and vertical increase of zonally averaged u w given by Stone and Yao (1987) (their Fig. 4) are in good agreement with ours outside the near-tropical region.Lau et al. (1978) published estimates of (bandpass filtered) f ( v v ag − u u ag ) in the Northern Hemisphere, including its geographic variation (Fig. 4 of Lau et al., 1978).This figure is in quantitative agreement with our estimates: Above the two ocean basins negative anomalies form at high latitude (east of Newfoundland and east of the Kamchatka Peninsula) and positive anomalies at lower latitudes (east of Florida and east of China).Also the magnitude of these anomalies is in good agreement, with f ( v v ag − u u ag ) ranging from −4 × 10 −3 m 2 s −2 to 4 × 10 −3 m 2 s −2 .The presented intermediate parameterizations are thus in qualitative (for K syn and u w ) and quantitative (for ageostrophicity) agreement with previously published data.
Early studies (e.g.Stone and Yao, 1987) based parameterizations upon the idea that the linear stage of baroclinic instability of the quasi-zonal wind is the dominant one shaping the structure of the energy-loading synoptic eddy ensemble characteristics.Here we argue that by accounting for the nonlinear baroclinic and barotropic forcings during the mature stage of synoptic eddies, one can satisfactorily reproduce the basic features of the geographical and vertical structure of the synoptic eddy-ensemble kinetic energy and momentum flux.This is in line with observational data and the modeling results of the individual baroclinic wave life cycle (Hoskins, 1983).Hence, we argue that in the equations for u 2 , v 2 and u v , the vertical terms, u w and v w , cannot be neglected (Fig. 1).Linear stochastic models (Farrell and Ioannou, 1996a,b) implicitly assume a strong upward transport of synoptic eddy activity in order for their mechanism of equilibration of the mean flow to be efficient.This agrees with our findings.The key difference is that our parameterizations explicitly describe both the vertical transport of synoptic activity and the nonlinearity of the processes driving synoptic activity, without prescribing the magnitude of these terms.In addition, we take into account the ageostrophic terms appearing in the equations, whose importance was first claimed in the classical papers by Lau et al. (1978) and Lau (1979).Orlanski and Chang (1993) and Chang and Orlanski (1993) emphasized the crucial role the ageostrophic component plays in the evolution of the individual baroclinic waves and eddies.In our study, we found that these terms are also important to correctly describe the ensemble-averaged u 2 , v 2 and u v .Finally, we do not use any a priory assumption on the stability of the basic flow, making our formulas rather general and therefore in principle applicable to a broad range of basic flow profiles.
With calculation times a fraction of those of GCMs, Aeolus can, when fully developed, be used to study Earths climate over multi-millennia and longer timescales.The computational expense of GCMs has hampered such studies so far.On the other hand, existing EMICs, though computationally efficient, lack accurate 3-D parameterizations for synoptic scale kinetic energy and momentum and thus cannot produce the climatic effect of storm tracks.This work thus prepares ground for the development of computationally efficient EMICs which can accurately resolve the climatic effects of changes in atmospheric dynamics, including its synoptic component.We acknowledge that we tested our parameterizations only for a moderate climatic range here, which does not guarantee proper solutions outside this range.However, it is exactly this climatic range for which we want to use Aeolus, as part of a new EMIC, to address scientific questions with respect to future climate stability and the origin of the glacial cycles (Ganopolski and Roche, 2009).Furthermore, useful applications will include ensemble simulations to reduce uncertainty in climate sensitivity (Schneider von Deimling et al., 2006) and other governing parameters of climate models (Eliseev, 2008).Based on this approach similar parameterizations for the synoptic scale heat-and moisture-flux can be developed.

Conclusions
We derived a closed set of 3-D statistical-dynamical equations describing the kinetic energy and momentum flux of synoptic scale eddies in the Earth's atmosphere by find-ing parameterizations for the vertical macro-turbulent diffusion coefficient and ageostrophic terms.We show that these statistical-dynamical equations fairly accurately reproduce the geographical distribution of storm tracks, including their variation with height and season.This works both for present-day climate as well as moderately colder (LGM) and warmer (2×CO 2 ) ones.Coupled to ocean and landsurface models, such atmospheric modules are likely to be efficient tools for studying the Earths climate over millennia and longer timescales.Such long-term, ensemble simulations are required to address key questions with respect to future climate stability, on the one hand, and the origin of the glacial cycles, on the other.

Appendix A Derivation of statistical-dynamical pde's A1 The statistical-dynamical approach
A key assumption of our statistical-dynamical (SD) approach to climate modeling is that any original large-scale atmospheric variable Y i can be written as a sum of a large-scale long-term (quasi-stationary) component ( Y i ) and a synoptic component (Y i ) (Petoukhov et al., 1998(Petoukhov et al., , 2003)): Here Y i has a characteristic time scale (or period) τ ≥ τ 0 Y ≈ 15 days, spatial scale (or wavelengths) L ≥ L0 Y ≈1500 km and phase speed C < C0 Y ≈3 m s −1 (Fraedrich and Böttger, 1978;Pedlosky, 1979).The dominant contribution to the perturbation term Y i is assumed to come from transient synoptic eddies (cyclones and anticyclones) with characteristic scales: τ ≤ τ 0 Y ≈ 6 days, L ≤ L0 Y ≈800 km, and C ≥ C0 Y ≈10 m s −1 (Eady, 1949;Fraedrich and Böttger, 1978;Pedlosky, 1979).Crucially important here is the clear separation between large-scale longterm and synoptic fields in 3-dimensional ( τ -L-C) atmospheric power spectra justifying our approach (see e.g. the power spectra plots given by Fraedrich and Böttger (1978) with phase velocity as the third axis).Notably, in twodimensional power spectra, without a phase velocity axis, the mentioned separation is not clearly marked (Fraedrich and Böttger, 1978).Even more, in one-dimensional power spectra, versus period or frequency only, the fast synopticscale and large-scale long-term components are practically non-separated ( Van der Hoven, 1957;Oort and Taylor, 1969;Vinnichenko, 1970;Mitchell, 1976).Only in 3-D phase space becomes the separation apparent.The averaging operator ... therefore denotes a temporal, spatial and phasespeed averaging over the 3-dimensional window bounded by equalities are direct consequences of the separation between Y i and Y i in the τ -L-C-phase space (Petoukhov, 1991;Petoukhov et al., 1998Petoukhov et al., , 2003)): for any i , where ∇ and are respectively the 3-D gradient and Laplacian operator, and F , ℵ and F i represent any type of algebraic functions or differential operators.Further, in the stratosphere, all synoptic moments Y i Y j are assumed to obey the Charney-Drazin formula (Charney and Drazin, 1961), describing exponential decay upwards in the stratosphere (Petoukhov et al., 2003).

A2 Synoptic scale kinetic energy
To derive an equation for the synoptic scale kinetic energy, E k , we start from the original equation for atmospheric wind in the hydrostatic approximation (Lorenz, 1967) Here, V H = {u,v} is the horizontal velocity vector, f = 2 sinϕ is the Coriolis parameter, with the Earth's angular velocity and ϕ latitude, k is the upward pointing unit-vector, ρ is air density, p is pressure, ∇ H is the horizontal gradient operator and F h is the friction force per unit mass acting on V H .In the zero-order inelastic approximation, ρ can be replaced by ρ 0 = ρ 0 (0)exp(−z/H 0 ) with constant reference density ρ 0 (0) and air density scale height H 0 .First, we insert V H = V H +V H into Eq.(A3).Next, we take the scalar product with V H , apply the averaging operator ... to both sides of Eq. (A3) and neglect all higher order terms given by Eq. (A2), giving: A2.1 Deducing the left hand side of Eq. (A4) Expanding the total derivatives in the left hand side of Eq. (A4) into their partial derivatives and using A2, we can rewrite it as follows: Here, E k = 1 2 ( u 2 + v 2 ) is the synoptic scale kinetic energy, with u and v the zonal and meridional components of V H .The last term in Eq. (A5) appears due to differentiating in spherical coordinates, referred to as the "metric" term.The second and third terms on the right hand side describe respectively the advection of synoptic scale kinetic energy by large-scale long-term motions and by ensembles of the synoptic eddies themselves (self-advection).Finally, the fourth and fifth term can be expanded as: The first three terms on the right-hand side of Eq. (A6) represent exchange of kinetic energy between eddies and the "basic" flow due to its horizontal shear, in particular feeding the eddy kinetic energy at the expense of that of the "basic" flow in case of barotropic instability of the "basic" flow (Pedlosky, 1979).The last two terms on the right-hand side of Eq. (A6) capture the exchange of energy between eddies and the "basic" flow owing to its vertical shear and therefore, due to a quasi-geostrophic character of the "basic" flow in the middle and high latitudes, owing to horizontal gradients of the large-scale long-term temperature field.The terms u w and v w together describe the direct transfer of kinetic energy from the "basic" flow to the ensembles of synoptic eddies (Stone, 1972a,b).This process is referred to as "synoptic-eddy friction" by Williams and Davis (1965), who express u w and v w in terms of ∂ u ∂z and ∂ v ∂z , respectively.Thus, the last two terms in the righthand side of Eq. (A6) capture the combined effect of "synoptic-eddy friction" and "basic" flow baroclinic instability, if the last occurs.In the general case, regardless of whether the "basic" flow is baroclinically stable or instable, these terms describe the exchange of energy between the "basic" flow and synoptic eddies, due to vertical shears of the former.Based on ideas of Williams and Davis (1965) and Petoukhov (1991), we derive a new expression for the vertical terms by parameterizing their sum: with K syn the vertical macro-turbulent diffusion coefficient (Petoukhov, 1991;Petoukhov et al., 2003).

A2.2 K syn parameterization
In line with "vertical mixing" of angular momentum due to synoptic scale eddies, introduced by Here K meso represents K syn from the background vertical mixing due to mesoscale eddies only and is taken constant here.With the use of the continuity equation, div(ρ 0 V ) = 0, we can estimate W as: where Ũ is the absolute value of the characteristic horizontal velocity in the synoptic transients, which is roughly proportional to u : This assumption is supported by observational data of the distribution of the characteristic magnitudes for winds for the "basic" flow and synoptic eddy ensembles (Oort and Rasmusson, 1971;Oort, 1983).Further, L ϕ Ro represents the characteristic meridional spatial scale of synoptic eddy displacement which we set to the internal Rossby radius of deformation in the meridional direction.Accounting for the latitudinal dependence of longitudinal-latitudinal asymmetry of synoptic eddies, we define: Hoskins et al. (1983) with ϕ 0 corresponding to 45 • latitude.F K (z) captures the vertical change in W and will be derived later.First we parameterize the characteristic synoptic-eddy vertical scale H , which we express in terms of the density scale height H 0 and the meridional Rossby number, or: Inserting these parameterizations for W and H into Eq.(A8), we obtain: The Rossby radius of deformation in the zonal direction ( L λ Ro ) can be written in terms of surface temperature ( T s ), the specific gas constant (R) and the Coriolis parameter (f ): with the static stability parameter α : with the non-dimensionalizing constant Ũ0 = 5 m s −1 and g the acceleration due to gravity. is the mid-troposheric lapse rate and a and wa are, respectively, the dry adiabatic and mid-troposphere moist adiabatic lapse rate.Further, V cl is the relative volume occupied by the middlelayer stratus and penetrative-convection cumulus clouds ("moist-convection towers") within the troposphere.The influence of the atmospheric moisture on the stability of the atmosphere is thus accounted for.Furthermore, the factor ( Ũ / Ũ0 ) m , with m ≈ 2/3, allows for the effect of the increase in the static stability parameter in the upper troposphere as compared to the middle troposphere (Gates, 1961).The difference between our parametrization of α with those of e.g.Dymnikov (1978), Emanuel et al. (1987) and Lapeyre and Held (2004), is that ours is relevant for climate time scales while theirs are more relevant for fast weather processes.
Using Eq. (A12), we can rewrite the denominator in Eq. (A11) as follows: |sinϕ| sinϕ 0 capturing both the effect of the latitudinal change in the Coriolis parameter and the synoptic eddies assymetry.We apply a coarse-resolution representation for F 0 (ϕ) based on the observational facts that synoptic eddies (1) possess a finite size in their mature stage, (2) seldom cross the equatorial plane and (3) seldom run over-the-pole (Palmen and Newton, 1969;Pedlosky, 1979).We treat F 0 (ϕ) separately in three latitudinal ranges: the near-equatorial belt (F 0 → 0), the mid-latitudes and the nearpolar regions (F 0 → 2).Based on observational facts (1) and (2), we set the width of the near-equatorial belt to the size of the eddies in this region, i.e.: ϕ e = L ϕ Ro (ϕ e )/a, where the horizontal bar indicates a zonal average.This is equivalent to the equatorial β-plane approximation for the synoptic eddy internal Rossby deformation radius.The value of F 0 within the near-equatorial belt is set to F 0 (ϕ e ).Analogously, motivated by observational fact (3), we prescribe the width of the near-polar region to the size of the eddies there, i.e. ϕ p = L ϕ Ro (ϕ p )/a and F 0 gets the value F 0 (ϕ p ) in the region from ϕ p to the pole.In the mid-latitudes, i.e. from ϕ e to ϕ p , we use a method which can be seen as a composite of a f 0plane and a β-plane approximation.We define an additional latitude ϕ max to be the latitude of maximum synoptic activity, with F 0 = F 0 (ϕ max ).F 0 (ϕ) in the mid-latitudinal belt from ϕ e to ϕ max is now defined as: where the subscript ag refers to the ageostrophic component of u and v , which play an important role as source/sink terms (Lau et al., 1978;Opsteegh and van den Dool, 1979).We assume that the vertical profile of ageostrophic wind can be described by trigonometric functions such that the horizontal mass flux over the total atmospheric column vanishes.We first give the full derivation of these equations, after which we show that our approach is conceptually in line with the atmospheric energy conversion cycle presented by Lorenz (1967).Following Petoukhov et al. (2000Petoukhov et al. ( , 2003)), the components u ag and v ag of an instantaneous synoptic scale ageostrophic velocity in the free troposphere are expressed as: u ag = e β 0 z 1 [A 1 (λ,ϕ)cos(β 0 z 1 ) − A 2 (λ,ϕ)sin(β 0 z 1 )]+ e −β 0 z 1 [A 3 (λ,ϕ)cos(β 0 z 1 ) + A 4 (λ,ϕ)sin(β 0 z 1 )] v ag = e β 0 z 1 [A 2 (λ,ϕ)cos(β 0 z 1 ) + A 1 (λ,ϕ)(sinβ 0 z 1 )]+ e −β 0 z 1 [A 4 (λ,ϕ)cos(β 0 z 1 ) − A 3 (λ,ϕ)(sinβ 0 z 1 )] (A28) Here 2K syn is analogous to the inverse Ekman height in the theory of the planetary boundary layer (PBL), but used here for the free troposphere.In Eq. (A28), z 1 = z − z 1 where z 1 stretches from the surface till the lower most "turning point" of the ageostrophic velocity at the PBL-free troposphere interface.This thus captures the change in sign of the monsoonal branch of the atmospheric circulation.At this "turning point", ageostrophic circulation vanishes, i.e. u ag (z 1 = 0) ≈ 0 and v ag (z 1 = 0) ≈ 0 and thus A 1 + A 3 ≈ 0 and A 2 + A 4 ≈ 0. We also impose that, over the total atmospheric column, the total horizontal mass flux due to ageostrophic circulation vanishes, i.e.: u ag e −z 1 /H 0 dz 1 = 0 v ag e −z 1 /H 0 dz 1 = 0 (A29) where u ag z 1 and v ag z 1 are the PBL-averaged zonal and meridional ageostrophic velocities and z up = H − z 1 with H the height of the upper boundary for the atmopsheric dynamics computational domain.With these imposed conditions, we can determine A 1 (λ,ϕ), A 2 (λ,ϕ), A 3 (λ,ϕ) and A 4 (λ,ϕ).The PBL-averaged ageostrophic components obey the equations (Petoukhov et al., 2003): Here, u z 1 ( v z 1 ) is an instantaneous geostrophic, zonal (meridional) component of the synoptic scale velocity averaged over z 1 , sgn(f ) equals +1 (−1) for the NH (SH), and C α is a dimensionless constant of the order unity.α 0 is the near-surface cross-isobar angle, which to a first order approximation based on findings of Hansen et al. (1983) obeys: Parameters CD and KPBL are characteristic values of the drag coefficient and the kinematic vertical small/mesoscale eddy diffusion coefficient within the PBL respectively (Hansen et al., 1983).Subscript 1 refers to variables in the PBL.Inserting Eqs.(A28) and (A30) into Eq.(A29) we get: where F (z) equals one within the PBL.Above it, it is given by: The full derivation of F is tedious but directly follows from Eqs. (A28)-(A30).From these equations it follows that the ageostrophic terms in the equation for the synoptic-scale kinetic energy crucially depend on atmospheric friction, which primarily acts at the surface and within the PBL.This is conceptually in line with the atmospheric energy conversion cycle presented in Lorenz (1967).According to Lorenz (1967),

Fig. 1 .
Fig. 1.Band-pass filtered ERA-40 data of terms capturing exchange of kinetic energy between eddies and the large scale flow (in m 2 s −3 ) arising from (a) zonal, (b) meridional and (c) vertical shear of the large scale wind field at 500 mb for January.Note the different scale in (c).

Fig. 10 .
Fig. 10.Sensitivity of E k (at 500 mb, December-January-February averaged) to parameters (a) τ E,min , (b) K H and K 2,z , (c) KPBL , (d) m, (e) U 0 and (f) K meso .Plotted is the change in E k when the value of the parameter is changed from minus 10 % to plus 10 % the value given in Table2, thus covering a range of 20 %.

Table 1 .
Description of used variables.
• S and 60 • S is well represented.Both the local minimum at the southern tip of South America (due to the influence of the Andes) and the area of maximum activity stretching from the region south of Africa till the south-western tip of Australia are reasonably resolved.At low altitude, the model overestimates E k , with the highest values at 850 mb nearly double those in the ERA-40 data.A reason for this could be the presence of Antarctic slope winds which could reduce synoptic activity but which are not captured by our model.However, also the large-scale wind field data used as input is likely less accurate in near-

Table 2 .
Description of used constants.
Prandtl-type description of K syn in terms of the Austausch coefficient.More specifically, K syn is treated as the product of the characteristic synoptic vertical velocity ( W ) and vertical scale ( H ):