Mathematical and Physical Ideas for Climate Science
Valerio Lucarini, Richard Blender, Corentin Herbert, Salvatore Pascale, Francesco Ragone, Jeroen Wouters
MMathematical and Physical Ideas for Climate Science
Valerio Lucarini, ∗ Richard Blender, Salvatore Pascale, Francesco Ragone, † and Jeroen Wouters ‡ Klimacampus, Meteorologisches Institut,University of Hamburg,20144 Grindelberg 5, Hamburg,Germany.
Corentin Herbert
National Center for Atmospheric Research,P.O. Box 3000, Boulder, CO, 80307,USA.
The climate is a forced and dissipative nonlinear system featuring non-trivial dynamics ofa vast range of spatial and temporal scales. The understanding of the climate’s structuraland multiscale properties is crucial for the provision of a unifying picture of its dynamicsand for the implementation of accurate and efficient numerical models. We presentsome recent developments at the intersection between climate science, mathematics, andphysics, which may prove fruitful in the direction of constructing a more comprehensiveaccount of climate dynamics. We describe the Nambu formulation of fluid dynamics,and the potential of such a theory for constructing sophisticated numerical models ofgeophysical fluids. Then, we focus on the statistical mechanics of quasi-equilibriumflows in a rotating environment, which seems crucial for constructing a robust theoryof geophysical turbulence. We then discuss ideas and methods suited for approachingdirectly the non-equilibrium nature of the climate system. First, we describe somerecent findings on the thermodynamics of climate and characterize its energy and entropybudgets, and discuss related methods for intercomparing climate models and for studyingtipping points. These ideas can also create a common ground between geophysics andastrophysics by suggesting general tools for studying exoplanetary atmospheres. Weconclude by focusing on non-equilibrium statistical mechanics, which allows for a unifiedframing of problems as different as the climate response to forcings, the effect of alteringthe boundary conditions or the coupling between geophysical flows, and the derivationof parametrizations for numerical models.
I. INTRODUCTION
The Earth’s Climate provides an outstanding exam-ple of a high-dimensional forced and dissipative complexsystem. The dynamics of such system is chaotic, so thatthere is only a limited time-horizon for skillful prediction,and is non-trivial on a vast range of spatial and temporalscales, as a result of the different physical and chemi-cal properties of the various components of the climatesystem and of their coupling mechanisms (
Peixoto andOort , 1992).Thus, it is extremely challenging to construct satis-factory theories of climate dynamics and it is virtuallyimpossible to develop numerical models able to describeaccurately climatic processes over all scales. Typically,different classes of models and different phenomenologi-cal theories have been and are still being developed by ∗ [email protected]; Also at: Department of Mathe-matics and Statistics, University of Reading, Reading, RG6 6AX,UK. † Now at: Klimacampus, Institut f¨ur Meereskunde, University ofHamburg, 20146 Bundestrasse 53, Hamburg, Germany. ‡ Now at: Laboratoire de Physique, ´Ecole Normale Sup´erieure deLyon, 69364 46, all´ee d’Italie, Lyon, France. focusing on specific scales of motion (
Holton , 2004;
Val-lis , 2006), and simplified parametrizations are developedfor taking into account at least approximately what can-not be directly represented (
Palmer and Williams , 2009).As a result of our limited understanding of and abilityto represent the dynamics of the climate system, it is hardto predict accurately its response to perturbations, werethey changes in the opacity of the atmosphere, in the so-lar irradiance, in the position of continents, in the orbitalparameters, which have been present for our planet dur-ing all epochs (
Saltzman , 2001). The full understandingof slow- and fast-onset climatic extremes, such as droughtand flood events, respectively, and the assessment of theprocesses behind tipping points responsible for the multistability of the climate system are also far from beingaccomplished.Such limitations are extremely relevant for problemsof paleoclimatological relevance such as the onset anddecay of ice ages or of snowball-conditions, for contingentissues like anthropogenic global warming, as well as in theperspective of developing a comprehensive knowledge onthe dynamics and thermodynamics of general planetaryatmospheres, which seems a major scientific challenge ofthe coming years, given the extraordinary developmentof our abilities to observe exoplanets (
Dvorak , 2008). a r X i v : . [ phy s i c s . a o - ph ] A ug Climate science at large has always been extremelyactive in taking advantage of advances in basic mathe-matical and physical sciences, and, in turn, in provid-ing stimulations for addressing new fundamental prob-lems. The most prominent cases of such interaction arerelated to the development of stochastic and chaotic dy-namical systems, time series analysis, extreme value the-ory, radiative transfer, and fluid dynamics, among oth-ers. At this regard, one must note that the year 2013has seen a multitude of initiatives all around the worlddedicated to the theme
Mathematics of Planet Earth (see http://mpe2013.org ), and, in this context, climate-related activities have been of great relevance.In this review we wish to present some interdisciplinaryresearch lines at the intersection between climate science,physics and mathematics, which are extremely promisingfor advancing, on one side, our ability to understand andmodel climate dynamics, and represent correctly climatevariability and climate response to forcings. On the otherside, the topics presented here provide examples of howproblems of climatic relevance may pave the way for new,wide-ranging investigations of more general nature.The literature related to the scientific interface men-tioned above is enormous, and the selection of the mate-rial we present here is partial and non-exhaustive. Weleave almost entirely out of this review very impor-tant topics such as extreme value theory (
Ghil et al. ,2011), multiscale techniques (
Klein , 2010), adjoint meth-ods and data assimilation (
Wunsch , 2012), partial differ-ent equations (
Cullen , 2006), linear and nonlinear sta-bility analysis (
Vallis , 2006), general circulation of theatmosphere (
Schneider , 2006), macroturbulence (
Love-joy and Schertzer , 2013), networks theory (
Donges et al. ,2009), and many relevant applications of dynamical sys-tems theory to geophysical fluid dynamical problems (
Di-jkstra , 2013;
Kalnay , 2003).Let us now mention what we are going to cover inthis review and give a motivation to the specific perspec-tive we have chosen. We are motivated by the desireof bridging the gap between some extremely relevant re-sults in mathematical physics, statistical mechanics, andtheoretical physics, and open problems and issues of cli-mate science, hoping to stimulate further investigationsand interdisciplinary activities. Our selection of topicswill focus on the concepts of energy, entropy, symmetry,coupling, fluctuations, and response.We will first concentrate on the properties of invis-cid and unforced flows relevant for geophysical fluid dy-namics (GFD). In section II, we provide an overview ofa very powerful formulation of hydrodynamics based onthe formalism introduced by
Nambu (1973) and presentits applications in a geophysical context, suggesting howthese ideas help clarifying somewhat hidden propertiesof fluid flows, and how the Nambu formulation of GFDcould lead to a new generation of numerical models, tobe used in a variety of weather and climate applications In section III, starting from the classical investigation by
Onsager (1949) of the dynamics of point vortices, we willshow how to develop an equilibrium statistical mechani-cal theory of turbulence for GFD flows and will discuss itsrelevance for interpreting observed climatic phenomena.Equilibrium methods allow investigating many proper-ties of GFD flows. Nonetheless, at this point we cannotignore anymore the elephant in the room , i.e. , the factthat the dynamics of the climate system cannot be as-similated to an inviscid and unforced GFD flow, becauseforcing and dissipative processes are of extreme relevance.Thus, we move towards the paradigm of non-equilibriumsystems. In section IV, taking inspiration from the pointsof view of Prigogine (1961) and of
Lorenz (1967), we ex-plore how through classical non-equilibrium thermody-namics one can construct tools for assessing the energybudget and transport of the climate system, define andestimate the efficiency of the climate machine , and studythe irreversible processes by evaluation the climatic ma-terial entropy production. This allows for characterizingthe large scale properties of climate, for developing toolsfor auditing climate models, for gathering information ontipping points, and for exploring the properties of gen-eral planetary atmospheres. In section V, we addressthe non-equilibrium statistical mechanics formulation ofclimate dynamics, and explore how the formalism of re-sponse theory allows for addressing in a rigorous frame-work the climatic response to perturbations, taking in-spiration from the work of
Ruelle (1997). We will showhow it is possible to construct operators useful for theprediction - in an ensemble sense - of climate change. Alast aspect of GFD we want to discuss in a statisticalmechanical setting is the derivation of parametrizationsproviding a surrogate description of the effect of fast,small scale variables, which are hard to represent explic-itly in numerical models, on the larger scale, slow vari-ables of more direct climatic relevance. Thus, in sectionVI, we present averaging and homogenization techniques,describe how projector operator methods due to
Mori (1965) and
Zwanzig (1961) provide powerful tools for de-riving parametrizations and firm ground to the inclusionof stochastic terms and memory effects, and discuss howresponse theory can be used to derive similar results.Finally, in section VII we draw our conclusions andpresent some perspectives of future research.
II. BEYOND THE HAMILTONIAN PARADIGM: NAMBUREPRESENTATION OF GEOPHYSICAL FLUIDDYNAMICS
Hamiltonian formalism constitutes the backbone ofmost physical theories. In the case of a discrete au-tonomous system, the basic idea is to provide a full de-scription of the degrees of freedom by defining a set ofcanonical variables q and of the related momenta p ( q , p ∈ R N , i.e. , they are N -dimensional vectors), and byidentifying the time evolution to a flow in phase spacesuch that the canonical Hamiltonian function H acts asa streamfunction, ˙ q = ∇ p H , ˙ p = −∇ q H , where H ( q, p )corresponds to the energy of the system, whose value isconstant in time. The flow is inherently divergence-free( solenoidal ), so that the phase space does not contractnor expands, as implied by the Liouville Theorem ( Lan-dau and Lifshits , 1996). The time evolution of any func-tion X ( q, p ) can be expressed as: ddt X = ˙ X = { X, H} P = ∇ q X · ∇ p H − ∇ p X · ∇ q H , (1)where { , } P are the so-called Poisson brackets and · indi-cates the usual scalar product. As suggested by Noether’stheorem, the presence of symmetries in the system im-plies the existence of so-called physically conserved quan-tities X i , such that ˙ X i = 0 = { X i , H} P . An autonomoussystem possesses time invariance and its energy is con-stant, while in a system possessing translational invari-ance, the total momentum M is also constant. A systemcan possess many constants of motions, called Casimirs ,apart from energy, but the Hamiltonian plays a specialrole as it is the only function of phase space appearing explicitly in the definition of the evolution of the system(
Landau and Lifshits , 1996).
Nambu (1973) presented a generalization of canonicalHamiltonian theory for discrete systems. The dynamicalequations are constructed in order to satisfy Liouville’sTheorem and are written in terms of two or more con-served quantities. The Nambu approach has been ex-tremely influential in various fields of mathematics andphysics and is viable to extension to the case of con-tinuum, so that it can be translated into a field the-ory. The construction of a Nambu field theory for geo-physical fluid dynamics went through two decisive steps.The first was the discovery of a Nambu representationof 2D and 3D incompressible hydrodynamics (
N´evir andBlender , 1993). The second important step was the find-ing that the Nambu representation can be used to designconservative numerical algorithms in geophysical models,and that classical heuristic methods devised by Arakawafor constructing accurate numerical models actually re-flected deep symmetries coming from the Nambu struc-ture of the underlying dynamics of the flow (
Salmon ,2005).The physical basis for the relevance of the Nambu the-ory for describing and simulating conservative geophysi-cal fluid dynamics comes from the existence of relevantconserved quantities apart energy when forcing and dis-sipative terms are disregarded from the evolution equa-tions. Such a property is found in several models rel-evant for studying geophysical flows, and are valid for2D and 3D hydrodynamics, Rayleigh-B´enard convection,quasi-geostrophy, shallow water model, and extends tothe fully baroclinic 3D atmosphere. In other terms, the Nambu representation provides the natural descriptionof geophysical fluid dynamics and is superior to the moretraditional approaches based essentially on Euler equa-tions, just like the action-angle representation of the dy-namics of a spring is superior to the simple descriptionprovided by the second Newton’s law of motion.
A. Hydrodynamics in 2D and 3D
In incompressible hydrodynamics enstrophy (in 2D)and helicity (3D) are known as integral conserved quan-tities besides energy (
Kuroda , 1991).
N´evir and Blender (1993) adapted Nambu’s formalism to incompressiblenonviscous hydrodynamics by using enstrophy and he-licity in the dynamical equations.
1. Two-dimensional hydrodynamics
The evolution of two-dimensional incompressible invis-cid and unforced flows described by the velocity field u is governed by the vorticity equation ∂ω∂t = ∂ t ω = − u · ∇ ω (2)where customary symbols are used for indicating partialderivatives, the vorticity ω can be expressed, in Cartesiancoordinates ( x, y ), as ω = v x − u y , and incompressibilityis described by ∇ · u = 0, where ∇ · U = ∂ x U x + ∂ y U y isthe divergence of the vector field U . As a result, we canwrite u = S ∇ ψ = ( − ∂ y ψ, ∂ x ψ ) , where S is the sym-plectic matrix [0 , −
1; 1 , ψ is the streamfunction, and ∇ φ = ( ∂ x φ, ∂ y φ ) is the gradient of the function φ . Notethat ω = ∇ ψ . In this section, we consider a compactdomain ( e.g. , a square of side L ) with periodic boundaryconditions.The Hamiltonian H is the kinetic energy H = 12 (cid:90) u d A = − (cid:90) ωψ d A, (3)and is a functional of velocity. In general, a functional F [ φ ] maps a function φ of the phase space into a num-ber. The functional derivative δ F /δφ the change of thefunctional F with respect to a change in the function φ .The functional derivative can be defined by consideringthe first term in the expansion F [ φ + δφ ] − F [ φ ] = δ F [ φ ] = (cid:90) δ F δφ ( x ) δφ ( x )d x + . . . (4)The functional derivative δ H /δω for (3) is explicitlycalculated by δ H = (cid:90) ∇ ψ · δ ∇ ψ d A = (cid:90) ∇ · ( ψδ ∇ ψ ) d A − (cid:90) ψδω d A Since the first integral vanishes due to the boundary con-ditions, and since ω = ∇ ψ , we obtain δ H /δω = − ψ .Equation (2) says that vorticity is transported acrossthe domain by a non-divergent flow. One can prove easilythat any functional of the vorticity is conserved C = (cid:90) s ( ω ) d A. (5)where the integration is performed over the whole domainof the system. The most familiar of such functional is thetotal enstrophy of the flow: E = 12 (cid:90) ω d A (6)The functional derivative of the enstrophy is simply δ E /δω = ω .Since u = S ∇ ψ = ( − ∂ y ψ, ∂ x ψ ), the 2D vorticity equa-tion can be expressed as ∂ω∂t = −J ( ψ, ω ) = −J (cid:18) δ E δω , δ H δω (cid:19) , (7)with the antisymmetric Jacobi operator J ( a, b ) = ∂ x a ∂ y b − ∂ y a ∂ x b = −J ( b, a ) . (8)Relating ψ and ω to the functional derivatives of twoconserved quantities amounts to expressing the evolutionequation in a Nambu form using the enstrophy E .The time-evolution of an arbitrary functional of vor-ticity F = F [ ω ] is determined by d F dt = − (cid:90) δ F ∂ω J (cid:18) δ E δω , δ H δω (cid:19) d A = {F , E , H} (9)which defines a Nambu bracket for the three functionalsinvolved. The bracket is anti-symmetric in all arguments, {E , H , F} = −{H , E , F} , etc. Using rearrangements ofthese functionals and partial integration it can be shownthat the Nambu bracket is cyclic {F , E , H} = {E , H , F} = {H , F , E} (10)The cyclicity of this bracket is a main ingredient inSalmon’s application of Nambu mechanics ( Salmon ,2005) to construct conservative numerical codes (see Sec-tion II.B.2).In the following the relationship between Nambu me-chanics and Hamiltonian theory of two-dimensional flowsis briefly summarized. As mentioned above, a Hamilto-nian description of the dynamics is obtained when we canwrite d F dt = {F , H} P (11)with an antisymmetric Poisson bracket, to be seen in gen-eral as an antisymmetric map in the space of functionals, such that {A , B} P = −{B , A} P . Deriving such a bracketamounts to defining the dynamics of the system.The Poisson bracket for 2D hydrodynamics ( Salmon ,1988;
Shepherd , 1990) is easily obtained from the Nambubracket if the dependency δ E /δω = ω is evaluated {F , H} P = {F , E , H} = (cid:90) ω J ( F ω , H ω ) d A (12)where we indicate H ω = δH/δω ; here cyclicity is used,see Eq. 10.The Poisson bracket used in Eulerian hydrodynamics isdegenerate because of the presence of an infinite numberof so-called Casimirs , i.e. , the functionals defined in Eq.5, which are automatically conserved so that {C , H } P =0. In this case, we talk about noncanonical Hamiltonianmechanics.The relationship (12) demonstrates that noncanonicalHamiltonian mechanics is embedded in Nambu mechan-ics. The main extension is that in Nambu mechanics twofunctionals acting as an Hamiltonian, the enstrophy andthe energy, are used (7), and that the Nambu bracket (9)is nondegenerate and void of Casimir functionals.
2. Three-dimensional incompressible hydrodynamics
The dynamics of incompressible unforced and inviscidfluid flows in three dimension is determined by the vor-ticity ω = ∇ × u evolution equation: ∂ ω ∂t = ω · ∇ u − u · ∇ ω (13)where u is the velocity field and ∇ · u = 0. Note that incartesian coordinates we have that the curl of U ( ∇ × U )can be expressed as ( ∇× U ) i = (cid:15) ijk ∂ j U k where (cid:15) ijk is thestandard totally antisymmetric Levi-Civita symbol and ∇ · U = ∂ x U x + ∂ y U + ∂ z U z is the divergence in threedimensions. Similarly to the two-dimensional case, thetotal energy H = 12 (cid:90) u d V = − (cid:90) ω · A d V (14)is conserved, where we have introduced A as the vectorpotential such that u = −∇ × A . Note that in deriv-ing the second identity we use integration by parts andconsider periodic boundary conditions. It is importantto note that the total helicity h = 12 (cid:90) ω · u d V (15)is also conserved, while e.g. , the enstrophy is not. Fol-lowing the procedure detailed in Eq. 4, we derive thatthe functional derivative of the energy with respect tothe vorticity is given by δ H /δ ω = − A and for helicity δh/δ ω = u (compare the 2D version (5)).The Nambu form of the vorticity equation is ∂ ω ∂t = K (cid:18) δhδ ω , δ H δ ω (cid:19) = − K ( u , A ) (16)with K ( U , U ) = −∇ × [( ∇ × U ) × ( ∇ × U )] (17)Considering that ω = ∇ × u and using some standardvector calculus identities, we obtain that Eq. 16 agreeswith Eq. 13. We can derive the evolution equations forfunctional F = F [ ω ] as follows: d F dt = − (cid:90) (cid:18) ∇ × δ F δ ω (cid:19) × (cid:18) ∇ × δhδ ω (cid:19) · (cid:18) ∇ × δ H δ ω (cid:19) d V = {F , h, H} (18)where the last equation defines the Nambu bracket for3D incompressible hydrodynamics based on the vorticityequation. Helicity is no longer a hidden conserved quan-tity but enters the dynamics on the same level as theHamiltonian. Therefore, the Nambu mechanics is ableto account explicitly for conservation laws of the system,and, correspondingly, to its symmetries. B. Geophysical fluid dynamics
A Nambu representation can be constructed also forsome of the most important mathematical models rele-vant for geophysical fluid dynamics on large scales: thequasi-geostrophic potential vorticity equation (
N´evir andSommer , 2009), the shallow water model (
Salmon , 2005;
Sommer and N´evir , 2009), and the baroclinic stratifiedatmosphere (
N´evir and Sommer , 2009). Other modelsof geophysical relevance can also be treated in this way,as, most notably, the Rayleigh-B´enard equations for two-dimensional convection, which have been studied in detailin
Bihlo (2008) and
Salazar and Kurgansky (2010). Wewill not treat this latter case in this review.
1. Quasi-geostrophic approximation
Quasi-geostrophic (QG) theory is one of the most im-portant and most studied pieces of geophysical fluid dy-namics and is of crucial relevance for studying the large-scale dynamics of the Earth’s atmosphere and ocean, and,more recently, of planetary atmospheres (
Holton , 2004;
Klein , 2010;
Pedlosky , 1987). QG dynamics is relevantwhen, within a good approximation, the fluid motions are1) hydrostatic and 2) the Coriolis acceleration balancesthe horizontal pressure gradients. This is typically real-ized, e.g. , in the atmospheric midlatitudes. In absenceof dissipative processes and of forcings, QG dynamics isdescribed by the material conservation of the QG poten-tial vorticity. We consider customary Cartesian coordi-nates plus time ( x, y, z, t ), where x indicates the zonal direction, y the meridional direction, and z the verticaldirection as defined by gravity as in Holton (2004). Theevolution equation reads as follows: ∂Q∂t + 1 f J (Φ , Q ) = 0 (19)where J is the Jacobian (8). Q is the QG approximationof Ertel’s potential vorticity Q = ω g + f N ∂ Φ ∂z + f (20)with the geostrophic vorticity ω g = 1 /f ∇ h Φ, geopoten-tial Φ, ∇ h is the Laplacian operator limited to the x and y directions, Brunt-V¨ais¨al¨a frequency N , and Cori-olis parameter f = f + βy . The geostrophic veloc-ity u g has nonzero components only along the x and y directions, so that we can write u g = ( u hg , u hg = 1 /f S ∇ h Φ = 1 /f ( − ∂ y Φ , ∂ x Φ), where ∇ h is thegradient operator limited to the x and y directions.The first conserved integral is the total energy of thesystem H = 12 (cid:90) (cid:34)(cid:18) ∇ h Φ f (cid:19) + (cid:18) N ∂ Φ ∂z (cid:19) (cid:35) d V (21)where the first term is the density of kinetic energy andthe second term is the density of potential energy. Ateach level z the geopotential acts as a stream function indefining the geostrophic velocity field, while the verticalderivative of the geopotential is proportional to the tem-perature fluctuations of the system ( Holton , 2004). Thesecond conserved integral is the potential enstrophy E = 12 (cid:90) Q d V (22)which is defined similarly to the enstrophy in Eq. 6. Onecan prove that QG dynamics can be written in a Nambuform as follows: ∂Q∂t = −J (cid:18) δ E δQ , δ H δQ (cid:19) (23)Thus, the mathematical structure is analogous to thetwo-dimensional vorticity equation (9). Moreover, wecan construct the evolution of any functional F [ Q ] bydefining the Nambu bracket as follows: d F dt = − (cid:90) δ F δQ J (cid:18) δ E δQ , δ H δQ (cid:19) d V = {F , E , H} (24)with δ E /δQ = Q and δ H /δQ = − Φ /f .
2. Shallow water model
Roughly speaking, shallow water equations are usefultwo-dimensional approximations of Navier-Stokes equa-tions often used for describing some fluid motions wherethe horizontal scale of motion is much larger than its ver-tical extent, such as in the case of tidal waves or tsunamiin the ocean, or Rossby and Kelvin waves in the atmo-sphere. Here the single layer model is summarized (
Som-mer and N´evir , 2009). The dynamics is given by theevolution of the vorticity ω and the divergence µ = ∇ · u of the horizontal velocity u ∂ t ω = −∇ · ( ω u ) (25) ∂ t µ = k · ∇ × ( ω a u ) + ∇ ( u / gh T ) (26) ∂ t h T = −∇ · ( h u ) (27)where h T is the total height of the fluid. The shallowwater model possesses two conserved integrals, the totalenergy, given by the sum of kinetic and potential energy H = 12 (cid:90) ρ (cid:0) h T u + gh T (cid:1) d A (28)and potential enstrophy E = 12 (cid:90) ρq d A (29)with the absolute potential vorticity q = ω a /h T , ω a = ω + f . The functional derivatives of the conserved inte-grals are δ H /δω = − ρψ , δ H /δµ = − ργ , δ H /δh T = ρ Ψ, δ E /δω = ρq , δ E /δµ = 0, δ E /δh T = − (1 / ρq , where ρ isthe density of the fluid, ψ is the streamfunction, γ the ve-locity potential for h v = S ∇ ψ + ∇ γ , Ψ = (1 / u + gh T is the Bernoulli function,The Nambu representation of the shallow water modelwas derived by Salmon (2005) and is a bit more cumber-some than in, e.g. , QG case.
Sommer and N´evir (2009)present a numerical simulation of these equations on aspherical grid, and
N´evir and Sommer (2009) publishedthe multilayer shallow water equations. In the case of asingle layer shallow water equations, the dynamics of anyfunctional F is determined by the sum of three Nambubrackets ddt F = {F , H , E} ω,ω,ω + {F , H , E} µ,µ,ω + {F , H , E} ω,µ,h T (30)The first bracket is {F , , H , E} ω,ω,ω = (cid:90) J ( F ω , , H ω ) E ω d A (31)where X ω = δ X /δω . Such first bracket is analogous tothe 2D Nambu bracket (9) (apart from the sign). Forthe other brackets we refer to ( Salmon , 2005;
Sommerand N´evir , 2009).
Salmon (2007) calculated the Nambubrackets based on the velocities instead of vorticity.
3. Baroclinic atmosphere
N´evir and Sommer (2009) published the equations de-termining the dynamics of a baroclinic dry atmospherein Nambu form (denoted as Energy-vorticity theory ofideal fluid mechanics). The Nambu representation en-compasses the Eulerian equation of motion in a rotatingframe, the continuity equation, and the first law of ther-modynamics. The Nambu dynamics uses three brack-ets for energy, helicity, energy-mass, and energy-entropy.Due to its special role in all three brackets the integral ofErtel’s potential enstrophy is coined as a super-Casimir.The Nambu form shows an elegant structure wherefundamental processes are combined by additive terms.Incompressible, barotropic or baroclinic atmospheres areassociated to additive contributions. Thus approxima-tions are simply attained by the neglect of terms.In absence of forcings and of dissipative processes, themomentum equation, the continuity equation and thefirst law of thermodynamics equation are (
Peixoto andOort , 1992) ∂ t u = − u · ∇ u − Ω × u − ρ ∇ p − ∇ Φ (32) ∂ t ρ = −∇ · ( ρ u ) (33) ∂ t s = − u · ∇ s (34)where u is velocity, Ω the angular velocity of the earth, Φis the sum of the gravitational and centrifugal potentialof the earth, ρ is density and s is the specific entropyper unit mass, determined by the equation of state of thegas.These equations possess four conservation laws. Thefirst is the total energy H = (cid:90) ρe d V ; e = 12 u + i + Φ (35)where e is the specific total enery and i is its internalenergy component. The absolute helicity is h a = (cid:90) u a · ω a d V (36)where the absolute velocity is u a = u + Ω × r and ω a = ∇ × v + 2 Ω . with the angular velocity of the earth Ω ,and r is the position vector. The total mass and entropyare given by M = (cid:90) ρ d V, S = (cid:90) ρs d V (37)and the total potential enstrophy is defined starting fromErtel’s potential vorticity Π E ρ = (cid:90) ρ Π d V, Π = ω a · ∇ θρ , (38)analogously to the definition in the QG context given inEq. 22. The functional derivatives of the conservationlaws are δ H /δ u = ρ u , δ H /δρ = (1 / u + i + p/ρ − T s + φ , δ H /δσ = T , δ M /δ u = 0, δ M /δρ = 1, δ M /δσ =0, δ S /δ u = 0, δ S /δρ = s , δ S /δσ = 1, δh a /δ u = ω a , δh a /δρ = 0, and δh a /δσ = 0, where T is temperature,and σ = ρs .An arbitrary functional F of u , ρ and σ evolves ac-cording to the sum of three brackets which are definedbelow ddt F = {F , h a , H} h + {F , M, H} m + {F , S, H} s (39)The three brackets are defined below. The first one isthe so-called helicity bracket: {F , h a , H} h = − (cid:90) (cid:20) ρ δ F δ u · (cid:18) δh a δ u × δ H δ u (cid:19)(cid:21) d V ; (40)the second is the so-called mass bracket: {F , M , H} m = − (cid:90) (cid:20) δ M δρ δ F δ u · ∇ δ H δρ + δ F δρ ∇ · (cid:18) δ M δρ δ H δ u (cid:19)(cid:21) d V + cyc ( F , M , H ); (41)where cyc indicates permutations in cyclic order of the ar-guments. The third one is the so-called entropy bracket: {F , S , H} s = − (cid:90) (cid:20) δ S δρ δ F δ u · ∇ δ H δσ + δ F δσ ∇ · (cid:18) δ S δρ δ H δ u (cid:19)(cid:21) d V + cyc ( F , M , H ) (42)For a barotropic flow the first law of thermodynamics isphysically not relevant and the entropy bracket is dis-carded in (39) because the functional derivatives withrespect to σ vanish. The continuity equation remains un-approximated and the pressure gradient term is replacedby the gradient of enthalpy. Note the different bracketsfor helicity (40) and vorticity (18) in 3D hydrodynamics. C. Conservative algorithms and numerical models
Salmon (2005, 2007) recognized that the existence ofa Nambu bracket with two conserved integrals allows thedesign of high-precision numerical algorithms for study-ing geophysical flows. The idea is in fact simple: just likein the usual case we aim at writing numerical codes ableto conserve energy when dissipation and forcing are ne-glected, Nambu mechanism provides encouragement andconceptual support for expanding this point of view byencompassing other important physical quantities. Theapproach is useful in GFD turbulence simulations be-cause these flows are characterized by the existence ofconservation laws besides total energy. In particular, the conservation of enstrophy inhibits spurious accumulationof energy at small scales.For the numerical design of conservative codes basedon a Nambu structure the following remarks are noted: • A Nambu form of the continuous physical systemis required. • The quantities used in the Nambu bracket are con-served. • The discrete form of the Jacobian needs to preserveantisymmetry (11). • The approach is applicable to any kind of dis-cretization, e.g. for finite differences, finite vol-umes, or spectral models. • Arbitrary approximations of the conservation lawsare possible; these approximations are conservedexactly. • For the barotropic vorticity equation the classicArakawa Jacobian could be retrieved by equallyweighting the cyclic permutations of the Nambubracket. In other terms, Arakawa found heuristi-cally a discrete Nambu representation of barotropicdynamics (
Dubinkina and Frank , 2007).In recent years, various authors have provided promisingexamples of actual implementations of GFD codes whichtake into explicit consideration the underlying Nambudynamics of the unforced and inviscid case.
Salmon (2007) presents the first numerical simulation of a shallowwater model derived from the Nambu brackets formalism.The simulation is on a square rectangular grid and thedesign on an unstructured triangular mesh is outlined.
Sommer and N´evir (2009) report the first simulationof a shallow water atmosphere using Nambu brackets.The authors use an isosahedric grid (as in the ICONmodel, ICOsahedric Non-hydrostatic model, of the Ger-man Weather Service and the Max Planck Institute forMeteorology, Hamburg). The construction of the algo-rithm is as follows (
Sommer and N´evir , 2009):1. First the continuous versions of the Nambu-brackets and conservation laws need to be obtained.2. On the grid, the following expressions need to becalculated: functional derivatives, discrete opera-tors (div and curl), discretization of the Jacobianand the Nambu brackets.3. Finally, the prognostic equations are obtained byinserting the variables in the brackets. Various op-tions are available for the time stepping is arbi-trary;
Sommer and N´evir (2009) use a leap-frogwith Robert-Asselin filter.
FIG. 1 Enstrophy tendencies in the enstrophy conservingICON shallow water model and the Nambu model of
Sommerand N´evir (2009) (courtesy of Matthias Sommer, Ludwig-Maximilians-Universit¨at M¨unchen). Note that the tendencyin the Nambu model is of the order of the numerical accuracy.
The authors find quasi-constant enstrophy and energycompared to a standard numerical design (Fig. 1).Along these lines,
Gassmann and Herzog (2008) sug-gest a radically new concept for a global numericalsimulation of the non-hydrostatic atmosphere using theNambu representation for the energy-helicity bracket {F , h a , H} ( N´evir , 1998). Their suggestion incorporatesa careful description of Reynolds averaged subscale pro-cesses and budgets.
Gassmann (2013) describes a globalnon-hydrostatic dynamical core based on an icosahedralnonhydrostatic model on a hexagonal C-grid. The modelconserves mass and energy in a noncanonical Hamil-tonian framework,even if some still unsolved numericalproblems occur when the non-hydrostatic compressibleequations are in a Nambu bracket form. The use of dy-namical cores constructed according to the sophisticatedversion of fluid dynamics discussed here might providecrucial for improving the ability of atmospheric modelsin representing correctly the global budgets of physicallyrelevant quantities also in the case when forcing and dis-sipative processes are taken into account. As discussedby
Lucarini and Ragone (2011) for the case of energy,this is far from being a trivial task.
D. Perspectives
Like Hamiltonian mechanics, the Nambu approach isa versatile tool for the analysis and simulation of dynam-ical systems. Here some possible research directions areoutlined.
Modular modeling and approximations:
In sev-eral applications a Nambu representation can be found byadding brackets which conserve a particular Casimir, thisis already mentioned by
Nambu (1973); see the baroclinicatmosphere (
N´evir and Sommer , 2009) and the classifi- cation by
Salazar and Kurgansky (2010). The dynam-ics is determined by these ’constitutive’ Casimirs (a no-tion coined in (
N´evir and Sommer , 2009)) which are notconserved in the complete system. Decomposition leadsto subsystems where the constitutive Casimirs are con-served. An example is helicity which is constitutive inthe baroclinic atmosphere and only conserved in the 3Dincompressible flow. The decomposition is directly asso-ciated with approximations (
N´evir and Sommer , 2009).Composition allows a process-oriented model design.
Statistical Mechanics:
The statistical mechanics offluids is characterized by the existence of conservationlaws besides total energy (
Bouchet and Venaille , 2012),see also section III in this review. Thus these conserva-tion laws have a two-fold impact: They determine thedynamics in a Nambu bracket and the canonical proba-bility distribution in equilibrium.
Dynamics of Casimirs:
Casimir-functions of a con-servative system are ideal observables to characterize thedynamics in the presence of forcing and dissipation. Thismight prove especially interesting when studying the re-sponse of a system to perturbations in the context ofthe Response theory proposed by (
Ruelle , 1997, 1998a,b,2009) and recently used in a geophysical context by var-ious authors with promising results (
Abramov and Ma-jda , 2008;
Eyink et al. , 2004;
Lucarini , 2009;
Lucariniand Sarno , 2011); see also section V in this review.As illuminating example, we mention the recent workof
Pelino and Maimone (2007) and
Gianfelice et al. (2012), who have used recurrence maps of extremes ofenergy and a Casimir in a Lorenz-like map to assess pre-dictability of the system and study the properties of theinvariant measure.
III. EQUILIBRIUM STATISTICAL MECHANICS FORGEOPHYSICAL FLOWS
We have seen in the previous section that differentmodels of geophysical flows have a specific mathemat-ical structure: they are Hamiltonian systems, and havean infinite number of conserved quantities - the
Casimirs .The previous section has shown how one could take ad-vantage of these features and construct theoretically richrepresentation of the dynamics and provide proposals forconstructing new numerical codes of GFD flows. Thissection goes in the direction of constructing a probabilis-tic description of GFD flows, basically taking the point ofview that due to the large amount of degrees of freedominvolved, one can consider the state of the atmosphereand the ocean as random variables. Here we shall reviewthe progress that has been made by using the simplestclass of possible probability distributions: the equilib-rium distributions which depend only on the conservedquantities. However, most of the standard applications ofequilibrium statistical mechanics deal with dynamics ona finite dimensional phase space (e.g., a gas with a finitenumber of molecules), with a finite number of dynamicalinvariants (often just the energy). The equations describ-ing the dynamics of geophysical flows violate both theseconstraints. Several solutions have thus been proposed:they are reviewed briefly in the next sections, going fromthe main fundamental ideas to selected geophysical ap-plications.
A. Finite-dimensional models: Point vortices
1. Negative temperature states and clustering of vortices
Onsager (1949) was the first to understand that the co-herent structures and persistent circulations that appearubiquitously in planetary atmospheres and in the Earth’soceans could be explained on statistical grounds. Hiswork focused on 2D incompressible, inviscid fluids givenin equation 2. To make the system tractable, he intro-duced an approximation of the vorticity field in terms of N point vortices with circulation γ i and position r i ( t ): ω ( r , t ) = (cid:80) Ni =1 γ i δ ( r i ( t ) − r ), where δ ( x ) is the usualDirac’s delta distribution. Introducing the Hamiltonian H = − (cid:80) i Yatsuyanagi et al. , 2005). For negative tempera-tures, we observe the clustering of same-sign vortices, whilefor positive temperatures, positive and negative vortices aredistributed homogeneously in the domain. 2. Mean-field equation The above argument is qualitative; to characterize thecoherent structures which are expected to emerge fromthe clustering of same-sign vortices, we introduce theprobability density ρ i ( r , t ) for a vortex with strength γ i to be found at point r at time t . It satisfies the normaliza-tion (cid:82) ρ i ( r , t ) d r = 1. We define a coarse grained vorticityfield ω ( r , t ) = (cid:80) i γ i ρ i ( r , t ). This probability density isexpected to converge towards its statistical equilibrium:the equilibrium distribution maximizes the statistical en-tropy S = − (cid:80) i (cid:82) ρ i ( r ) ln ρ i ( r ) d r . The solution of thisvariational problem is given by ρ i ( r ) = e β ( γ i ψ ( r )+ µ i ) / Z where β and βµ i are the Lagrange parameters associ-ated with conservation of global energy and normaliza-tion of each ρ i , respectively, and ψ = ∆ − ω is the coarse-grained stream function, while the normalization factor Z is called the partition function . Averaging over thisequilibrium distribution gives the coarse-grained vortic-ity field, which satisfies the mean-field equation : ω ( r ) = 1 Z (cid:88) i γ i e β ( γ i ψ ( r )+ µ i ) . (45)This is an equation of the form ω = F ( ψ ), characteris-tic of the steady-states of the 2D Euler equations. Awell-known particular case is that of N vortices withcirculation 1 /N and N vortices with circulation − /N .In that case, the mean-field equation can be recast as ω = A sinh( β Ψ), with Ψ = ψ − ( µ + − µ − ) / Montgomeryand Joyce , 1974).0The theory can be generalized in a straightforwardmanner to quasi-geostrophic (QG) flows ( Miyazaki et al. ,2011). DiBattista and Majda (2001) have given solu-tions of the mean-field equation for a two-layer model- i.e. , a QG model where the stream function is de-fined only at two discrete value of the vertical coordi-nate and the temperature is defined at the interface be-tween such level ( Holton , 2004) - where the point vor-tices stand for hetons , introduced by Hogg and Stommel (1985) as a model of individual convective towers in theocean. They have shown that a background barotropiccurrent (the barotropic governor ) confines potential vor-ticity and temperature anomalies, thereby suppressingthe baroclinic instability, in agreement with numericalsimulations ( Legg and Marshall , 1993).The point vortex model suffers from a number of lim-itations inherent to the approach. First of all, when welet the number of vortices tend to infinity (the thermo-dynamic limit ), we have to introduce an ad-hoc scalingof the Lagrange parameters to retain the organized, neg-ative temperature states. Besides, there is no uniqueway to approximate a vortex patch by a finite number ofvortices. A consequence is also that the area of vorticitypatches cannot be conserved in this singular formulation.We shall see in section III.C that dealing directly withthe vorticity field will solve these issues, while predict-ing a relation between vorticity and stream function verysimilar to the one obtained above. B. Finite-dimensional models: truncated Fourier modes 1. 2D Turbulence Rather than a discretization in physical space, one mayconsider a finite number of modes in Fourier space, asproposed by Lee (1952) and Kraichnan (1967) in thecontext of the Euler equations. For 2D flows — for sim-plicity, we consider here a rectangular geometry with pe-riodic boundary conditions; the case of a spherical ge-ometry can be found in Frederiksen and Sawford (1980)— writing the vorticity field as a truncated Fourier series ω ( x ) = (cid:80) k ˆ ω ( k ) e i k · x , the evolution in time of the Fouriercoefficients follows an equation of the form ∂ t ˆ ω ( k ) = (cid:80) p , q A kpq ˆ ω ( p )ˆ ω ( q ), where the summation is restrictedto a finite set of wave vectors B = { k ∈ π/L Z , k min ≤ k ≤ k max } and A kpq takes care of the quadratic non-linearity terms. This dynamics preserves two quadraticquantities: the energy E = (cid:80) k | ˆ ω ( k ) | / (2 k ) and the en-strophy Γ = (cid:80) k | ˆ ω ( k ) | . Kraichnan (1967) suggested toconsider the canonical probability distribution: ρ ( { ˆ ω ( k ) } k ∈B ) = e − ˜ βE − α Γ Z , (46) In particular, the average energy at absolute equilibrium is given by (cid:104) E (cid:105) = − ∂ ln Z ∂ ˜ β = 12 (cid:88) k ∈B β + 2 αk , (47)which corresponds to an equipartition spectrum for thegeneral invariant ˜ βE + α Γ : E ( k ) = πk/ ( ˜ β + 2 αk ). In-viscid numerical runs indeed relax to this spectrum ( Bas-devant and Sadourny , 1975; Fox and Orszag , 1973). Notethat the Lagrange parameters α and ˜ β cannot take ar-bitrary values; they are constrained by the realizability condition — for the Gaussian integral defining Z to con-verge. Here, this condition reads: ˜ β + 2 αk > β + 2 αk > 0. In particular, when α > 0, negativetemperatures can be attained. In this regime, which cor-responds to (cid:104) Γ (cid:105) / (2 (cid:104) E (cid:105) ) small enough ( Kraichnan andMontgomery , 1980), the energy spectrum is a decreasingfunction of k . When ˜ β → − αk , a singularity appearsat k → k min , which means that the energy is expected toconcentrate in the largest scales. Hence, statistical me-chanics for the truncated system predicts that when theenstrophy is small enough compared to the energy, weexpect the energy to be transferred to the large scales. Kraichnan (1967) gives other arguments to support andrefine this view; in particular he shows the existence oftwo inertial ranges, with a constant flux of energy andenstrophy, respectively, with the energy spectrum scal-ing as E ( k ) ∼ Cε / k − / and E ( k ) ∼ C (cid:48) η / k − re-spectively, where ε and η are the energy and enstrophyfluxes. In particular, the equilibrium energy spectrum atlarge scales is shallower than the energy inertial rangespectrum. Assuming a tendency for the system to relaxto equilibrium — although the equilibrium is never at-tained in the presence of forcing and dissipation — wethus expect the flux of energy to be towards the largescales; a process referred to as the inverse cascade of 2Dturbulence. Similarly, the transfer of enstrophy in thecorresponding inertial range should be towards the smallscales. The dual cascade scenario has been confirmedboth by numerical simulations ( Boffetta , 2007) and lab-oratory experiments ( Paret and Tabeling , 1997). 2. Quasi-Geostrophic Turbulence The dynamical equations of QG flow are very similarto the Euler equations, replacing vorticity by potentialvorticity (see section II.B.1). In particular, they con-serve similar quadratic invariants, and the theory can beextended in a straightforward manner ( Holloway , 1986; Salmon , 1998). We will discuss in this section the effectof stratification and β effect.Perhaps the simplest framework to consider the roleof stratification is the two-layer QG case. As in sectionIII.B.1, a canonical probability distribution can be con-structed, taking into account the three invariants: the1total energy E and the potential enstrophies of eachlayer, Z and Z . The corresponding partition functioncan be computed, and the spectrum studied in the var-ious regimes, with similar results. In particular, nega-tive temperature states are accessible, which correspondto condensation of the energy on the largest horizon-tal scales and the Fofonoff (1954) solutions mentionedbelow. Maybe more interestingly, although the vari-ous forms of energy (kinetic energy K , K in each layerand potential energy P ) are not individually conserved,we can compute their average value at equilibrium, as Salmon et al. (1976) did. Alternatively, the standarddecomposition in terms of the barotropic and baroclinicmodes (constructed by taking the average and the differ-ence of the streamfunctions in the two layers), with theirkinetic energies K T and K B , can be used. As Salmonet al. (1976) highlighted, the Rossby deformation scale k D = 2 π/R D plays an important role. R D = N H/f ,where H is the vertical extent of the domain and f isthe reference Coriolis parameter, defines the typical hor-izontal scale of perturbations of vertical extent equal to H , with N/f (cid:29) k (cid:29) k D ), the two layers behave essentially as two in-dependent copies of 2D turbulence; the energy spectrumin each layer is the same as in the 2D case, the correla-tion at statistical equilibrium is low, and there is aboutas much energy in the barotropic mode and the baro-clinic mode: (cid:104) K T ( k ) (cid:105) / (cid:104) K B ( k ) (cid:105) ∼ 1. Besides, the po-tential energy is small compared to the kinetic energy: (cid:104) P ( k ) (cid:105) / (cid:104) K T ( k ) (cid:105) = O ( k D /k ). At scales larger than thedeformation radius ( k (cid:28) k D ), the system rather behavesas a unique barotropic layer: the amount of energy inthe two layers is about the same, but the energy is es-sentially in the barotropic mode, with negligible energyin the baroclinic mode, and a statistical correlation be-tween the two layers of order 1. This theoretical anal-ysis goes in strong support of the standard picture oftwo-layer QG turbulence, developed on phenomenolog-ical grounds ( Rhines (1979); Salmon (1978), see also Vallis (2006, chap. 9) and Fig. 3), and is in agree-ment with numerical simulations ( Rhines , 1976). Theseresults have been extended to an arbitrary number oflayers and to continuously stratified flows by Merryfield (1998). Although the equilibrium mean, vertically inte-grated stream function remains similar to the two-layercase, the distribution of the statistics on the vertical dif-fers as higher-order moments are considered. The ratioof potential to kinetic energy for instance, can becomesignificantly underestimated, especially in the limit ofstrong stratification ( k D → 0) where the two-layer modeldoes not capture well the possibility that an importantfraction of the energy may be trapped near the bottom.The second dominant effect in geophysical flows, inaddition to stratification, is rotation. The Coriolis forceintroduces a linear term in the equations, which does not Baroclinic EnergyBarotropic EnergyWind orSolar InputLoss toBoundaryLayer Friction Scattering Into3D Turbulence k k D k D FIG. 3 Energy (solid arrows) and potential enstrophy (dashedarrows) flux diagram for two-layer quasi-geostrophic turbu-lence, taking inspiration from Salmon (1978). The energyinjected in the baroclinic mode at large scales is cascadeddownscale until the deformation scale is reached, then it istransferred to the barotropic mode and cascaded upscale likein 2D turbulence, in agreement with the predictions of equi-librium statistical mechanics. affect directly the previous analysis of the nonlinear en-ergy transfers: the conserved quantities remain the sameand the statistical theory is easily extended by replacingrelative vorticity with absolute vorticity. However, thevariation of the Coriolis force with latitude is responsi-ble for the appearance of Rossby waves, which modify thephysical interpretation of the predicted cascade of energy.As anticipated by Rhines (1975) and verified numerically(e.g., Vallis and Maltrud , 1993), the Rossby waves deflectthe inverse energy cascade: they dominate over nonlineareffects in a part of Fourier space and prevent access tolow wavenumbers along one direction in Fourier space.This leads to the preferential formation of zonal flows. 3. Beyond balanced motion Although the large-scale motions of the atmosphereand oceans of the Earth are very close to geostrophicand hydrostatic balance, these relations break up whenmoving down to the mesoscale, and the transfers of en-ergy due to turbulence, or the non-linear interaction ofinertia-gravity waves, might not follow the inverse cas-cade scenario described in sections III.B.1-III.B.2. As amatter of fact, a downscale transfer of energy is neededin the ocean to feed enhanced vertical mixing (e.g., Led-well et al. , 2000) or small-scale dissipation in the oceaninterior ( Nikurashin et al. , 2013). Such processes are nec-essary to close the energy budget of the ocean ( Wunschand Ferrari , 2004). It is therefore natural to ask howequilibrium statistical mechanics can help understand-ing how energy is exchanged by nonlinear interactionsbetween the slow, balanced motions and the fast, wavemotions. Errico (1984) first observed a tendency for unforced in-viscid flows described by hydrostatic primitive equationsto reach an energy equipartition state, in which the en-ergy in the fast wave modes is comparable to that in theslow balanced modes. The study by Warn (1986), in the2context of the shallow water equations, essentially con-firms that QG flows are not equilibrium states, and thata substantial part of the energy may end up in the fast(surface) wave modes at statistical equilibrium, implyinga direct cascade of energy to the small scales. Bartello (1995) has obtained analytically the equilibrium energyspectrum for the Boussinesq equations (neglecting thenonlinear part of potential vorticity), in the presence ofrotation, confirming the direct cascade of energy. In par-ticular, there is no negative temperature states in thiscase, due to the presence of the inertia-gravity waves. Infact, numerical simulations ( Pouquet and Marino , 2013)indicate that turbulence with rotation and stratificationmight have at the same time an inverse and a direct cas-cade of energy. A natural interpretation would be thatvortical modes are responsible for the inverse cascadewhile waves cascade energy downscale simultaneously. Bartello (1995) had discussed the possibility of a wave-vortical mode decoupling on the basis of resonant triadicinteractions. Without any assumptions on the dynamics,another interpretation in the statistical mechanics frame-work uses an analogy with metastable states: restrictingthe equilibrium probability distribution to the slow man-ifold yields an inverse cascade, while taking into accountthe whole phase space including the waves results in adirect cascade ( Herbert et al. , 2014). C. The mean-field theory for the continuous vorticity field 1. Mean-field theory Above, we have considered finite-dimensional modelsconserving at most two quadratic quantities, the energyand the enstrophy. In fact, the majority of the flows con-sidered above — and in particular 2D and QG flows —conserve an infinite family of invariants, called Casimirinvariants: for any function s , (cid:82) s ( ω ) d r is conserved (seeEq. 5). The specific case s n ( x ) = x n corresponds to themoments of the vorticity distribution. Instead, the con-servation of s σ ( x ) = δ ( x − σ ) implies that the area γ ( σ )where the vorticity takes value σ is conserved. This is dueto the absence of a vortex stretching term, in contrastwith full 3D flows; here the vorticity (or potential vortic-ity, in the QG case) patches are stirred in such a way thattheir area remains conserved. The theory developed by Miller (1990) and Robert and Sommeria (1991) (see also Bouchet and Venaille (2012) for a review) introduces acoarse-grained vorticity field ω , which corresponds to themacroscopic state of the flow. This coarse-grained vor-ticity field can be predicted based on the invariants usingstatistical mechanics. To do so, we introduce ρ ( σ, r ), theprobability density for the vorticity field to take value σ at point r . The coarse-grained vorticity field is given by ω ( r ) = (cid:82) ∞−∞ σρ ( σ, r ) dσ . The invariants of the system are the energy E [ ρ ] = (cid:90) D d r d r (cid:48) (cid:90) R dσdσ (cid:48) σσ (cid:48) G ( r , r (cid:48) ) ρ ( σ, r ) ρ ( σ (cid:48) , r (cid:48) ) , (48)with G the Green function of the Laplacian, and theCasimir invariants G n [ ρ ] = (cid:90) D d r (cid:90) R dσσ n ρ ( σ, r ) , (49)or equivalently, the vorticity levels D σ [ ρ ] = (cid:90) D d r ρ ( σ, r ) . (50)The idea of the theory is to select the probability dis-tribution ρ which maximizes a mixing entropy S [ ρ ] = − (cid:82) D (cid:82) R d r dσρ ( σ, r ) ln ρ ( σ, r ), under the constraints ofconservation of the invariants, and pointwise normaliza-tion N [ ρ ]( r ) = (cid:82) R dσρ ( σ, r ) = 1. Hence, we are inter-ested in the variational problem:max ρ, N [ ρ ]( r )=1 { S [ ρ ] | E [ ρ ] = E, ∀ n ∈ N , G n [ ρ ] = Γ n } , (51)or equivalently,max ρ, N [ ρ ]( r )=1 { S [ ρ ] | E [ ρ ] = E, ∀ σ ∈ R , D σ [ ρ ] = γ ( σ ) } . (52)The solutions of this variational problem correspond tothe most probable states for a given set of conservedquantities.The critical points of the variational problem (52)are simply given by δ S − (cid:82) d r ζ ( r ) δ N ( r ) − ˜˜ βδ E − (cid:82) dσα ( σ ) δ D σ = 0, where ˜ β and α ( σ ) are the Lagrangemultiplier associated with the conservation constraints.Easy computations yield the solution ρ ( σ, r ) = 1 Z e ˜ βσψ ( r ) − α ( σ ) , (53)so that the coarse-grained vorticity is given by ω = F ( ψ ) , with F ( ψ ) = 1˜ β δ ln Z δψ , (54)and Z ( ψ ) = (cid:82) R dσe ˜ βσψ − α ( σ ) . To compute the equilib-rium states of the system, one should solve the partialdifferential equation (54), referred to as the mean-fieldequation , and check afterwards that the obtained criticalpoints are indeed maxima of the constrained variationalproblem by considering the second derivatives. This willautomatically ensure that the equilibrium states are non-linearly stable steady states ( Chavanis , 2009).3 FIG. 4 Maximum entropy states as a function of the aspectratio for a rectangular domain, in the linear ( strong mixing ) ω - ψ limit ( Chavanis and Sommeria , 1996). For τ < τ c , theequilibrium is a monopole, while for τ > τ c , it is a dipole. 2. Equilibrium states for 2D and barotropic flows The mean-field equation (54) is in general difficult tosolve; one issue is that the ω - ψ relation is in generalnonlinear. Most of the analytical solutions have been ob-tained in the linear case, by decomposing the fields ona basis of eigenfunctions of the Laplacian on the domain D . This technique was first introduced in a rectangulardomain by Chavanis and Sommeria (1996), who showedthat the statistical equilibrium is either a monopole or adipole, depending on the aspect ratio (Fig. 4). Thesame method was extended to the case of barotropicflows, replacing vorticity by potential vorticity. Tak-ing into account the β effect, Fofonoff (1954) flows areobtained as statistical equilibria in a rectangular basin( Naso et al. , 2011; Venaille and Bouchet , 2011). Such so-lutions correspond to flows with two gyres (anticyclonicin the northern basin, cyclonic in the southern basin) ina rectangular basin, see Fig. 5. The relative vorticity isconfined to a boundary layer, whose with decreases withthe total energy or when the β effect ( i.e. , the relativestrength of the gradient of the planetary vorticity) in-creases. The flow is westward in the interior of the basin,with an eastward compensating flow near the boundaries.Different geometries can be studied: in a rotatingsphere, the equilibria, in the linear limit, can be eithersolid-body rotations, dipole flows ( Herbert et al. , 2012)or quadrupoles, taking into account conservation of an-gular momentum ( Herbert , 2013). In the latter case, aperturbative treatment of the nonlinearity in the ω - ψ relationship leads to the same flow topology, but sharpervortex cores ( Qi and Marston , 2014). Bouchet and Si-monnet (2009) have also considered the role of a smallnonlinearity in the ω - ψ relationship for a rectangular domain of aspect ratio close to 1, with periodic boundaryconditions, thereby obtaining two topologies for the equi-librium states: dipole and unidirectional flows. Addinga small stochastic forcing generates transitions from oneto the other equilibrium. 3. Stratified flows In addition to the 2D and quasi-2D cases mentionedabove, the theory has also been applied to stratified flu-ids (essentially in the quasi-geostrophic regime). Herbert (2014) has obtained and classified the statistical equi-libria of the two-layer QG model in the framework ofthe Robert-Miller-Sommeria theory, and updated the dis-cussion of the vertical distribution of energy at statisti-cal equilibrium (see section III.B.2): in particular, it isshown that even at statistical equilibrium, there will re-main some residual energy in the baroclinic mode, unlessthe initial vertical profile of fine-grained enstrophy is uni-form. In the context of continuously stratified flows, Ve-naille (2012) has taken up the thread initiated by Merry-field (1998) (see section III.B.2) and shown that bottom-trapped currents are indeed statistical equilibria of theRobert-Miller-Sommeria theory. Still in the continuouscase, Venaille et al. (2012) have also studied the verticaldistribution of energy at statistical equilibrium, focusingon the tendency to reach barotropic equilibrium states;as also observed in the two-layer model, the constraintof conservation of fine-grained enstrophy prevents com-plete elimination of energy in the baroclinic mode. Asthe β effect increases, barotropization is facilitated, untilwe enter a regime dominated by waves. It is well knownthat baroclinic dynamics is hindered by strong values β ( Holton , 2004). D. Subgrid scale parameterization Results from equilibrium statistical mechanics havefound practical applications in the development of pa-rameterization methods. Holloway (1992) suggested toreplace the usual sub-grid scale parameterizations inocean models, where, e.g. , viscous forces are representedwith terms of the form ν ∗ ∆ u ), where ν ∗ is the eddyviscosity. He proposed to replace such formula with ν ∗ ∆( u − u ∗ ), so that viscosity relaxes the system towardsthe statistical equilibrium state u ∗ . Such a parameteri-zation has been implemented, tested and commented in anumber of studies ( Cummins and Holloway , 1994, e.g.).For more perspective on this type of subgrid-scale param-eterizations, the reader is referred to ( Holloway , 2004)and ( Frederiksen and O’Kane , 2008).Along similar lines, Kazantsev et al. (1998) have pro-posed more generally to treat the subgrid scales so asto maximize the entropy production, inspired by the4 Journal of Marine Research a ----~:~:::~~~~~ ------------------------- izz:.-;- ________----------------------- _.-- --.. _A .- -____________-------------------------- ._ --._ /- ___---- _____--------------________ ..- --. __.- --.. -.._ _--- --- ________----------------------- __--- --__ ---_ _____-- _________---------------------------------- b d 2 -.6 Figure 8. Experiment using the smoothed average 5 field in ISRDl instead of random field as initial condition. Now the field has very little transient component. The fields are averaged from 10 T,, and then smoothed. Shown here are (a) relative vorticity 5, (b) absolute vorticity q, (c) streamfunction IJJ and (d) scatter-plot of q - I). consistent with the theory. Roughly measuring from the scatter-plots, we obtain Apl CL- -0.37 where E is the average of p’s in the two experiments, A~J, is the increase in t,r, from 32 x 32 to 128 x 128 experiments. Therefore, we have: AE 3Ak -- --- E - - 2 CL - 55-5%* Both energies are the components contained in the time mean fields. Hence, when R, = 3.18 x 10e3 in the initial random flows, then as resolution is increased from FIG. 5 Convergence towards the statistical equilibrium in inviscid truncated barotropic flow on a β -plane ( Wang and Vallis ,1994). Left: Stream function. Right: scatterplot of the q - ψ relation. relaxation equations formulated in the Robert-Miller-Sommeria theory as an algorithm to construct equilib-rium states ( Chavanis and Sommeria , 1997). Note alsothat it has been shown in direct numerical simulationsof ideal 3D turbulence that the small scales thermalizeprogressively, and act as a sort of effective viscosity inthe ideal system, leading to the appearance of transientKolmogorov scaling laws ( Cichowlas et al. , 2005). Thisseems to be consistent with the above suggestions forsubgrid scale parameterizations. IV. CLIMATE AS A FORCED-DISSIPATIVETHERMODYNAMIC SYSTEM In the previous sections the focus has been on identify-ing symmetry properties and conservation laws of GFDflows and relate these to dynamical features and statis-tical mechanical properties. Neglecting forcing and dis-sipation has led us to study reversible equations whosestatistical properties can be described using equilibriumstatistical mechanics.Indeed, this provides the backbone of the properties ofGFD flows and are of great relevance for studying morerealistic physical conditions. Nonetheless, at this stage,a reality-check is necessary. The atmosphere and theoceans are out-of-equilibrium systems, which exchangeirreversibly matter and energy from their surroundingenvironment and re-export it in a more degraded form athigher entropy. For example, Earth absorbs short-waveradiation (low-entropy solar photons emitted at a tem-perature of ≈ ≈ 255 K). In addition to that,spatial gradients in chemical concentrations and temper-ature as well as their associated internal matter and en-ergy fluxes can be established and maintained for longtime within non-equilibrium systems (e.g. the temper- ature contrast between the polar and equatorial regionsand the associated large-scale, atmospheric and oceaniccirculation). In this and in the next sections we will takesuch a point of view.The basis of the physical theory of climate was estab-lished in a seminal paper by ( Lorenz , 1955), who eluci-dated how the mechanisms of energy forcing, conversionand dissipation are related to the general circulation ofthe atmosphere. Oceanic and atmospheric large scaleflows results from the conversion of available potentialenergy - coming from the differential heating due to theinhomogeneity of the absorption of solar radiation- intokinetic energy through different mechanisms of instabil-ity due to the presence of large temperature gradients( Charney , 1947; Eady , 1949). Such instabilities createa negative feedback, as they tend to reduce the tem-perature gradients they feed upon by favoring the mix-ing between masses of fluids at different temperatures.Furthermore, in a forced and dissipative system like theEarth’s climate, entropy is continuously produced by irre-versible processes ( deGroot and Mazur , 1984; Prigogine ,1961). Contributions to the total material entropy pro-duction , which is related to the non-radiative irreversibleprocesses ( Goody , 2000; Kleidon , 2009), come from: dis-sipation of kinetic energy due to viscous processes, tur-bulent diffusion of heat and chemical species, irreversiblephase transitions associated to various processes relevantfor the hydrological cycle, and chemical reactions relevantfor the biogeochemistry of the planet.It is important to note that the study of the climate en-tropics has been revitalized after Paltridge (1975, 1978)proposed a principle of maximum entropy production(MEPP) as a constraint on the climate system. Whilethe scientific community disagrees on the validity of sucha point of view - see, e.g. , Goody (2007) - the discussionrevolving around MEPP has led the scientific communityto refocus on the importance of a thermodynamical ap-5proach – as complementary to the dynamical one – inproviding physical insights for student the climate sys-tem. In this paper we will not discuss MEPP and othernon-equilibrium variational principles (for an updated re-view see Dewar et al. , 2013). A. Climatic energy budget and energy flows 1. Energy Budget We first focus on developing equations describing theenergy budget of the climate system. The total specific(per unit mass) energy of a geophysical fluid is given bythe sum of internal, potential, kinetic and latent energy.This can be expressed as e = u / i + φ + Lq for the at-mosphere, where u is the velocity vector, i = c v T is thenertanl energy, with c v is the specific heat at constantvolume for the gaseous atmospheric mixture and T is itstemperature, Φ is the gravitational (plus centrifugal) po-tential, L is the latent heat of evaporation, and q is thespecific humidity. In this formula, we neglect the heatcontent of the liquid and solid water and the heat asso-ciated to the phase transition between solid and liquidwater. The approximate expression for the specific en-ergy of the ocean reads e = u / i + Φ, where i = c W T is the specific heat at constant volume of water (we ne-glect the effects of salinity and of pressure), while we canconsider e = c S T + φ as the specific energy for solid earthor ice. The conservation of energy and the conservationof mass imply that ( Peixoto and Oort , 1992): ∂ρe∂t = −∇ · ( J h + F R + F S + F L ) − ∇ ( τ · u ) (55)where ρ is the density; p is the pressure; J h = ( ρe + p ) u is the total enthalpy transport; F R , F S , and F L are thevectors of the radiative, turbulent sensible, and turbulentlatent heat fluxes, respectively; and τ is the stress tensor.By expressing Eq. (55) in spherical coordinates ( r, λ, ϕ ),and assuming the usual thin shell approximation r = R + z , z/R (cid:28) 1, where R is the Earth’s radius and z isthe vertical coordinate of the fluid, we have ( Peixoto andOort , 1992): [ ˙ E ] = − R cos ϕ ∂T T ∂ϕ + [ F T OAR ] (56)where [ X ]( ϕ, t ) ≡ (cid:82) X ( λ, ϕ, t ) dλ , F T OAR is the net radi-ation at the top of the atmosphere (with the conventionthat the value is positive when there is an excess of incom-ing over outgoing radiation) and the meridional enthalpytransport has been defined as: T T ( ϕ, t ) ≡ (cid:90) (cid:90) J hϕ ( ϕ, λ, z, t ) R cos ϕdzdλ. (57)Equation (56) relates the rate of change of the verticallyand zonally integrated total energy to the divergence of the meridional transport by the atmosphere and oceansand the zonally integrated radiative budget at the top-of-the-atmosphere. Integrating along ϕ ( { X } = (cid:82) Xdϕ ),the expression for the time derivative of the net globalenergy balance is straightforwardly derived: { [ F T OAR ] } = { [ ˙ E ] } . (58)Similar relationships can be written for the atmosphere, below. Uncertainties are discussed in Trenberth et al. (2009) and are mostly not dealt withhere, except to highlight some sources of discrepancy with other studies.In Fig. 1, use has been made of conservation of energy and the assumption that, on atime scale of years, the change in heat storage within the atmosphere is very small.Accordingly, the net radiation at TOA R T is the sum of the ASR minus the OLR: R T = ASR - OLR. In turn, the ASR is the difference between the incoming solar radi-ation and the reflected solar radiation. At the surface, the ASR has to be offset by thesensible heat and latent heat fluxes plus the net longwave radiation. The latter is made up oftwo large terms: the emitted radiation from the surface and the downwelling longwaveradiation coming back from the atmosphere. Both at the surface and TOA the imbalance isthe same and, as noted above, is estimated to be 0.9 W m - .Updates in Trenberth et al. (2009) included revised absorption in the atmosphere bywater vapor and aerosols, since Kim and Ramanathan (2008) found that updated spec-troscopic parameters and continuum absorption for water vapor increased the absorption by4–6 W m - . The sensible heat has values of 17, 27, and 12 W m - for the globe, land, andocean (just over 70% of the Earth), and, even with uncertainties of 10%, the errors are onlyorder 2 W m - . There is widespread agreement that the global mean surface upwardlongwave (LW) radiation is about 396 W m - , which is dependent on the skin temperatureand surface emissivity (Zhang et al. 2006).Global precipitation should equal global evaporation for a long-term average, andestimates are likely more reliable of the former. However, there is considerable uncertaintyin precipitation over both the oceans and land (Trenberth et al. 2007b; Schlosser andHouser 2007). The latter is mainly due to wind effects, undercatch, and sampling, while the Fig. 1 The global annual mean earth’s energy budget for 2000–2005 (W m - ). The broad arrows indicatethe schematic flow of energy in proportion to their importance. Adapted from Trenberth et al. (2009) withchanges noted in the textSurv Geophys FIG. 6 Global annual mean earths energy budget for 2000-2005 ( W m − ). From Trenberth and Fasullo (2012). ocean and land provided that energy fluxes of sensible, la-tent heat as well as radiative fluxes are taken into accountat the surface ( Peixoto and Oort , 1992). A schematicview of the surface and TOA energy fluxes for presentday Earth ( Trenberth and Fasullo , 2012) can be see inFig. 6. Under steady state conditions, the long termaverage ˙ E = 0. Therefore from equation (58) the sta-tionarity condition implies that { [ F R ] toa } = 0 . (59)Equation (59) describes the basic fact that he climatesystem, at steady state, does not on the average receivesnor emits energy.These constraints can be used for auditing climatemodels. At observational level non-zero energy balancesare found at TOA and at the surface ( Trenberth and Fa-sullo , 2012; Wild et al. , 2013), due to the fact the theactual Earth is not at a stationary state, most notablybecause of the ongoing greenhouse gas forcing. However,a physically consistent climate model should feature avanishing net energy balance when its parameters areheld fixed and statistical stationarity is eventually ob-tained. Lucarini and Ragone − . − , with a few ones6 a)b)c)d)FIG. 7 Mean and standard deviation of globally averagedtop-of-the-atmosphere radiative budget (a), atmosphere en-ergy budget (b), ocean (c) and land energy budget (d) forinter-comparable CMIP3 (red) and CMIP5 (blue) climatemodels control simulations. Updated from Lucarini andRagone (2011). featuring imbalances larger than 3 W m − . The anal-ysis of similar budgets for the last generation of cli-mate models (CMIP5 intercomparison project, Tayloret al. (2012)) does not show a significant improvement(Fig. 7). Spurious energy biases may be associated withnon-conservation of water in the atmospheric branch ofthe hydrological cycle ( Liepert and Lo , 2013; Liepert and Previdi , 2012) and in the water surface fluxes ( Hassonet al. , 2013; Lucarini et al. , 2008), with the fact thatdissipated kinetic energy is not re-injected in the systemas thermal energy ( Becker , 2003; Lucarini and Fraedrich ,2009), as well as with nonconservative numerical schemes( Gassmann , 2013). 2. Meridional enthalpy transport The next step in constructing the energetics of the cli-mate system is the study of the large scale transport ofvarious forms of energy. The meridional distribution ofthe radiative fields at the top-of-the-atmosphere poses astrong constraint on the meridional general circulation( Stone , 1978). As clear from equation (56), the station-arity condition (59) leads to the following indirect rela-tionship for T T : T T ( ϕ ) = − π (cid:90) π/ ϕ R cos ϕ (cid:48) ( F R ) toa ( ϕ (cid:48) ) . (60)In other terms, the flux T T transports enthalpy from thelow-latitudes, which feature a positive imbalance betweenthe net input of solar radiation – determined by planetaryalbedo, determined mostly by i.e. , clouds ( Donohoe andBattisti , 2012) and by surface properties – and the outputof longwave radiation, to the high-latitudes, where a cor-responding negative imbalance is present. Atmosphericand oceanic circulations act as responses needed to equi-librate such an imbalance ( Peixoto and Oort , 1992).The climatic meridional enthalpy transport T T ( ϕ ) re-duces the temperature difference between the low andhigh latitude regions with respect to what imposed by theradiative-convective equilibrium picture. Stone (1978)showed that T T depends essentially on the mean plan-etary albedo and on the equator-to-pole contrast of theincoming solar forcing, while being mostly independentfrom dynamical details of atmospheric and oceanic cir-culations. As emphasized by Enderton and Marshall (2009), if one assumes drastic changes in the meridionaldistributions of planetary albedo differences emerge withrespect to Stone’s theory. A comprehensive thermody-namic theory of the climate system able to predict thepeak location and strength of the meridional transport,the partition between atmosphere and ocean ( Rose andFerreira , 2013) and to accommodate the variety of pro-cesses contributing to it, is still missing.Besides theoretical difficulties, observational estima-tions of T T , T A and T O also poses non-trivial challenges.For simplicity, we here refer to T T . There is still notan accurate estimate of such a fundamental quantity fortesting the output of climate models, despite the effortsof several authors ( Fasullo and Trenberth , 2008; Mayerand Haimberger , 2012; Trenberth and Caron , 2001; Tren-berth and Fasullo , 2010; Wunsch , 2005). The precisionof the estimates relies on the knowledge of the boundary F R , F S , and F L and on the reanalysis datasets. Wunsch (2005), by using measurements of the radiativefluxes at the top of the atmosphere and previous esti-mates of the oceanic enthalpy transport, gave a rangeof values of 3 . − . W) for themaximum of the total poleward transport in the North-ern Hemisphere (NH) and 4 . − . Trenberth and Fasullo (2010), by com-bining measurements of top-of-the-atmosphere radiativefields with different reanalyses and ocean datasets, foundthe range to be 4 . − . . − . Mayer and Haimberger (2012),using two reanalysis datasets (ERA-40 and the more re-cent ECMWF reanalysis ERA-Interim), constrained thetwo peaks in narrower confidence intervals: 5 . − . . − . . − . . − . . − at the top-of-the-atmosphere and+6 . − over oceans in ERA- Interim, ( Mayer and V OLUME 14J O U R N A L O F C L I M A T EF IG . 7. The required total heat transport from the TOA radiationRT is compared with the derived estimate of the adjusted ocean heattransport OT (dashed) and implied atmospheric transport AT fromNCEP reanalyses (PW). cient. This result is especially so at 24 N at which lat-itude no adjustments have been applied to either derivedOT estimate. Further, the adjustments applied south of30 S are miniscule for NCEP but amount to 0.7 PW at68 S for ECMWF OT or 0.4 PW for ECMWF AT (thedifference being integrated effects from the north vs thesouth), suggesting also that the latter are less reliablein absolute values. In the Tropics, the problems withchanges in the observing system, particularly satellitedata, adversely influence the ECMWF results (Trenberthet al. 2001b), which are not within the error bounds ofthe other estimates at several latitudes.Aside from the North Atlantic, for which the coupledmodel results are high for reasons that are beginning tobe understood (section 3c), the largest discrepancyamong the results is in the SH Tropics, where the NCEP-derived values imply a larger southward transport thando the direct ocean estimates or the coupled climatemodels. This is in a region where Ekman transports playa key role and surface wind specifications are very un-certain. For instance, at 11 S a change in mean windalters the direct estimate from 0.48 to 0.63 PW in theAtlantic (Holfort and Seidler 2001), and in coupled-model simulations tropical convergence zones are oftendislocated in some seasons when a spurious ITCZ formsin the SH, potentially corrupting values in the Tropics(Boville and Gent 1998).We have inferred the surface fluxes and thus the zonalmean ocean heat transports assuming no changes inocean heat storage except for those associated with glob-al warming. The local ocean heat storage change is notneglible from year to year (Sun and Trenberth 1998),although it is a reasonable assumption for zonal meansfor the four years or so we used here provided that theglobal warming trend is factored in, as we have done.Future computations should factor in the changes inocean heat storage to examine the ocean heat transportslocally using this method. Variability in AT from sam-pling this particular interval is mostly less than 0.05 PW and is not a major factor, although interannual variabilityis not very reproducible between ECMWF and NCEPreanalyses. TOA radiation fluxes also contain some un-certainty, and adjustments for expected imbalances atthe TOA may be a refinement worth considering in fu-ture.It is important to note that although the ocean heattransports and surface fluxes derived from the TOA ra-diation plus the atmospheric transports (the indirectmethod) have improved substantially and mostly agreewith the independent estimates the same cannot be saidfor the atmospheric NWP model surface fluxes com-puted with bulk parameterizations, which contain sub-stantial biases. The NWP models have not yet beenimproved to satisfy the global energy budgets in thesame way that the best coupled climate models have,highlighting the facts that weather prediction is con-strained by the specification of the SSTs and the modelsdo not have to get the SST tendencies correct to produceexcellent weather forecasts.Shortcomings in the hydrological cycle in the NCEPreanalyses in the Tropics (Trenberth and Guillemot1998) suggest that they have limitations, although, be-cause there is huge compensation between the budgetsfor dry static energy and the moist component, the totalenergy transport is more robustly computed (Trenberthand Solomon 1994). The discrepancies between the at-mospheric transports in the two reanalyses suggest thatfurther revisions will occur, especially regionally. Nev-ertheless, the results for the ocean heat transports de-rived from the NCEP reanalyses are in good agreementwith those from the other approaches, suggesting thatthe coupled models, the atmospheric transports, and theindependently estimated ocean transports are converg-ing to the correct values. Acknowledgments. This research was sponsored byNOAA Office of Global Programs Grant NA56GP0247and the joint NOAA–NASA Grant NA87GP0105. Manythanks to Dave Stepaniak for computing the atmosphericenergy transports. REFERENCESAagaard, K., and P. Greisman, 1975: Toward new mass and heatbudgets for the Arctic Ocean. J. Geophys. Res., J. Phys. Oceanogr., Deep-Sea Res. I, Mon. Wea.Rev., J. Climate, J. Cli-mate, N latitude. In-teractions between Global Climate Subsystems, The Legacy of FIG. 8 Annual meridional enthalpy transports of ocean(dashed), atmosphere (dash-dotted) and total (solid) esti-mated from satellite and reanalysis data ( P W ). From Tren-berth and Caron (2001). Haimberger , 2012)) problems that may potentially biasthe transport estimates. Furthermore, these estimatesare dependent on the analysis method and the modelused – Trenberth and Caron (2001), using other reanaly-sis dataset (NCEP), found a value of the maxima 0 . Trenberth and Caron (2001) areshown in Fig. 8.The use of numerical climate model does not help toreduce such uncertainties Lucarini and Ragone . Tayloret al. (2012)) show little improvement in terms of mu-tual agreement (Fig. 9). Donohoe and Battisti (2012)attributed such a large spread in T T to intermodel dif-ferences in the meridional contrast of absorbed solar ra-diation, which, in turn, is mainly due to the inter-modeldifference in the shortwave optical properties of the atmo-sphere; for an intercomparison of the cloud distribution isdifferent climate models see Probst et al. (2012). Figure9 also shows that, while the disagreement among modelsfor the peak of the atmospheric transport is compara-ble to that for the peak of the total transport, enormousdifferences emerge when comparing oceanic transports.Interesting information emerge when looking at the po-sition of the peaks of the transport. Stone (1978) pre-dicted that the position of the maximum of T T is wellconstrained by the geometry of the system and weakly de-pendent of longitudinal homogeneities, and, accordinglyin Fig. 9 both CMIP3 and CMIP5 models feature smallspread in the position of the peak of T T , with minute dif-ferences between the two hemispheres, except one outlier.Similarly, the spread among models is small with respectto the position of the peak of T A in both hemispheresand of T O in the Northern Hemispheres, while a largeruncertainty exists in the position of the peak of T O in theSouthern Hemisphere. B. The maintenance of thermodynamicaldisequilibrium The basic understanding of the maintenance of the at-mospheric general circulation was achieved nearly sixtyyears ago by Lorenz (1955, 1967) through the conceptsof available potential energy and atmospheric energy cy-cle. The concept of available potential energy, first intro-duced by Margules (1905) to study storms, is defined as A = (cid:82) c p ( T − T r ) dV , where T r is the temperature fieldof the reference state, obtained by an isentropic redis-tribution of the atmospheric mass so that the isentropicsurfaces become horizontal and the mass between the twoisentropes remains the same. By its own definition, thisstate minimizes the total potential energy at constantentropy. Such a definition is somewhat arbitrary anddifferent definitions lead to different formulations of at-mospheric energetics ( Tailleux , 2013). For example, thechoice of a reference state maximizing entropy at con-stant energy ( Dutton , 1973) leads in a natural way tothe concept of exergy. Exergy is the part of the internalenergy measuring the departure of the system from itsthermodynamic and mechanical equilibrium, i.e. ,a stateof maximum entropy at constant energy., and is a com-monly used concept in heat engines theory ( Rant , 1956). Lorenz (1967) proposed the following picture of thetransformation of energy in the atmosphere. We de-fine E = P + K , where K = (1 / (cid:82) dV ρ u repre-8 a)b)c)FIG. 9 Value and position of the peak of the poleward merid-ional enthalpy transport in the pre-industrial scenario for thewhole climate (a), atmosphere (b) and ocean (c) for the someof the CMIP3 (red) and CMIP5 (blue) general circulationmodels. Updated from Lucarini and Ragone (2011). sents the total kinetic energy and P = (cid:82) dV ρ ( c V T + Φ)the dry static energy and V is the atmospheric do-main. Under hydrostatic approximation one can showthat (cid:82) dzρ ( c v T + gz ) = (cid:82) dzρ ( c v T + RT ) = (cid:82) dzρc p T (see e.g. Lorenz , 1967).. In the Lorenz framework oneconsiders the hydrological cycle as a forcing to the at-mospheric circulation. This amounts to separating thebudget of the moist static energy and of the part relatedto the phase changes of water. See Peixoto and Oort (1992), Chap. 13. We obtain:˙ P = − W ( P, K ) + ˙Ψ + D, (61)˙ K = − D + W ( P, K ) , (62)where D = (cid:82) dV ρ(cid:15) > W ( P, K ) = − (cid:82) dV ρ u · ∇ p is the potential-to-kineticenergy conversion rate, and ˙Ψ = (cid:82) dV ρ ˙ q nf is the non-frictional diabatic heating due to the convergence of tur-bulent sensible heat fluxes, condensation/evaporation in-side the atmosphere, and convergence of radiative fluxes.The conversion term W can be interpreted as the in-stantaneous work performed by the system. In this re-spect, Eq. (61) represents the statement of the first lawof thermodynamics for the atmosphere. Equations (61)-(62) imply that ˙ E = ˙ P + ˙ K = ˙Ψ and therefore the fric-tional heating D does not increase the total energy sinceit is just an internal conversion between kinetic and po-tential energy. Stationarity implies that ˙ P = ˙ K = 0 andtherefore D = W , which is referred to as the intensityof the Lorenz energy cycle. One has to note that thelatter can be expressed as the average rate of generationof available potential energy, G = (cid:82) dV ρ ˙ q nf (1 − T /T r ),where T r is the temperature field of the reference state( Lorenz , 1967).The strength of the Lorenz energy cycle is a fundamen-tal non-equilibrium property of the atmosphere, which,just as the meridional enthalpy transport (Sect. IV.A.2),is known with a certain degree of uncertainty for thepresent climate. Reanalysis datasets (with all associ-ated problems, see Sect. IV.A.2) constrain D in the range1 . − . − ( Li et al. , 2008). On the other hand,general circulation models feature values of D rangingfrom 2 to 3 . − ( Marques et al. , 2011). Numericalsimulations show that a CO doubling causes a decreaseof G of nearly 10% ( Lucarini et al. , 2010a). Warm-ing patterns can alter G either by affecting the grossstatic stability (stronger stability implies a weaker en-ergy cycle, as clear from the theory of baroclinic insta-bility) or the meridional temperature/diabatic heatingdistribution. Hernandez-Deckers and von Storch (2012)show that the decrease in G is mostly associated withchanges in the gross static stability changes rather thanwith meridional temperature gradient changes.Another aspect to be considered is that the intensityof the Lorenz energy cycle is formulated assuming hy-drostatic conditions. Therefore, the Lorenz energy cy-cle in itself neglects any systematic transfer of potentialinto kinetic energy occurring through non-hydrostatic,small scale motions ( Steinheimer et al. , 2008). Alongthese lines, Pauluis and Dias (2012) suggest that smallscales processes such as precipitation may significantlycontribute to D , which might therefore be considerablyunderestimated when computed for models that do not9treat explicitly convection.In the case of the ocean, available potential energyis generated through thermohaline forcings due to thecorrelation of density inhomogeneities and density forc-ings (e.g. through heat and freshwater fluxes) at sur-face. In addition to that, mechanical energy enters theocean through direct transfer of kinetic energy by sur-face winds (and though tidal effects). Kinetic energy isdissipated through a variety of frictional processes, occur-ring mostly at the bottom of the ocean, and, similarly,available potential energy is lost through diffusion mostlydue to small scale eddies ( ? ). The understanding of thedetails of the oceanic Lorenz energy cycle is still at arelatively early stage. Estimates of dissipation and gen-eration terms range within (1 − × − W m − ( Oortet al. , 1994; Storch et al. , 2012; Tailleux , 2013). 1. Atmospheric heat engine and efficiency Johnson (2000) proposed an interesting constructionfor further elucidating the idea the the climate can beseen as a heat engine. We define the total diabatic heat-ing ˙ q = ˙ q nf + ρ(cid:15) and splitting the atmospheric domain V into the subdomain V + in which ˙ q = ˙ q + > , and V − ,where ˙ q = ˙ q − < 0, it can be seen from equation (61)that: W = (cid:90) V + ˙ q + ρdV + (cid:90) − ˙ q − ρdV ≡ Φ + + Φ − , (63)with Φ + > − < + and Φ − are the net heat gain and loss, and W themechanical work. The efficiency of the atmospheric heatengine, i.e. ,the capability of generating mechanical workgiven a certain heat input, can therefore be defined as: η = ( ˙Φ + + ˙Φ − )( ˙Φ + ) = W ˙Φ + . (64)The analogy between the atmosphere and a (Carnot)heat engine can be pushed further if we introducethe total rate of entropy change of the system, ˙ S = (cid:82) dV ρ ˙ q/T = ˙ S + + ˙ S − . In a steady state the followingexpression holds: ˙ S = ˙Φ + T + + ˙Φ − T − = 0 , (65)where T ± ≡ ˙Φ ± / (cid:82) V ± dV ρ ˙ q ± /T from which it followsthat η = 1 − T − /T + . Johnson’s approach provides aself-consistent treatment of the heat engine of a geophys-ical fluid and extends closely related thermodynamicaltheories of hurricane dynamics ( Emanuel , 1991).In Emanuel’s theory a mature hurricane is depicted asan ideal Carnot engine driven by the thermal disequi-librium between the sea-surface temperature T s and the cooling temperature T with an efficiency 1 − T /T s ≈ / 3. A similar approach was extended also to moist con-vection ( Emanuel and Bister , 1996; Renn`o and Ingersoll ,1996) for determining the the wind speed reached by theconvective system for a certain rate of heat input F in from the sea, W = F in (1 − T /T s ). Such an approachhas been used to study large scale, open systems like theHadley cell ( Adams and Renn`o , 2005) and the monsoonalcirculation ( Johnson , 1989). 2. Entropy production in the Climate System We wish now to emphasize a different aspect of theclimate’s thermodynamics, namely the study of its ir-reversibility by the investigation of its material entropyproduction, i.e. ,the entropy produced by the geophysi-cal fluid, neglecting the change in the properties of theradiative fields ( Goody , 2000; Ozawa et al. , 2003). Theentropy budget of the fluid can be rewritten as:˙ S = − (cid:90) dV ∇ · F R T + ˙ S mat , (66)so that we separate the contribution coming from the ab-sorption of the radiation from other effects related to theother irreversible processes occurring in the fluid. Notethat, in the previous formula, we refer to the entropybudget of the whole climate, not of the atmosphere, asdone, instead, in the previous section.The material entropy production, ˙ S mat can be ex-pressed as ˙ S diff + ˙ S fric + ˙ S hyd , i.e. , the sum of contri-butions associated with heat diffusion, frictional heatingand the hydrological cycle (due to diffusion of water andphase-changes) respectively. Detailed estimates of theentropy budget of the climate system and of the mate-rial entropy production ( ˙ S mat ≈ 50 mW m − K − ) canbe found in ( Goody , 2000; Pascale et al. , 2011). Oceanicentropy production due to small-scale mixing in the in-terior gives a small contribution ( ≈ − K − ) to˙ S mat ( Pascale et al. , 2011). Therefore we will limit thediscussion to processes occurring in the interior and atthe boundaries of the atmosphere.Entropy production due to heat diffusion ˙ S diff = − (cid:82) dV ∇ · J S /T is generally small ( ≈ − K − ,( Kleidon , 2009)) and associated mostly with dry atmo-spheric convection occurring nearby the surface and withvertical mixing in the mixed layer of the ocean. Theentropy production due to frictional heating - ˙ S fric = (cid:82) dV ρ(cid:15) /T ≈ 10 mW m − K − ( Fraedrich and Lunkeit ,2008; Pascale et al. , 2011) - is associated with turbu-lent energy cascades bringing kinetic energy from largescales down to scales (millimeters or less for geophys-ical flows) where viscosity can efficiently operate. Fi-nally, ˙ S hyd is due to irreversible processes associated withthe hydrological cycle – evaporation of liquid water in0unsaturated air, condensation of water vapor in super-saturated air and molecular diffusion of water vapour( Pauluis and Held , 2002a,b) and requires the knowledgeof relative humidity H and the molecular fluxes of watervapor J v = F L /L :˙ S hyd = (cid:82) dV ( C − E ) R (ln H + J v · ∇ p w ) − (cid:82) z = surf dAJ v,z R ln H. (67)where C and E indicate condensation and evaporation,respectively, and p w is the partial pressure of the watervapor. The importance of ˙ S diff and ˙ S hyd in the contextof thermodynamic theories of moist convection is exten-sively discussed in Pauluis and Held (2002b) and Pauluis (2010). The impact of water vapor on the production ofkinetic energy in deep convection can be described asa steam engine and it is to lower the maximum possi-ble amount of work which can produced by an equiva-lent Carnot cycle ( Emanuel and Bister , 1996; Renn`o andIngersoll , 1996) acting between the same temperaturereservoirs. An indirect estimate of (67) can be obtainedfrom the entropy budget for water˙ S W = (cid:90) V W dV ρ W ˙ s w = (cid:90) V W dV ρ W ˙ q W T + ˙ S hyd as discussed in Pauluis and Held (2002b), where ˙ S W isthe rate of change of entropy of water and ˙ q W the neatheating amount of heat per time that the water sub-stance receives from its environment (i.e. through evap-oration and condensation). At steady state ˙ S W = 0 andso ˙ S hyd = (cid:82) V W dV ρ W ˙ q W /T ≈ 37 mW m − K − ( Pas-cale et al. , 2011). Therefore, it is possible to computethe material entropy production by considering the ex-clusively heat exchanges and the temperature at whichsuch exchanges take place, thus bypassing the need forlooking into the complicated details of phase separationprocesses.Furthermore, in climate models aphysical entropysources due to diffusive/dispersive numerical advectionschemes and parameterizations are also present ( Egger ,1999; Johnson , 1997). In particular, Woollings andThuburn (2006) showed that dispersive dynamical corescan lead to negative numerical entropy production. Moregenerally, it has been argued that parameterizations ofsub-grid turbulent fluxes of heat, water vapor and mo-mentum should conform to the second law of thermo-dynamics, and therefore should lead to locally positivedefinite entropy production, this being generally not thecase ( Gassmann , 2013). C. Applications and future perspectives 1. Auditing Climate Models At steady state, we have that ˙ S = 0. Hence, from Eq.66 we derive: ˙ S mat = (cid:90) − dV ∇ · F R T . (68)Usually, this is referred to as indirect formula for com-puting the material entropy production ( Goody , 2000),because it provides an alternative way for estimating thematerial entropy production of the geophysical fluid byonly looking at the correlation between radiative heatingrates and temperature fields). Therefore, this formula al-lows for computing the material entropy production dueto fluid motions bypassing all the complex fluid dynam-ical behavior of the system. See Lucarini and Pascale (2014) for an in-depth discussion of different ways forcomputing the material entropy production and of theeffect of coarse graining the thermodynamic fields. Start-ing from Eq. (68) it is possible to derive for Earth condi-tions an approximate formula for the long term averageof the material entropy production, and to disentanglethe contributions due to horizontal and vertical processes( Lucarini et al. , 2011) as ˙ S mat ≈ ˙ S vmat + ˙ S hmat , where˙ S hmat = − (cid:90) A dA ∇ h · Υ T E = − (cid:90) A dA F T OAR T E (69)where Υ = (cid:82) dzρ ( z ) J h is the vertically integrated atmo-spheric enthalpy flux introduced in Eq. (55), F T OAR = F T OA,SWR − F T OA,LWR , where SW and LW refer tothe short- and long-wave contributions, respectively, and T E = (cid:16) F T OA,LWR /σ (cid:17) / is the emission temperature at agiven location. The contribution to the material entropyproduction coming from vertical processes can instead bewritten as: ˙ S vmat = (cid:90) A dA (cid:16) F surfR (cid:17) (cid:18) T s − T E (cid:19) (70)where F surfR = F surf,SWR + F surf,LWR is the net radi-ation at surface (defined as positive when the there isa net incoming radiation into the atmosphere) , SW and LW refer to the short- and long wave components,and, T s is the surface skin temperature defined as T s = (cid:16) F surf,LWR /σ (cid:17) / ∼ T surf ( Lucarini et al. , 2011). Equa-tions (69)-(70) allow one to compute the material en-tropy production due to internal irreversible processesmaking use only of 2D radiative fields at the boundariesof the relevant planetary fluid envelope (surface and topof the atmosphere). This makes Eqs. (69)-(70) suitablefor the post-processing of data hosted in publicly avail-able archives of GCMs output, intercomparison studies,1and studies of observational datasets of the Earth andother planets (where radiative data are the only avail-able source of information). Instead, direct computa-tions of ˙ S mat require the knowledge of the full 3-D time-dependent heating and temperature fields, making theirapplicability nontrivial for numerical models and unfea-sible for observations.Figure 10 shows a scatter-plot of the globally averagedannual mean values of the vertical and horizontal compo-nents of the material entropy production computed fromthe outputs of several GCMs from the CMIP3 datasetin pre-industrial and post-industrial conditions (updatedfrom ( Lucarini et al. , 2011), limiting to the models forwhich the data availability made possible the compari-son). The post-industrial case corresponds to the first100 years after the stabilization of the CO in the A1Bclimate change scenario. Issues related to the effectivenon-stationarity of the system have been treated as in( Lucarini and Ragone , 2011).Comparing with the direct computation ( Pascale et al. ,2011) of ˙ S mat for the case of Had-CM3 (model 13) in thepre-industrial case shows that the relative error on theestimate is less than 5% ( Lucarini et al. , 2011). The typ-ical values of the annual material entropy production inpre-industrial conditions range for most models between47 and 53 mW m − K − , matching well the approximateestimate by ( Ambaum , 2010). The contribution due tovertical processes is dominant by about one order of mag-nitude with respect to the contribution due to horizon-tal processes. This suggests that from the point of viewof the entropy production, the climate system approxi-mately behaves as a collection of weakly coupled verticalcolumns where mixing takes place ( Lucarini and Pascale ,2014).In increased CO concentration conditions, the rate ofmaterial entropy production increases for all the modelsbetween 10% and 20%. Such a change is dominated bythe increase in the vertical component, while the hori-zontal component sees in most cases a reduction of upto 10%, despite the fact that large scale horizontal en-thalpy transports increase for all models ( Lucarini andRagone , 2011). This implies - see Eq. (69) - a projectedstrong reduction in the large scale gradients of emissiontemperature, thus suggesting that in warmer conditionsthe climate system becomes more homogeneous in termsof meridional and zonal temperature differences. Thisfits well with what reported in ( Lucarini et al. , 2010a,b)in terms of climate response to global warming-like con-ditions, and hints to a dominant role of the latent heatrelease due to convective processes in the response to theclimate change ( Lucarini and Ragone , 2011).Figure 11a shows the spatial distribution of the in-tegrand of Eq. (70) for the Had-CM3 model in pre-industrial conditions ( Lucarini et al. , 2011), i.e. ,the localcontribution to the vertical component of material en- a) ˙ S hormat ( mW m − K − ) ˙ S v e r m a t ( m W m − K − ) b) −1 −0.5 0 0.5 1−10−8−6−4−20246810 47 89 13 1415 16 18 ∆ ˙ S hormat ( mW m − K − ) ∆ ˙ S v e r m a t ( m W m − K − ) FIG. 10 a) Scatter plot of contributions to the rate of materialentropy production due to horizontal (x-axis) and vertical (y-axis) processes. Each point corresponds to a GCM from theCMIP3 dataset in pre-industrial (black) and post-industrial(red) scenarios (updated from ( Lucarini et al. , 2011)). b)Difference between the SRESA1B scenario run (average ofthe last 30 years of the XXIII century and the pre-industrialclimatology. tropy production. Overall, the local material entropyproduction due to vertical processes seems to be a goodindicator of the geographical distribution of convectiveactivity: the highest values are observed in the warmpool of the western Pacific and Indian Ocean and in landareas characterized by warm and moist climates, whilerelatively low values are instead observed in the coldtongue of the eastern Pacific, near western boundary cur-rents, and in the temperate and cold oceans, as well ason deserts and middle and high latitudes of terrestrial2 a) o E o E o W o W o W o o S o S o o N o N b) −0.03−0.0275−0.025−0.0225−0.02−0.0175−0.015−0.0125−0.01−0.0075−0.005−0.0025 0 0.0025 0.005 0.0075 0.01 0.0125 0.015 0.0175 0.02 0.0225 0.025 0.0275 0.03 o E o E o W o W o W o o S o S o o N o N FIG. 11 (a) Spatial distribution of the contribution to therate of material entropy production due to vertical processesin pre-industrial scenario for Had-CM3 (model 13 in figure10). (b) Anomalies in the post-industrial scenario with re-spect to the pre-industrial case for the same model (updatedfrom ( Lucarini et al. , 2011)). areas. Note that also in this case the role of latent heatreleases is fundamental in determining the characteris-tics of the system, showing how the hydrological cycle isa crucial component of a thermodynamically consistentrepresentation of the climate system.Figure 11b shows the difference between the post-industrial and pre-industrial cases. The local verticalcomponent of material entropy production increases al-most everywhere, with negative anomalies confined topolar regions and to limited areas of the Southern Hemi-sphere, with very small values. The positive anomaliesare extremely high in the tropical regions, particularlyin the eastern and western Pacific ocean. Note that thepattern of increase does not strictly follow the patternof the absolute value in the pre-industrial case. In par-ticular the maximum of the increase is located eastwardto the maximum of the entropy production in the pre-industrial case, a signature of a shifting of the warm pooland a modification of the Walker circulation ( Bayr et al. ,2014; DiNezio et al. , 2013; IPCC , 2013). High values arealso found in the Indian ocean, suggesting an increaseof the convective activity connected with the Monsoon( IPCC , 2013; Turner and Annamalai , 2012). Significantlocal maxima are also observed in the Gulf of Mexico andalong the Gulf Stream, and in the Mediterranean Sea. The local entropy production due to vertical processesbehaves as a robust indicator of the impact of the cli-mate change on large-scale features connected to convec-tive activity. The pattern of increase is correlated to thepattern of variation of the surface temperature only to aminor extent. The reason is that this indicator containsin a synthetic way the information of the change in thesurface temperature, in the vertical stability of the atmo-sphere, and in the intensity of the energy fluxes connectedto the vertical processes. Therefore, it could be used inorder to define robust indexes for large-scale processesfor which strong convection is an important component.Moreover, the range of variation due to climate changeof the local vertical entropy production is rather high ifcompared to the range of variation of standard fields likesurface temperature or pressure. Therefore, one couldexpect a better signal-to-noise ratio and a more distinc-tive signature of climate change from indicators basedon this quantity compared with what obtained with in-dicators based on simpler observables, similarly to whatdiscussed by Lucarini et al. (2010a,b) and by Boschi et al. (2013) in the context of the identification of multi-stableregimes of the climate system. 2. Bistabiliy and tipping points Based on the evidence supported by Hoffman andSchrag (2002) and from numerical models ( Budyko , 1969; Ghil , 1976; Sellers , 1969), it is expected that the Earth ispotentially capable of supporting multiple steady statesfor the same values of some parameters such as, for ex-ample, the solar constant. Such states are the presentlyobserved warm state (W), and the entirely ice covered Snowball Earth state (SB). This is due to the presenceof two disjoint strange (chaotic) attractors. The W → SBand SB → W transitions are due, mathematically, to thecatastrophic disappearance of one of the two strange at-tractors ( Arnold , 1992) and, physically, to the positiveice-albedo feedback. The SB condition, which might bea common feature also of Earth-like planets, hardly al-lows for the presence of life, so this issue is of extremerelevance for defining habitability condition in extrater-restrial planets.PLASIM ( Fraedrich et al. , 2005), a general circulationmodel of intermediate complexity, was used by Boschiet al. (2013) and by Lucarini et al. (2013) to reconstructan extensive portion of the region of multistability inthe plane described by the parameters ( S ∗ , [CO ]). Thesurface temperature T s ( S ∗ , [CO ]) is shown in Fig. 12.The boundary of the domain in the parametric spacewhere two states are admissible correspond to the tippingpoints of the system.The thermodynamical and dynamical properties of theW and SB states are largely different. In the W states,surface temperature are 40 − 60 K higher than in the3 FIG. 12 Material entropy production (mW m − K − ) as afunction of solar constant S ∗ and the CO concentration. Thetransition SB → W and W → SB are marked with dashed arrowsstarting from the tipping point regions (courtesy of RobertBoschi, Universit¨at Hamburg). corresponding SB state and the hydrological cycle dom-inates the dynamics. This leads to a material entropyproduction (Fig. IV.C.2) larger by a factor of 4 – orderof (40 − × − W m − K − vs. (10 − × − Wm − K − – with respect to the corresponding SB states( Boschi et al. , 2013). The SB state is eminently a dryclimate, with entropy production mostly due to sensibleheat fluxes and dissipation of kinetic energy.The response to increasing temperatures of the twostates is rather different: the W states feature a decreaseof the efficiency of the climate machine, as enhanced la-tent heat transports reduces energy availability by damp-ening temperature gradients, while in the SB states theefficiency is increased, because warmer states are asso-ciated to lower static stability, which favors large scaleatmospheric motions (Fig. IV.C.2). The entropy produc-tion increases for both states, but for different reasons:the system become more irreversible and less efficient inthe case of W states, while stronger atmospheric motionslead to stronger dissipation and stronger energy trans-ports in the case of SB states. A general property whichhas been found is that, in both regimes, the efficiency η increases for steady states getting closer to tippingpoints and dramatically drops at the transition to thenew state belonging to the other attractor (Fig. IV.C.2).In a rather general thermodynamical context, this canbe framed as follows: the efficiency gives a measure ofhow far from equilibrium the system is. The negativefeedbacks tend to counteract the differential heating dueto the stellar insolation pattern, thus leading the systemcloser to equilibrium. At the bifurcation point, the nega-tive feedbacks are overcome by the positive feedbacks, sothat the system makes a global transition to a new state,where, in turn, the negative feedbacks are more efficient in stabilizing the system ( Boschi et al. , 2013). a)b)FIG. 13 a): Rate of material entropy production (10 − W m − K − ) vs. emission temperature T E ( K ) for Ω = Ω earth (magenta) and Ω = 0 . earth (black). b): as in a) figure butfor efficiency (courtesy of Robert Boschi, Universit¨at Ham-burg). Another interesting aspect is the determination of em-pirical functional relations between the main thermody-namical quantities and globally averaged emission tem-perature T E = ( LW toa /σ ) / , as shown in Fig. IV.C.2.This would permit to express non-equilibrium thermody-namical properties of the system in terms of parameterswhich are more directly accessible through measurements( Lucarini et al. , 2013). 3. Applications to planetary sciences The discovery of hundreds of planets outside the solarsystem (exoplanets) ( Seager and Deming , 2010) is ex-tending the scope of planetary sciences towards the studyof the so-called exoclimates ( Heng , 2012a). A large num-ber of the exoplanets discovered so far are tidally lockedto their parental star, experiencing extreme stellar forc-ing on the dayside where temperature up to 2000 K can4be reached. Starlight energy, deposited within the at-mosphere at the planet’s dayside, is then transported byatmospheric circulation to the night side. Such a system,similarly to the Earth’s climate, works like a heat engine(Sect. IV.A.2, Sect. IV.B).The strength of the day-to-night enthalpy flux con-trols the ratio of outgoing longwave energy fluxes fromthe night and day side ξ = LW night /LW day , called ef-ficiency of heat redistribution in the astrophysical liter-ature. Observations through infrared light curves showthat the hotter the planet, the more inefficient is theatmospheres at redistributing stellar energy leading tolarger day-night temperature differences. Numerical sim-ulations ( Perna et al. , 2012) show that ξ varies between0 . Am-baum , 2010; Johnson , 2000; Perna et al. , 2012; Schubertand Mitchell , 2013) and understanding their differenceswould be useful to provide a link between energy conver-sion and energy transport in planetary atmospheres.A thorough understanding of dissipative processesis fundamental for dealing with planetary atmospheres( Goodman , 2009; Pascale et al. , 2013). Dissipative pro-cesses are poorly known on Solar System planets and onexoplanets. Let us make some examples. In hot Jupiterstemperatures may be very high ( ≥ Saha equation ) andthus fast-moving (in hot Jupiters winds ∼ − ) elec-tric charges. This induces an electric current towards theinterior of the planet, where energy is then converted intoheat by ohmic dissipation.. Another dissipative mecha-nism believed to be a common feature in planetary atmo-spheres is shock wave breaking ( Batygin and Stevenson ,2010; Heng , 2012b). Note that the indirect method (Eq.68) could, in principle, be applied in order to infer in-formation about the dissipative processes in the interiorof exoplanets ( Schubert and Mitchell , 2013), where ra-diative fluxes are the only piece of information we canaccess. V. CLIMATE RESPONSE AND PREDICTION In the previous section, we have investigated the cli-mate as a non-equilibrium physical system and have em-phasized the intimate relation between forcing, dissipa-tion, energy conversion, and irreversibility. The sameapproach can be brought to a more theoretical level bytaking the point of view of non-equilibrium statisticalmechanics.Non-equilibrium statistical mechanics provides thenatural setting for investigating the mathematical prop-erties of forced and dissipative chaotic systems, which live in a non-equilibrium steady state (NESS). In this state,typically, the phase space contracts, entropy is gener-ated, and the predictability horizon is finite. Deviationsfrom this behavior are possible, but extremely unlikely.Conceptually, non-equilibrium steady states are gener-ated when a system is put in contact with reservoirs atdifferent temperatures or chemical potentials, and onedisregards the transient behaviors responsible for the re-laxation processes ( Gallavotti , 2006). This fits well thedescription of the non-equilibrium properties of the cli-mate system given in section IV.The science behind non-equilibrium statistical me-chanical systems is still in its infancy, so that, as op-posed to the equilibrium case, we are not able to predictthe properties of a system given the parameters describ-ing its internal dynamics and the boundary conditions,except in special cases where the dynamics is trivial.It is then important to choose a suitable mathemat-ical setting for being able to state some useful generalresults and compare numerical experiments with theory.The mathematical paradigm we will consider is the one ofso-called Axiom A systems ( Eckmann and Ruelle , 1985; Ruelle , 1989), which, according to the Chaotic Hypothe-sis ( Gallavotti , 1996), can be considered as good effective models of chaotic systems with many degrees of freedom.In general, we can say that an (time-continuous) Ax-iom A system ( Eckmann and Ruelle , 1985; Ruelle , 1989)obeys an evolution equation of the form ˙ x = F ( x ), x ∈ R n , and possesses an invariant measure ρ ( dx ) sup-ported on its attractor, which is, roughly speaking, theset of points where the system is asymptotically attractedto.If forcing and dissipation are present, the attractor is strange , i.e. , it does not look locally at all like a smoothmanifold, so that we cannot write ρ ( dx ) = ρ ( x ) dx , where ρ ( x ) is the density. Instead, in the very intuitive lan-guage of Lorenz, it looks like the Cartesian product ofa smooth manifold and a fractal set. The smooth man-ifold corresponds to the unstable directions of the flow,which make the system chaotic, while the Cantor set cor-responds to the contracting directions, which result fromdissipation. The invariant measure ρ ( dx ) gives the weightto be used in phase space to compute the expectation ofany observable A , which agrees, thanks to ergodicity, tothe long-time average, so that (cid:104) A (cid:105) = ρ ( A ) = (cid:90) ρ ( dx ) A ( x ) = lim T →∞ T (cid:90) T dtA ( x ( t ))with probability 1 with respect to the choice of the initialconditions.The invariant measure of an Axiom A systemis of Sinai-Ruelle-Bowen (SRB) type ( Eckmann and Ru-elle , 1985; Ruelle , 1989; Young , 2002). This has manyconsequences, including the fact that the measure is sta-ble against weak stochastic forcing, see also the discussionin Lucarini (2012).5 Ruelle (1997, 1998a,b, 2009) recently proved that inthe case of an Axiom A system, its SRB measure, despitethe geometrical complexity of the attractor supporting it,has also an extremely fascinating degree of regularity. Infact, there is a smooth dependence of the SRB measureto small perturbations of the flow, and it is possible toderive corresponding explicit formulas. This approach isespecially useful for studying the impact of changes inthe internal parameters of a system or of small modu-lations to the external forcing, and various studies havehighlighted the practical relevance of Ruelle theory forstudying what we may call the sensitivity of the systemto small perturbations. We will here recapitulate somefeatures of the Ruelle response theory and argue that itis a potentially useful tool for studying various classes ofGFD problems, and, most notably for addressing rigor-ously and in an unified perspective climate change pre-diction, climate response, and climate sensitivity. A. Response formulas and Fluctuation-Dissipation Theorem Let us consider an Axiom A dynamical system whoseevolution equation can be written as ˙ x = F ( x ) andlet’s assume that it possesses an invariant SRB measure ρ (0) ( dx ). Ruelle (1997, 1998a,b, 2009) has shown thatif the system is weakly perturbed so that its evolutionequation can be written as:˙ x = F ( x ) + Ψ( x ) T ( t ) (71)where Ψ( x ) is a weak time-independent forcing and T ( t )is its time modulation, it is possible to write the mod-ification to the expectation value of a general smoothobservable A as a perturbative series: ρ ( A ) t = ∞ (cid:88) n =0 ρ ( n ) ( A ) t , (72)where ρ (0) ( A ) t = ρ (0) ( A ) is the expectation value of A according to the unpertubed invariant measure ρ relatedto the dynamics ˙ x = F ( x ), while ρ ( n ) ( A ) t with n ≥ n th order processes( Lucarini , 2008).Limiting our attention to the linear case we have: ρ (1) ( A ) t = (cid:90) + ∞−∞ d τ G (1) A ( τ ) T ( t − τ ) , (73)where the first order Green function can be expressed asfollows: G (1) A ( τ ) = (cid:90) ρ ( dx )Θ( τ )Ψ( x ) · ∇ A ( x ( τ )) , (74)where Θ is the usual Heaviside distribution (Θ( x ) = 1 if x > 0, Θ( x ) = 0 if x < Kubo , 1957).In systems possessing a smooth invariant measure(which, as discussed above, is not typically the casefor Axiom A systems), like when equilibrium conditionsapply or stochastic forcing is imposed, we can write ρ ( dx ) = ρ ( x ) dx , where ρ ( x ) is the so-called density .In this case, we can rewrite Eq. (74) as follows: ρ (1) ( A ) t = (cid:90) + ∞−∞ d τ Θ( τ ) × (cid:90) dxρ ( x ) B ( x ) A ( x ( τ )) T ( t − τ ) , (75)where B ( x ) = −∇ · (cid:0) ρ ( x )Ψ( x ) (cid:1) /ρ ( x ). In other terms,one can predict the response at any time horizon t fromthe knowledge of the lagged correlation between the cho-sen observable A and the observable B , which depends onthe invariant measure ρ and on the perturbation vectorfield Ψ. See Colangeli and Lucarini (2014) for a detaileddiscussion on the physical meaning of B . Equation (75)provides a very general form of the FDT ( Lacorata andVulpiani , 2007; Ruelle , 1998a), which extends the resultsby ( Kubo , 1957). Recently, the FTD for system possess-ing a smooth invariant measure result has been extendedto the nonlinear case ( Lucarini and Colangeli , 2012).The more common forms of the FDT can be obtainedby taking one or more of the following assumptions: : • the perturbation flow is the form Ψ( x ) = (cid:15) ˆ x i ; • the observable is of the form A ( x ) = x j .where x k is the k th component of the x vector and ˆ x k is the corresponding unit vector. In this case, Eq. (75)takes the form: ρ (1) ( x j ) t = − (cid:15) (cid:90) + ∞−∞ d τ Θ( τ ) × (cid:90) dxρ ( x ) ∂ i log[ ρ ( x )] x j ( τ ) T ( t − τ ) , (76)If one takes the additional simplifying assumption thatunperturbed invariant measure has a Gaussian form, sothat ρ ( x ) = 1 /Z exp( − ˜ β (cid:80) Nj =1 x j / β > Z is a normalizing factor, we obtain: ρ (1) ( x j ) t = (cid:15) ˜ β (cid:90) + ∞−∞ d τ Θ( τ ) (cid:90) dxρ ( x ) x i x j ( τ ) T ( t − τ )= (cid:15)β (cid:90) + ∞−∞ d τ Θ( τ ) C i,j ( τ ) T ( t − τ ) , (77)where C i,j is the lagged correlation between x i and x j inthe unperturbed state.Unfortunately, the link between linear response of thesystem to external perturbations and its internal fluctua-tions seems more elusive when the unperturbed state has6a singular invariant measure. Ruelle (2009) shows thatsince the unperturbed invariant ρ (0) ( dx ) is singular, theresponse of the system contains two contributions, suchthat the first may be expressed in terms of a correlationfunction evaluated with respect to the unperturbed dy-namics along the space tangent to the attractor (unstablemanifold) and is formally identical to what given in Eq.(75). This part of the response decays rapidly due todecay of correlations due to chaos. On the other hand,the second term, which has no equilibrium counterpart,depends on the dynamics along the stable manifold, and,hence, it may not be determined from the unperturbeddynamics and is also quite difficult to compute numer-ically. These properties suggest the basic fact, alreadysuggested heuristically by Lorenz (1979), that in the caseof non-equilibrium systems internal and forced fluctua-tions of the system are not equivalent, the former beingrestricted to the unstable manifold only.Despite such a serious mathematical difficulty, the ap-plication of FDT, even in extremely simplified, quasi-Gaussian, approximation, has enjoyed a good success inclimate ( Gritsun and Branstator , 2007; Langen and Alex-eev , 2005) even if it is clear that the ability of FDT inpredicting the response to perturbation depends criticallyon the choice of the observable of interest, on the lengthof the integrations needed for constructing the approxi-mation of the invariant measure, and, of course, on thevalidity of the linear approximation ( Cooper and Haynes ,2011; Cooper et al. , 2013).There are, in fact, various ways to circumvent the prob-lem of the rigorous non-equivalence between forced andfree fluctuations. Apart from the obvious smoothing ef-fect due to unavoidable physical or numerical noise, whenconsidering smooth, coarse-grained observables(like thisof climatic interest), one expects to see little influence ofthe fine structure of the invariant measure of chaotic de-terministic systems, as projections from high-dimensionalspaces to lower dimensional ones are involved( Marconiet al. , 2008) and coarse-graining effects can be invoked( Wouters and Lucarini , 2013). One expects that theFDT will perform better in predicting the response of thesystem if one considers as observable A quantities like theglobally averaged surface temperature rather than, e.g. ,the surface temperature in an individual grid point. Fur-ther comments can be found at the end of section VI B. Computing the Response 1. Spectroscopic method If we select T ( t ) = (cid:15) cos( ω t ) = (cid:15)/ − iω t +exp( iω t )) as modulating factor of the perturbation field Ψ( x ), from equation Eq. (73) we derive:˜ ρ (1) ( A ) t = (cid:15)/ (cid:90) + ∞−∞ d τ G (1) A ( τ ) exp( − iω ( t − τ ))+ (cid:15)/ (cid:90) + ∞−∞ d τ G (1) A ( τ ) exp( iω ( t − τ ))= (cid:15)/ − iω t ) χ (1) A ( ω ) + c.c. (78)where χ (1) A ( ω ) is the Fourier Transform of G (1) A ( t ), usu-ally referred to as linear susceptibility, evaluated at fre-quency ω = ω , and c.c. indicates complex conjugate.Therefore, under the hypothesis of linearity, by perform-ing an ensemble of experiments where the forcing is of theform T ( t ) = (cid:15) cos( ω t ), we can extract the linear suscep-tibility at frequency ω by selecting the ω component ofthe Fourier transform of the signal ˜ ρ (1) ( A ) t obtained bytaking the ensemble average of the difference between thetime series of A in the perturbed and unperturbed case.By changing systematically the frequency ω of the forc-ing, one can reconstruct the susceptibility χ (1) A ( ω ) on achosen interval of frequencies. It is useful to recapitulatesome useful features of the susceptibility: • Resonances in the susceptibility function corre-spond to spectral ranges where the system is ex-tremely sensitive to forcings. In Fig. 14 we showthe real and imaginary part of the susceptibilityfor z variable of the ( Lorenz , 1963) model for the classical values of the parameters ( m = 1, σ = 10, r = 28, β = 8 / 3) and a given choice of the forc-ing (Ψ( x ) = [0; x ; 0] (cid:62) , T ( t ) = 2 (cid:15) cos( ωt )). We findthat for ω ∼ . 3, a very peaked spectral feature isapparent. Such a resonance is due to the UnstablePeriodic Orbits (UPO) of the system with the cor-responding period ( Eckhardt and Ott , 1994). UPOspopulate densely the attractors of chaotic systemsand constitute the so-called skeleton of the dynam-ics. In the case geophysical flows, UPOs have beenassociated to modes of low-frequency variability( Gritsun , 2008). One can, more qualitatively, asso-ciate resonance to positive feedbacks acting on timescales corresponding to the resonant frequency. • While | χ (1) A ( ω ) | measures the amplitude of the re-sponse of the system to perturbation at frequency ω , arctan( (cid:61){ χ (1) A ( ω ) } / (cid:60){ χ (1) A ( ω ) } ) gives the phasedelay between the forcing and the response, be-cause (cid:60){ χ (1) A ( ω ) } ( (cid:61){ χ (1) A ( ω ) } ) gives the compo-nent of the response that is in phase (out of phase)with the forcing. Depending on the forcing, on thesystem, and on the observable, this angle can varysignificantly even in a relatively small range of fre-quencies, as a result of resonances. • The two components (cid:61){ χ (1) A ( ω ) } and (cid:60){ χ (1) A ( ω ) } are connected by integral equations, the so-called7Kramers-Kronig relations ( Lucarini , 2008, 2009; Lucarini et al. , 2005; Ruelle , 2009; ? ). Such re-lations have their foundation in the causality of theGreen function (due to the presence of the Heavi-side distribution in Eq. 74) and establish a funda-mental connection between the response at differ-ent time scales: (cid:60){ χ (1) A ( ω ) } = 2 π P (cid:90) dω (cid:48) ω (cid:48) (cid:61){ χ (1) A ( ω (cid:48) ) } ω (cid:48) − ω ; (79) (cid:61){ χ (1) A ( ω ) } = − ωπ P (cid:90) dω (cid:48) ω (cid:48) (cid:60){ χ (1) A ( ω (cid:48) ) } ω (cid:48) − ω . (80)where P indicates that the integral is taken in prin-cipal part ( ? ). In particular one finds that: (cid:60){ χ (1) A (0) } = 2 π (cid:90) dω (cid:48) (cid:61){ χ (1) A ( ω (cid:48) ) } ω (cid:48) , (81)which provides a link between the static response- the sensitivity - and the out-of-phase responseat all frequencies. A large literature exists in op-tics, acoustics, condensed matter physics, particlephysics, signal processing on the theory and on themany applications of Kramers-Kronig relations andon the related sum rules , which provide integralconstraints related to the asymptotic behavior ofthe susceptibility ( Lucarini et al. , 2005).In Fig. 15 we present the real and imaginary part ofthe susceptibility of the mean energy e of the celebrated Lorenz (1996) model: dx i dt = x i − ( x i +1 − x i − ) − x i + F (82)where i = 1 , , ....., N , and the index i is cyclic so that x i + N = x i − N = x i , and e = 1 /N (cid:80) Nj =1 x j / 2. Thequadratic term in the equations simulates advection, thelinear one represents thermal or mechanical dampingand the constant one is an external forcing. See detailson the experiments in Lucarini and Sarno (2011), per-formed using N = 40 and F = 8. The system is per-turbed by the vector field Ψ( x ) = [1; . . . ; 1] (cid:62) modulatedby T ( t ) = 2 (cid:15) cos( ωt ). The resulting real and imaginarypart of χ (1) e ( ω ) are reported in Fig. 15, together withthe output of the data inversion performed via Kramers-Kronig relations. Once we obtain the susceptibility, asdiscussed in ( Lucarini and Sarno , 2011), it is possible toderive the corresponding Green function by applying theinverse Fourier Transform. This is the first applicationof the Kramers-Kronig theory in a geophysical context. 2. Broadband forcing If, instead, we select T ( t ) = δ ( t ), we derive from Eq.(74) that ρ (1) ( A ) t = G (1) A ( t ), i.e. ,the Green function cor- −15−10−505101520 ω R e [ χ ], I m [ χ ] FIG. 14 Measured real (blue line) and imaginary (red line)part of the susceptibility z variable of the Lorenz 63 model.Data from Lucarini (2009) . −2 −1 −0.500.511.52 ω R e [ χ ], I m [ χ ] ω l ω h FIG. 15 Measured real (blue line) and imaginary (red line)part of the susceptibility for the average energy of the Lorenz96 model. The rigorous extrapolation of the susceptibilityobtained via Kramers-Kronig analysis is reported (real part:black line; imaginary part: magenta line). Data from Lucariniand Sarno (2011) . responds to the relaxation of an ensemble of trajecto-ries of the system after a finite displacement along Ψ( x ).Obviously, we have that ˜ ρ (1) ( A ) ω = χ (1) A ( ω ), where the˜ symbol, indicates, as customary, that a Fourier Trans-form has been applied, so that the Fourier Transform ofthe signal is the linear susceptibility. Therefore, usingjust one ensemble of experiments where the perturbationis described by an impulsive forcing, we can gather thesame information on the response of the system which,in the previous case required an accurate sampling of dif-ferent frequencies.Let us look at the problem from a slightly more generalpoint of view. We apply the Fourier Transform to bothsides of Eq. (73) and obtain:˜ ρ (1) ( A ) ω = χ (1) A ( ω ) ˜ T ( ω ) (83)Choosing a sine or cosine function with argument ω t forthe function T ( t ) amounts to selecting as ˜ T ( ω ) the sumof two δ ’s centered in ω = ± ω . Therefore, the input8(forcing) allows only a small portion of the informationto derived on the system from the output (response). Letus assume that we choose the modulation T ( t ) such that˜ T ( ω ) is not vanishing for any ω , so that we have a broad-band modulation, where e.g. | ˜ T ( ω ) | for large values of ω decreases like a power law. If we perform an ensemble ofsimulations of the forced system, measure ˜ ρ (1) ( A ) ω , wecan invert Eq. 73 and readily derive: χ (1) A ( ω ) = ˜ ρ (1) ( A ) ω ˜ T ( ω ) (84)Therefore, one single set of experiments is, in fact all weneed to do to learn about the linear response propertiesof the system for the observable A . If we want to predictthe response at finite and infinite time of the system toforcing with the same spatial pattern Ψ( x ) but with dif-ferent time modulation R ( t ), we can derive G (1) A ( t ) from χ (1) A ( ω ) obtained via Eq. (84), and then plug it into Eq.(73). Alternatively, one can write:˜ ρ (1) ( A ) Rω = ˜ ρ (1) ( A ) Tω ˜ R ( ω )˜ T ( ω ) (85)where the upper indices R and T have been inserted forclarity, and then compute the inverse Fourier transformto derive the response at all times, or, if we apply theinverse Fourier transform to Eq. 84, we can compute theresponse to the R perturbation as: ρ (1) ( A ) Rt = (cid:90) + ∞−∞ d τ G (1) A ( τ ) R ( t − τ ) . (86) C. Prediction via Response theory The real test of the quality of an experimentally de-rived linear Green function G (1) A is the assessment ofits ability to support predictions about the system’s re-sponse to any temporal pattern of forcing R ( t ). The realbenefit of the broadband approach described here relieson exploiting linearity, and so deriving G (1) A from just oneensemble of simulations, each performed with the samemodulation T ( t ). Computing the G (1) A per se might be,in fact of little relevance.At this regard, we have performed additional experi-ments on the Lorenz (1996) model mirroring what pre-sented in section V.B.1. In this case, we have chosen astime modulation T ( t ) = (cid:15) Θ( t ), whose spectrum is indeedbroadband ( ˜ T ( ω ) /(cid:15) = πδ ( ω )+ i P[1 /ω ], where P indicatesthe principal part ) ( ? ). In this case, we have: G (1) e ( t ) = ddt ρ (1) ( e ) t . (87)Using about 1/100 of the computing time needed in Lu-carini and Sarno (2011), we have produced an estimate −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.500.511.522.53 t G ( ) e ( t ) FIG. 16 Linear Green function G (1) e ( t ) for the average energy e of the Lorenz (1996) model obtained by considering a step-like perturbation and using Eq. (87). Compare with Fig. 4in Lucarini and Sarno (2011). t ρ ( ) ( e ) t , R d τ G ( ) e ( t − τ ) T ( τ ) FIG. 17 Prediction of finite-time response of the average en-ergy e of the Lorenz (1996) model to a forcing with modula-tion T ( t ) = (cid:15) sin(2 πt ) ( (cid:15) = 0 . ρ (1) ( e ) t (blue line) vs. prediction obtained using the linear Greenfunction G (1) e ( t ) shown in Fig. 16. of the Green function of comparable quality; see Fig. 16.Additionally, we decided to check the predictive powerof the reconstructed Green function given in Fig. 16 bytesting its performance in predicting, through Eq. (73),the response of the system to a perturbation having tem-poral pattern given by T ( t ) = (cid:15) sin(2 πt ) ( (cid:15) = 0 . ρ (1) ( e ) t and the value predictedusing (cid:82) dτ G (1) e ( τ ) T ( t − τ ) is remarkable. One must em-phasize that the agreement is comparable if one selects (cid:15) = 1, thus moving away from the linear regime. D. Climate Response, Climate Change prediction Let us take inspiration from the previous example inorder to get some results of stricter geophysical rele-vance: we want to perform predictions on the impact9 a)b)FIG. 18 Studying climate change using response theory. a)Change in T S after an instantaneous doubling of the CO concentration. The lightly colored band indicates the two-standard deviation range around the ensemble mean. Insert:Green function of the T S . b) Comparison between GCM simu-lations (blue) and response theory prediction (red) for 1% peryear increase of the CO concentration up doubling. Lightlycolored band as in a). of increases in the CO concentration on the globallyaveraged surface temperature as simulated by a climatemodel, the simplified yet Earth-like PLASIM ( Fraedrichet al. , 2005). In what follows we present some new results(see also discussion in Ragone et al. (2014)), with the goalof proving the feasibility of the proposed methodology. • ˙ x = F ( x ) is the system of equations describing thediscretized version of a given model of the contin-uum PDEs describing the evolution of the climatein a baseline scenario with set boundary conditionsand values for, e.g. , the CO concentration and thevalue of the solar constant. We assume, for sim-plicity, that system model does not feature daily orseasonal variations in the radiative input at the top of the atmosphere. • Let us choose for the observable A the globally av-eraged surface temperature of the planet T S . • We study the perturbed system ˙ x = F ( x ) + f ( t )Ψ( x ). Let us choose as perturbation field Ψ( x )the convergence of radiative fluxes due to changein the logarithm of the atmospheric CO concen-tration. We want to be able to predict at finiteand infinite time the response of the the systemto one of the standard CO forcing scenario givenby the IPCC by performing an independent set ofperturbed model integrations. • The test perturbation is modulated by the function f ( t ) = (cid:15) Θ( t ), where (cid:15) is such that we double theamount of CO concentration in the atmosphere.Our goal is to predict the climate response to thecustomary 1% increase of CO concentration fromthe baseline value to its double. We select as base-line concentration [ CO ] = 360 ppm . • We perform 200 simulations, each lasting 200 yearsfor both scenarios of CO forcing. Our exper-iments are performed using PLASIM ( Fraedrichet al. , 2005) with a T21 spatial resolution, 10 ver-tical layers in the atmosphere, and swamp oceanhaving depth of 50 m.From the time series of the ensemble mean of the changeof T S - ρ (1) ( T S ) t - resulting from the sudden increase inthe CO , we derive the Green function G (1) T S ( t ) using Eq.(87). See Fig 18 a). Climate sensitivity is, in fact, definedby Eq. (81). Given the chosen pattern of forcing, we canrewrite is as follows:∆ T = (cid:60){ χ (1) T S (0) } = 2 π (cid:90) dω (cid:48) (cid:60){ ˜ ρ (1) ( T S ) ω (cid:48) } , (88)which relates climate response at all frequencies to itssensitivity.In order to test the predictive power of the responsetheory, we then convolute the Green function with thetemporal pattern of forcing of the second set of experi-ments.We choose as test experiment the classical IPCC sce-nario of 1% per year exponential increase of CO concen-tration up to doubling of the initial concentration (real-ized in τ ∼ 70 years, and constant concentration after-wards. Since our relevant control parameter is the log-arithm of the CO concentration, the second pattern offorcing g ( t ) is, in fact, a ramp increasing linearly withtime from 0 to (cid:15) in τ ∼ 70 years, with constant valueequal to (cid:15) for larger times. The results are presentedin Fig. 18 b), where we compare the predicted patternof increase (blue line), obtained using Eq. 86, with the0measured one (black line). The agreement is remarkable,both on the short and on the long time scales, while asome discrepancy exists between 20 and 50 years leadtime, where strong nonlinear effects due to ice-albedofeedback are dominant (not shown).Apparently, despite all the nonlinear feedbacks of theclimate model, the response to changes in the logarithmof CO concentration can be accurately described by lin-ear response theory at all time scales. Nonlinearity in theunderlying equations and presence of strong positive andnegative feedbacks do not rule out the possibility of con-structing accurate methods for predicting the response.In fact, the methods described here could be extended tothe nonlinear case by looking at the response in the fre-quency domain ( Lucarini , 2008, 2009), even if the dataquality requirement is obviously stricter.The result presented here suggests that many of thescenarios of greenhouse gases concentration included inthe IPCC reports ( IPCC , 2001, 2007, 2013) may in factbe partly redundant, as for certain variables might beaccurately described by linear response theory startingfrom just one scenario. Equations 84 - 85 constitute thebasis for predicting climate response at all scales.Obviously, with a given set of forced experiments, itis possible to derive the sensitivity to the the given forc-ing for as many climatic observables as desired. It isimportant to note that, for a given finite intensity (cid:15) ofthe forcing, the accuracy of the linear theory in describ-ing the full response depends also on the observable ofinterest. Moreover, the signal to noise ratio and, con-sequently, the time scales over which predictive skill isgood may change a lot from variable to variable. Theresults presented in this section extend to a more generalsetting and with stronger foundation the excellent intu-ition by Hasselmann et al. (1993) on the use of the linearresponse for addressing the problem of the so-called coldstart of coupled atmosphere-ocean models.Here we have shown results from just one observableprimary climatic interest. The analysis of other observ-ables will shed light on the mechanisms determining theclimate response to the forcing due to changes n the at-mospheric composition. As an example, the analysis ofthe response of large-scale meridional gradients of tem-perature at surface and in the middle troposphere willprovide information on changes in the midlatitude circu-lation. The existence of approximate functional relation-ship between the susceptibilities of different observables( Lucarini , 2009) would provide the key for defining rigor-ously the so-called emergent constraints ( Bracegirdle andStephenson , 2012).In practical terms, the applicability of response the-ory corresponds to having smooth dependence of climateproperties with respect to some given parameters. In-deed, this is not the case in the vicinity of tipping points(see Fig. 12). Response theory, may, nonetheless, suggestrigorous ways for defining and detecting tipping points, because one expects that these are associated to a diver-gence of the linear response.Finally, in order to talk about predictability, we needto specify what are the time scales over which we expectto have satisfactory predictive skills. In fact, linear re-sponse theory allows for deriving some scaling laws foraddressing this matter. The main obstacle for achievinga good degree of predictability is the uncertainty on theestimate of response signal given in Eq. (83) from theoutcomes of the numerical experiments because of thefiniteness of the ensemble and of the duration of each nu-merical simulation. See a detailed discussion of this issuein Ragone et al. (2014). VI. MULTISCALE SYSTEMS ANDPARAMETRIZATIONS The climate system features non-trivial behavior on alarge range of temporal and spatial scales ( Fraedrich andBttger , 1978; Lucarini , 2013; Peixoto and Oort , 1992; Vallis , 2006) When representing such a complex systemin a numerical simulation, the ratio of smallest to largesttime scale determines the number of required time steps,and the number of interactions between scales that haveto be calculated at each step can increase exponentiallywith the range of spatial variables. It is therefore clearthat, no matter which are the available computing re-sources, we are able to simulate explicitly only the vari-ables relevant for given ranges of spatial and temporalscales. Different choices of such ranges correspond to dif-ferent approximate theories of geophysical fluid dynamicsaimed at describing specific phenomenologies, a promi-nent case being that of quasi-geostrophic theory ( Klein ,2010).A manifestation of the inability to treat ultraslow vari-ability can be found in the usual practice in climate mod-eling of choosing fixed or externally driven boundary con-ditions, such as done when assuming a fixed extent for theland-based glaciers, and, consequently, for the sea-level,or imposing a specific path of CO concentration for theatmosphere. Instead, the impossibility of treating accu-rately fast processes requires the construction of so-calledparametrizations able to account, at least approximately,for the effect of the small scales on the large scales, asa function of the properties of the large scale variables,such in the case of several important physical processes,such as, e.g. , deep and shallow atmospheric convection,gravity wave drag, clouds, mixing in the ocean.Parametrizing small scale processes is important be-cause such unresolved processes impact the dynamics oflarger scales in terms of error growth, predictability, andclimatic biases. Presently, most of the parametrizationsused in climate models are deterministic, i.e. , for givenstate of the resolved variables, the effect of the unresolvedscale on the resolved scales is uniquely determined. We1often refer to these as bulk parametrizations . More re-cently, it has been emphasized that such a point of viewshould be modified for taking into account the fact thatthat many different states of the unresolved variables arecompatible with a given state of the resolved variables.This leads to considering the possibility of using stochas-tic parametrizations ( Palmer and Williams , 2009), whichshow promising abilities in reducing biases and reproduc-ing more effectively the uncertainties associated to per-forming mode reduction.When large time-scale separation exists between theresolved and non-resolved variables, the problem ofparametrization can be cast as follows. We consider asystem of the form ˙ Z = F ( Z ), Z ∈ R N and we dividethe state vector Z = ( X, Y ), where X are the slow com-ponents we are interested into and Y are the fast compo-nents we want to parametrize. We rewrite the evolutionequation as follows: dXdt = G X ( X, Y ) = F X ( X ) + Φ X ( X, Y ) dYdt = G Y ( X, Y ) = F Y ( Y ) + Φ Y ( X, Y ) (89)where we have split the dynamics of each set of variablesinto the autonomous part and into the coupling terms.The basic goal is to be able to write as an equation ofthe form dXdt = F X ( X ) + M X ( X ) + η ( X ) (90)where M X and η corresponds to the deterministic andstochastic components of the parametrization, respec-tively. A now classic example of empirical constructionand of testing of stochastic parametrizations is given by Wilks (2005), see Fig. 19.It must be emphasized that many of the approachesused so far have been based on the existence of a timescale separation between microscopic and macroscopicprocesses, following, conceptually, the pioneering pointof view proposed by Hasselmann (1976). If one does as-sume such a vast time-scale separation between the slowvariables X and the fast variables Y , averaging and ho-mogenization methods ( Arnold , 2001; Kifer , 2004; Pavli-otis and Stuart , 2008) allow for deriving an effective au-tonomous dynamics for the X variables, able to encom-pass the impact of the dynamics of the Y variables. Themotivation most often stated for the applicability of thistheory to climate science is the setting considered by Has-selmann, where fast weather systems influence slow cli-mate dynamics.Unfortunately, in many practical cases of interest ingeophysical fluid dynamics, such a scale separation doesnot exist - see, e.g. , the classical study by Mitchell (1976)- so that there is no spectral gap able to support univo-cally the identification of the X and Y variable. In fact,when the resolution of a numerical model is changed, all FIG. 19 Diagram describing how to parametrize the effect ofthe fast variables on the tendency of the slow variables X .The solid line U - see y − axis - corresponds to M X in Eq. 90,while the variability associated to the cloud of points needsto be represented via a stochastic term like η in Eq. 90. From Wilks (2005). the parametrizations have to be re-tuned, because the setof resolved variables has changed.Here we will focus on analytical methods that allowone to derive reduced models from the dynamical equa-tions of a full model. Projector operator techniques havebeen introduced in statistical mechanics with the goal ofeffectively removing the Y variables. In particular, con-siderable interest has been raised by the Mori-Zwanzigapproach, through which a formal - albeit practically in-accessible - solution for the evolution of the X variables isderived ( Mori , 1965; Zwanzig , 1961, 2001). These equa-tions in general contain both a correlated noise term anda memory term. Some attempts have been made to makeapproximation to the Mori-Zwanzig projected equationsto obtain practically useful equations. In applications ofstochastic mode reduction in climate science, the mem-ory term is usually not taken into account. This termcould however be very relevant in systems without a time-scale separation, as for example in the parametrizationof cloud formation in an atmospheric circulation model.The presence of memory in such systems has been dis-cussed by Bengtsson et al. (2013); Davies et al. (2009); Piriou et al. (2007). Note that, when we consider cou-pled systems where asynchronous is used, memory effectsare implicitly present in the dynamics.Besides considering the limit of infinite time scale sep-aration, another point of view can be taken, namely con-sidering the limit of weak coupling between the dynam-ical processes occurring at different scales. In this limit,the dynamics retains the correlated noise and memory2dependence that appeared in the Mori-Zwanzig reducedequations. The advantage of looking at this limit is how-ever that the noise autocorrelation function and memorykernel can now be written as simple correlation and re-sponse functions of the unresolved dynamics. A. Averaging and homogenization When applying averaging and homogenization tech-niques, one considers dynamical systems where a smallparameter (cid:15) controls the time scale separation between aslow and fast evolution in the system. The prototypicalset of equations for such a problem is dXdt = G X ( X, Y ) dYdt = G Y ( X, Y ) = 1 (cid:15) ˜ G Y ( X, Y )The parameter (cid:15) controls the time scale separation be-tween the variables X and Y , which becomes infinite as (cid:15) → Y , X will remain almostconstant. The fast variable Y will obey an evolution de-fined by ( X, Y ) for the current fixed value of X . On themuch longer time scale connected to the slow system, theevolution of X integrates out the rapid fluctuation of Y .As in the law of large numbers, the overall effect of allthese integrated fluctuations can be substituted by onesingle value. It can be shown that for finite time T , thefollowing applies: • the trajectory X ( t ) converges to a solution of: d ¯ Xdt = ¯ G ¯ X ( ¯ X )where ¯ G X ( ¯ X ) = ρ ¯ X ( G X ( ¯ X, Y )) is the averagedvalue of the tendency; • the average is taken over the invariant measure ρ X of the Y variable of the dynamical system dYdt = G Y ( ¯ X, Y )resulting when ¯ X is considered as a fixed forcingparameter.Let us consider a simple example system. dXdt =(1 − Y ) xdYdt = − (cid:15) Y + (cid:114) (cid:15) dWdt The Y system is here independent of X . The invariantmeasure of the fast y system is a Gaussian distribution with zero mean and unit variance. Taking the average of G X ( X, Y ) = (1 − Y ) x under the invariant measure of Y , we see that the averaged equation in this case is theuninteresting equation ˙ X = 0.This simple example immediately motivates the use ofhomogenization methods. Here one scales the equationto a longer time scale θ = (cid:15)t , the so called diffusive timescale and then performs the asymptotic expansion. Sim-ilarly to how correctly rescaling the sums of the law oflarge number leads to the more interesting central limittheorem, which describes the fluctuations around the av-erage value, also in the setting of time scale separatedsystems, we get stochastic behavior on the diffusive timescale. For the example considered above, we get a weakconvergence to a reduced stochastic differential equationfor the X variable instead of the trivial dynamical systemobtained before ( Pavliotis and Stuart , 2008).The theory for averaging and homogenization in time-scale separated stochastic differential equations is wellunderstood, with results for both one way and two waycouplings between the levels ( Bakhtin and Kifer , 2004).As usual, the theory is more complicated for determin-istic systems. Examples of dynamical systems can beconstructed where for a large set of initial conditions of Y , the solution for X does not converge to the averagedsolution ( Kifer , 2008). Furthermore, if the Y system haslong time correlations, such as in a system with regimebehavior, the homogenized system may converge badlyand an extension based on a truncation of the transferoperator has been proposed ( Sch¨utte et al. , 2004). Abramov (2012) has recently presented a study of un-certainty and predictability of the slow dynamics for asystem of geophysical relevance. A study of averagingand homogenization for idealized climate models, with arange of examples, can be found in Monahan and Culina (2011). Another rather successful attempt in this direc-tion is given in Majda et al. (2001). In Strounine et al. (2010) stochastic mode reduction is applied to a three-level quasi-geostrophic model whereas in Arnold et al. (2003) the authors perform mode reduction on a simplecoupled atmosphere-ocean model. Another applicationof homogenization to a toy model for the large-scale dy-namics of the atmosphere can be found in ( Frank andGottwald , 2013). Averaging for the case where one dealswith partial differential equations, as is relevant for cli-mate modeling, is discussed by Dymnikov (2012).A study of homogenization for geophysical flows wasperformed in Bouchet et al. (2013). The slow system isconsidered to be the evolution of zonal jets of a barotropicflow, which is forced by noise. The fast degrees of freedomare those representing the fast non-zonal turbulence. Ho-mogenization has also been applied in ( Dolaptchiev et al. ,2012) to the Burgers equation, where the slow variablesare taken to be averages over large grid boxes and thefast variables are the subgrid variables.When one wants to consider very large time scales (for3examples times of the order of exp(1 /(cid:15) )), one needs tolook beyond the central limit type theorems of homog-enization and consider so called large deviation results.These describe for example the transitions between dis-connected attractors of the averaged equations ( Kifer ,2009) and are of great relevance for studying tippingpoints ( Lenton et al. , 2008), going beyond simple one-dimensional approximate theories (see, e.g. , discussionin Lucarini et al. (2012)). B. Projection operator techniques Projection operator techniques do not constitute amode reduction per se, but are a way to rewrite the dy-namical equation of a multi-level equations to dependonly on a subset of variables. A projection is carried outon the level of the observables to remove unwanted, irrel-evant and usually fast degrees of freedom. The price onehas to pay for this apparent reduction is the appearanceof additional terms that are as difficult to compute as theoriginal system. It can however be a useful starting pointfor further approximations. These of techniques are alsoknown as the Mori-Zwanzig approach ( Mori et al. , 1974; Zwanzig , 1960, 1961).If a dynamical system is defined on a manifold M , onedefines a projection P from the space of observable func-tions on the full phase space M to a space of observableswhich are considered to contain only the interesting dy-namics. Many different choices are possible; if the mani-fold M consists for example of a product of submanifolds K of relevant and L of irrelevant variables, one can takea conditional expectation with respect to a measure on M , given the value of the relevant variables X ∈ K :( P A )( X ) = (cid:82) N A ( X, Y ) ρ ( X, Y ) dY (cid:82) N ρ ( X, Y ) dY . Another possible choice is a projection onto a set of func-tions on M , such as linear functions of the coordinatesin a Euclidean phase space. In general, one can think atvarious ways of performing coarse graining .Let us go back to our general formulation of a dy-namical system of the form ˙ Z = F ( Z ), Z ∈ R N .The evolution of an observable A ( Z ) can be written as˙ A ( Z ) = F ( Z ) ·∇ Z A ( Z ), which can be written as ˙ A = LA ,often referred to as Liouville equation. The evolution op-erator L is split into its projection P L onto the relevantspace of observables and the complement Q L := (1 −P ) L .As described by Zwanzig (2001), a generalized Langevinequation can then be derived based on Dyson’s formulafor operator exponentials e tL = e t Q L + (cid:90) t e ( t − s ) L P Le s Q L ds (91) We write the Liouville equation for an observable A as dA ( t ) dt = LA ( t ) = e tL LA = e tL P LA + e tL Q LA The factor exp( tL ) in the second term can be furtherexpanded by making use of Eq. (91). This gives thefollowing equation dA ( t ) dt = e tL P LA + ( e t Q L + (cid:90) t ds e ( t − s ) L P Le s Q L ) Q LA Zwanzig (2001) proposes the following interpretation ofthis equation. The first term on the right hand side cor-responds to the regular, deterministic dynamics of thesystem. The second term can be seen as describing acontribution from correlated noise, dependent on the ini-tial conditions of the irrelevant degrees of freedom. Thethird term (under integral) represents the memory of thesystem due to the presence of irrelevant variables thathave interacted with the relevant ones in the past. Inother term, the price we pay by separating somewhatarbitrarily relevant from irrelevant degrees of freedom sthat the irrelevant degrees of freedom act as a stochasticcomponent and, somewhat counter-intuitively, as proxiesfor the past state of the relevant degrees of freedom. Notethat we have done nothing more than manipulating theoriginal evolution equation ˙ A = LA . Correspondingly,the Mori-Zwanzig equation in itself does not simplify theproblem. In order to derive a set of equations that areuseful for numerical simulations, assumptions need to bemade about the dynamical system.Several approximations to the Mori-Zwanzig equationshave been proposed in the literature. There are the shortand long memory approximations made in the method ofoptimal prediction ( Bernstein , 2007; Chorin and Stinis ,2006; Chorin and Hald , 2013; Chorin et al. , 1998, 2000,2002, 2006; Defrasne , 2004; Hald and Kupferman , 2001; Park et al. , 2007).In the limit of an infinite time-scale separation be-tween the relevant and irrelevant variables, the stochas-tic component of the parametrization can be representedas a white noise term, while the memory (also knownas non-Markovian ) term vanishes, as the irrelevant vari-ables decorrelate quickly. Therefore, in such a limit theMori-Zwanzig decomposition is equivalent to the homog-enization method of section VI.A. For a comparison of theshort memory approximation of Mori-Zwanzig to homog-enization for climate-relevant models, see ( Stinis , 2006).We also refer the reader to recent results of Chekrounet al. (2013a,b), where general mathematical results forthe procedure of mode reduction, with thorough geomet-rical and dynamical interpretations, are given.Applications of the Mori-Zwanzig approach to fluid dy-namics can be found in ( Chandy and Frankel , 2009; Haldand Stinis , 2007; Hou , 2007; Stinis , 2007). A simpleapproximation to Mori-Zwanzig has been applied to jetformation on a β plane in ( Tobias and Marston , 2013).4 C. Weakly coupled systems We now consider dynamical systems consisting of twosystems with a weak coupling. In this case an expansionof the dynamics can be made in orders of the coupling,giving insight into what properties of the coupled systemsdetermines the memory kernel and correlated noise thatappeared in the Mori-Zwanzig approach ( Wouters andLucarini , 2012, 2013), because no assumptions are takenregarding time-scale separation.A possible application of this theory in climate sciencecan be found in the interaction between cloud forma-tion and large scale atmospheric flow, where there is nodistinct time scale separation, but instead the couplingcould be considered as weak. The weak coupling limit ofa tropical ocean-atmosphere model has also been consid-ered in the literature ( Neelin and Jin , 1993).Let us go back to Eq. 89. In this setting, the back-ground vector field F consists of a Cartesian product( F X , F Y ) (cid:62) of the vector fields F X and F Y defining theautonomous X and Y dynamics. The perturbing vec-tor field δF is a coupling (Ψ X , Ψ Y ) (cid:62) between the twosystems. We rewrite the full dynamical system as: dXdt = F X ( X ) + (cid:15) Ψ X ( X, Y ) dYdt = F Y ( Y ) + (cid:15) Ψ Y ( X, Y ) (92)where (cid:15) is added in order to clarify what kind of pertur-bative expansion we consider. For simplicity of presen-tation, for now we consider the case where Ψ X ( X, Y ) =Ψ X ( Y ) and Ψ Y ( X, Y ) = Ψ Y ( X ). We will come back tothe general case later.Given that the coupling term (cid:15) Ψ can be seen as a smallperturbation to the uncoupled system, one can make useof response theory to study the change of long time meansunder a change in the coupling parameter (cid:15) . We cantherefore use the response formalism described in SectionV. After lengthy calculations, one obtains the explicitexpression for ρ ( A ) t = ρ ( A ) t + ρ (1) ( A ) t + ρ (2) ( A ) t + O (Ψ ) (93) FIG. 20 Diagram describing the mean field effect of the Y variables on the X variables. Term M in Eq. (94). FIG. 21 Diagram describing the impact of fluctuations of the Y variables on the X variables. Term σ in Eq. (94).FIG. 22 Diagram describing the non-Markovian effect of theon the X variables on themselves, mediated by the Y vari-ables. Term h in Eq. (94). As shown by Wouters and Lucarini (2012), if one col-lects these first and second order responses to the cou-pling Ψ, an identical change in expectation values fromthe unperturbed ρ up to third order in Ψ can be ob-tained by adding a Y -independent forcing to the ten-dency of the X variables as follows: dX ( t ) dt = F X ( X ( t )) + M + σ ( t )+ (cid:90) ∞ dτ h ( τ, X ( t − τ )) (94)where M = ρ ,Y (Ψ X ) is an averaged version of the Y to X coupling, σ is a stochastic term, mimicking thetwo time correlation properties of the unresolved vari-ables and h is a memory kernel that introduces the non-Markovianity. A diagrammatic representation of pro-cesses responsible that these three additional terms areparametrizing is given in Figs. 20-22. Figure 20 refersto the mean field effect, which is captured by the firstorder correction, and corresponds to the deterministicparametrization. Figure 21 describes the effect of thefluctuations of the unresolved variables, which resultsinto an effective stochastic term in the parametrization.Finally, figure 22 clarifies how memory effects enter intothe picture of the parametrization: the resolved variablesat a given time impact the resolved variables at a latertime through a transfer of information mediated by theunresolved variables. The memory effect is present due tothe finite time scale difference between resolved and un-resolved variables, which also ensures that the stochastic5contribution shown in Fig. 21 cannot be represented by awhite noise process. In Wouters and Lucarini (2013) thisreduced equation was shown to be related to an expan-sion in the coupling strength of a Mori-Zwanzig equation.If the coupling functions Ψ X and Ψ Y are allowed tobe dependent on both X and Y , the above analysiscan still be carried out. In practical terms, this ac-counts for the possibility that the coupling terms arefunction of both the variables we want to parametrizeand of those we want to keep explicitly represented in ourmodel. For the case of separable couplings Ψ X ( X, Y ) =Ψ X, ( X )Ψ X, ( Y ) and Ψ Y ( X, Y ) = Ψ Y, ( X )Ψ Y, ( Y ) theaverage term becomes X dependent and the noise termbecomes multiplicative instead of additive. An expres-sion for more general couplings can be derived be de-composing the coupling functions into a basis of sepa-rable functions, see ( Wouters and Lucarini , 2012), andthen the same procedure can be applied. VII. SUMMARY AND CONCLUSIONS The goal of this review paper is the provision of anoverview of some ideas emerging at the interface betweenclimate science, physics, and mathematics, with the ob-jective of contributing to bridging the gap between dif-ferent scientific communities. The topics have been se-lected by the authors with the goal of covering (at leastpartially) relevant aspects of the deep symmetries of geo-physical flows, of the processes by which they convertand transport energy, and generate entropy, and of con-structing relevant statistical mechanical models able toaddress fundamental issues like the response of the cli-mate system to forcings, the representation of the in-teraction across scales, the definition of relevant physicalquantities able to describe succinctly the dynamics of thesystem. This review also informs the development andtesting of climate models of various degrees of complex-ity, by analyzing their physical and mathematical well-posedness and for constructing parametrizations of unre-solved processes, and by putting the basis for construct-ing diagnostic tools able to capture the most relevantclimate processes.The Nambu formulation of geophysical fluid dynam-ics explored in section II emphasizes the existence, inthe inviscid and unforced case, of non-trivial conservedquantities that are embedded in the equations of motion.Such quantities play a fundamental role, analogous to en-ergy’s, in the description of the state and of the dynamicsof the system and can be regarded as observables of greatrelevance also in the case where dissipation and forcingare present. Moreover, the Nambu formalism suggests usways for devising very accurate numerical schemes, whichdo not have spurious diffusive behavior.The symmetry properties of the flow in the inviscidlimit allow the construction of the ensembles describing the equilibrium statistical mechanical properties of thegeophysical flows (section III), where the vorticity - inthe two dimensional case - plays the role of the mostimportant physical quantity. Starting from the classi-cal construction due to Onsager of the gas of interact-ing vortices, the theory leads us to construct a theory ofbarotropic and baroclinic QG turbulence.Taking the point of view of non-equilibrium systems,we have that thanks to the presence of gradients of phys-ical quantities like temperature and chemical concentra-tions - in first instance due to the inhomogeneity of theincoming solar radiation, of the optical properties of thegeophysical fluids, and of the boundary conditions - theclimate system can transform available potential energyinto kinetic energy via internal instabilities, resulting inorganized fluid motions. In section IV the analysis ofthe energy and entropy budgets of the climate system isshown to provide a comprehensive picture of climate dy-namics, new tools for testing and auditing climate mod-els and measuring climate change, for investigating of theclimate tipping points, and for studying the properties ofgeneral planetary atmospheres.Section V introduces some basic concepts of non-equilibrium statistical mechanics, connecting the macro-scopic properties described in the previous section to thefeatures of the family of chaotic dynamical systems whichconstitute the backbone of the mathematical descriptionof non-equilibrium systems. For such systems, the re-lationship between internal fluctuations and response toforcings is studied with the goal of developing methodsfor predicting climate change. After clarifying the condi-tions under which the FDT is valid, we present some newresults such as a successful climate prediction for decadaland longer time scales. In this sense, we show that theproblem of climate change is mathematically well-posed.Non-equilibrium statistical mechanics is also the sub-ject of section VI, where we show how the Mori-Zwanzigformalism supports the provision of rigorous methods forconstructing parametrizations of unresolved processes. Itis possible to derive a surrogate dynamics for the coarsegrained variable of interest for climatic purposes, incor-porating, as result of the coupling with the small scale,fast variables, a deterministic, a stochastic, and a non-Markovian contribution, corresponding to memory ef-fects, which add to the unperturbed dynamics. The sameresults can be obtained using the response theory de-scribed in section V, thus showing that the constructionof parametrizations for weather and for climate modelsshould have common ground.Among the many topics and aspects left out of thisreview, we need to mention recent developments aimedat connecting the complementary, rather than oppos-ing ( Lorenz , 1963) and ( Hasselmann , 1976) perspectiveson complex dynamics dynamics, which focus on deter-ministic chaos and stochastic perturbations to dynami-cal systems, respectively. We refer in particular to the6idea of constructing time-dependent measures for nonautonomous dynamical systems ( Chekroun et al. , 2011)through the introduction of the so-called pullback attrac-tor , which is the geometrical object the trajectories ini-tialized in a distant past tend to at time t with proba-bility 1 as a result of the contracting dynamics. Such anobject is not invariant with time, as a result of the time-dependent forcing, but, under suitable conditions on theproperties of the dynamical system, the supported mea-sure has at each instant properties similar to those of the(invariant) SRB measure one can construct for, e.g. au-tonomous Axiom A dynamical ( Ruelle , 1989). Such anapproach allows for treating in a coherent way the pres-ence of modulations in the dynamics of the system, with-out the need of applying response formulas or of assumingtime-scale separations, and in particular allows for ana-lyzing the case where the forcing is stochastic, leading tothe concept of random attractor ( Arnold , 1988). On adifferent line of research, it is instead possible to use Ru-elle response theory for computing the impact of addingstochastic noise on chaotic dynamical systems ( Lucarini ,2012). One finds the rate of convergence of the stochas-tically perturbed measure to the unperturbed one, anddiscovers the general result that adding noise enhancesthe power spectrum of any given observables at all fre-quencies. The difference between the power spectrum ofthe perturbed and unperturbed system can be used, mir-roring an FDT, for computing the response of the systemto deterministic perturbations.The methods, the ideas, the perspectives presented inthis paper are partially overlapping, partially comple-mentary, partly in contrast. In particular, it is not obvi-ous, as of today, whether it is more efficient to approachthe problem of constructing a theory of climate dynam-ics starting from the framework of hamiltonian mechan-ics and quasi-equilibrium statistical mechanics or takingthe point of view of dissipative chaotic dynamical sys-tems, and of non-equilibrium statistical mechanics, andeven the authors of this review disagree. The former ap-proach can rely on much more powerful mathematicaltools, while the latter is more realistic and epistemolog-ically more correct, because, obviously, the climate is,indeed, a non-equilibrium system. Nonetheless, the ex-perience accumulated in many other scientific branches(chemistry, acoustics, material science, optics, etc.) hasshown that by suitably applying perturbation theory toequilibrium systems one can provide an extremely ac-curate description of non-equilibrium properties. Such alack of unified perspective, of well-established paradigms,should be seen as sign of the vitality of many researchperspectives in climate dynamics. ACKNOWLEDGMENTS The authors acknowledge interactions with U. Achatz,M. Ambaum, F. Cooper, M. Ghil, J. Gregory, K. Heng,T. Kuna, P. N´evir, O. Pauluis. D. Ruelle, A. Seifert, andR. Tailleux. The authors wish to thank F. Lunkeit and F.Sielmann for helping in data analysis, and R. Boschi, A.M¨uller and M. Sommer for providing useful figures. VL,RB, SP, and JW wish to acknowledge the financial sup-port provided by the ERC-Starting Investigator GrantNAMASTE (Grant no. 257106). The authors wish to ac-knowledge support by the Cluster of Excellence CliSAP.The National Center for Atmospheric Research is spon-sored by the National Science Foundation. Appendix A: Glossary For the benefit of the reader, we report here the mostrelevant symbols used in this paper, indicating if the samesymbol is used with different meaning. H Hamiltonian functional {• , •} P Standard Poisson Brackets u Velocity vector (two- or three-dimensional) u a Absolute velocity vector (including planetary rota-tion) ∇· Divergence operator (two- or three-dimensional) ∇ h · Horizontal divergence operator for three-dimensional vectors ∇ Gradient operator (two- or three-dimensional) ∇ h Horizontal gradient operator for three-dimensionalfields ψ streamfunction ω Vorticity function (two-dimensional dynamics) ω Vorticity vector (three-dimensional dynamics) ω a Absolute vorticity vector (including planetary vor-ticity) S Symplectic matrix [0, -1; 1, 0]; J Jacobian operator X Generic functional δ X /δa Functional derivative of X with respect to thefunction a . {• , • , •} Nambu brackets µ Horizontal divergence of the velocity field h T Total thickness of the fluid (shallow water equa-tions) h Helicity h a Absolute helicity (including planetary rotation) f = f + βy Planetary vorticity in β − plane approxi-mation ( y indicates the South-North coordinate)Φ Geopotential Q Quasi-geostrophic potential vorticity q Potential vorticity for shallow water equations (Sec-tion II); specific humidity (Section IV) N Brunt-V¨ais¨al¨a frequency7 ρ Density of the fluid (Sections II and IV); Invariantmeasure of the system (Sections III, V, and VI).Ω Earth’s angular velocity vectorΠ Ertel’s potential vorticity E Enstrophy functional E ρ Potential enstrophy functional S Entropy functional M Mass functional Z partition function R D = 2 π/k D Rossby deformation radius S Mixing entropy e Specific energy per unit mass i Specific internal energy per unit mass p pressure H relative humidity g Gravity T Temperature C v , C W Specific heat at constant volume p pressure s ( σ = ρs ) specific entropy per unit mass (per unitvolume)Ω( E ) structure function L Latent heat of vaporization (Section IV); Liouvilleoperator (Section VI) F R Vector of radiative flux. F S Vector of turbulent sensible heat flux. F L Vector of turbulent latent heat flux. τ Surface stress tensor E Total energy P Total static potential energy K Total kinetic energy W Conversion rate between potential and kinetic en-ergy D Rate of dissipation of the kinetic energy˙Φ + Net positive heating rate taking place at averagetemperature T + ˙Φ − Net negative heating rate taking place at averagetemperature T − η Climate efficiency˙ S mat Rate of material entropy production of the cli-mate system˙ S fric Rate of material entropy production due to fric-tion˙ S diff Rate of material entropy production due to dif-fusion˙ S hyd Rate of material entropy production due to thehydrological cycle˙ S vmat Rate of material entropy production of the cli-mate system due to vertical processes˙ S hmat Rate of material entropy production of the cli-mate system due to horizontal processes ρ (1) ( A ) First order correction to the expectation valueof the observable A G (1) A ( t ) First order Green function for the observableA χ (1) A ( ω ) First order susceptibility function for the ob-servable A∆ T Climate Sensitivity P = 1 − Q Projection operator performing coarsegraining on the dynamics and eliminates irrelevant de-grees of freedom; Q is the complementary operator REFERENCES Abramov, R. V. (2012), Suppression of chaos at slow variablesby rapidly mixing fast dynamics through linear energy-preserving coupling, Communications in Mathematical Sci-ences , (2), 595–624.Abramov, R. V., and A. Majda (2008), New approximationsand tests of linear fluctuation-response for chaotic nonlin-ear forced-dissipative dynamical systems, Journal of Non-linear Science , , 303–341, 10.1007/s00332-007-9011-9.Adams, D. K., and N. O. Renn`o (2005), Thermodynamic ef-ficiencies of an idealized global climate model, Clim. Dyn. , , 801–813.Ambaum, M. H. P. (2010), Thermal Physics of the Atmo-sphere , 256 pp., Wiley, New York.Arnold, L. (1988), Random Dynamical Systems , Springer,New York.Arnold, L. (2001), Hasselmann’s program revisited: Theanalysis of stochasticity in deterministic climate models, Stochastic climate models , , 141–158.Arnold, L., P. Imkeller, and Y. Wu (2003), Reduction of de-terministic coupled atmosphere-ocean models to stochas-tic ocean models: a numerical case study of the lorenz-maas system, Dynamical Systems , (4), 295–350, doi:10.1080/14689360310001607979.Arnold, V. (1992), Catastrophe theory, 3d edition , Springer,Berlin.Bakhtin, V., and Y. Kifer (2004), Diffusion approximation forslow motion in fully coupled averaging, Probability Theoryand Related Fields , (2), 157–181, doi:10.1007/s00440-003-0326-7.Bartello, P. (1995), Geostrophic adjustment and inverse cas-cades in rotating stratified turbulence., J. Atmos. Sci. , ,4410–4428.Basdevant, C., and R. Sadourny (1975), Ergodic proper-ties of inviscid truncated models of two-dimensional in-compressible flows, J. Fluid Mech. , , 673–688, doi:10.1017/S0022112075001620.Batygin, K., and D. J. Stevenson (2010), Inflating hot jupiterswith ohmic dissipation, The Astrophysical Journal Letters , (2), L238–L24.Bayr, T., D. Dommenget, T. Martin, and S. Power (2014),The eastward shift of the walker circulation in responseto global warming and its relationship to enso variability, Clim. Dyn. , doi:10.1007/s00382-014-2091-y.Becker, E. (2003), Frictional heating in global climate models, Mon. Weather Rev. , , 508–520.Bengtsson, L., M. Steinheimer, P. Bechtold, and J.-F. Ge-leyn (2013), A stochastic parametrization for deep con-vection using cellular automata, Quarterly Journal of theRoyal Meteorological Society , (675), 15331543, doi:10.1002/qj.2108.Bernstein, D. (2007), Optimal prediction of burgerss equa-tion, Multiscale Modeling & Simulation , (1), 27–52, doi:10.1137/060651720. Bihlo, A. (2008), Rayleigh-B´enard convection as a Nambu-metriplectic problem, J. Phys. A , , 292,001.Boffetta, G. (2007), Energy and enstrophy fluxes in the doublecascade of two-dimensional turbulence, J. Fluid Mech. , ,253–260.Boschi, R., V. Lucarini, and S. Pascale (2013),Bistability of the climate around the habitablezone: a thermodynamic investigation, Icarus , doi:http://dx.doi.org/10.1016/j.icarus.2013.03.017.Bouchet, F., and E. Simonnet (2009), Random changesof flow topology in two-dimensional and geophysi-cal turbulence, Phys. Rev. Lett. , , 94,504, doi:10.1103/PhysRevLett.102.094504.Bouchet, F., and A. Venaille (2012), Statistical mechanics oftwo-dimensional and geophysical flows, Phys. Rep. , ,227, doi:10.1016/j.physrep.2012.02.001.Bouchet, F., C. Nardini, and T. Tangarife (2013), KineticTheory of Jet Dynamics in the Stochastic Barotropic and2D Navier-Stokes Equations, Journal of Statistical Physics , , 572–625, doi:10.1007/s10955-013-0828-3.Bracegirdle, T. J., and D. B. Stephenson (2012), On the ro-bustness of emergent constraints used in multimodel cli-mate change projections of arctic warming, Journal of Cli-mate , (2), 669–678, doi:10.1175/JCLI-D-12-00537.1.Budyko, M. (1969), The effect of solar radiation variations onthe climate of the earth, Tellus , , 611–619.Chandy, A. J., and S. H. Frankel (2009), The t-model as alarge eddy simulation model for the navier-stokes equa-tions, Multiscale Modeling & Simulation , (2), 445–462,doi:10.1137/090760787, WOS:000277582900005.Charney, J. G. (1947), The dynamics of longwaves in a baroclinic westerly current, Journalof Meteorology , (5), 136–162, doi:10.1175/1520-0469(1947)004¡0136:TDOLWI¿2.0.CO;2.Chavanis, P.-H. (2009), Dynamical and thermodynamical sta-bility of two-dimensional flows: variational principles andrelaxation equations, Eur. Phys. J. B , , 73–105, doi:10.1140/epjb/e2009-00196-1.Chavanis, P.-H., and J. Sommeria (1996), Classification ofself-organized vortices in two-dimensional turbulence: thecase of a bounded domain, J. Fluid Mech. , , 267–297,doi:10.1017/S0022112096000316.Chavanis, P.-H., and J. Sommeria (1997), Thermody-namical approach for small-scale parametrization in2D turbulence, Phys. Rev. Lett. , , 3302–3305, doi:10.1103/PhysRevLett.78.3302.Chekroun, M. D., E. Simonnet, and M. Ghil (2011),Stochastic climate dynamics: Random attractorsand time-dependent invariant measures, Physica D:Nonlinear Phenomena , (21), 1685 – 1700, doi:http://dx.doi.org/10.1016/j.physd.2011.06.005.Chekroun, M. D., H. Liu, and S. Wang (2013a), On stochasticparameterizing manifolds: Pullback characterization andNon-Markovian reduced equations, ArXiv e-prints .Chekroun, M. D., H. Liu, and S. Wang (2013b), Non-Markovian Reduced Systems for Stochastic Partial Differ-ential Equations: The Additive Noise Case, ArXiv e-prints .Chorin, A., and P. Stinis (2006), Problem reduction, renor-malization, and memory, Communications in AppliedMathematics and Computational Science , (1), 1–27, doi:10.2140/camcos.2006.1.1.Chorin, A. J., and O. H. Hald (2013), Stochastic Tools inMathematics and Science , no. 58 in Texts in Applied Math-ematics, 3rd ed., Springer, New York. Chorin, A. J., A. P. Kast, and R. Kupferman (1998), Optimalprediction of underresolved dynamics, Proceedings of theNational Academy of Sciences , (8), 4094–4098, PMID:9539695.Chorin, A. J., O. H. Hald, and R. Kupferman (2000), Opti-mal prediction and the MoriZwanzig representation of ir-reversible processes, Proceedings of the National Academyof Sciences , (7), 2968–2973, doi:10.1073/pnas.97.7.2968,PMID: 10737778.Chorin, A. J., O. H. Hald, and R. Kupferman (2002), Optimalprediction with memory, Physica D: Nonlinear Phenomena , (34), 239–257, doi:10.1016/S0167-2789(02)00446-3.Chorin, A. J., O. H. Hald, and R. Kupferman (2006), Pre-diction from partial data, renormalization, and averaging, Journal of Scientific Computing , (2-3), 245–261, doi:10.1007/s10915-006-9089-5.Cichowlas, C., P. Bona¨ıti, F. Debbasch, and M. E. Brachet(2005), Effective Dissipation and Turbulence in SpectrallyTruncated Euler Flows, Phys. Rev. Lett. , (26), 264,502,doi:10.1103/PhysRevLett.95.264502.Colangeli, M., and V. Lucarini (2014), Elements of a unifiedframework for response formulae, Journal of Statistical Me-chanics: Theory and Experiment , (1), P01,002.Cooper, F. C., and P. H. Haynes (2011), Climate sensitivityvia a nonparametric fluctuation-dissipation theorem, Jour-nal of the Atmospheric Sciences , (5), 937–953.Cooper, F. C., J. G. Esler, and P. H. Haynes (2013), Es-timation of the local response to a forcing in a high di-mensional system using the fluctuation-dissipation theo-rem, Nonlinear Processes in Geophysics , (2), 239–248,doi:10.5194/npg-20-239-2013.Cullen, M. J. P. (2006), A Mathematical Theory of Large-ScaleAtmosphere/Ocean Flow , World Scientific Singapore.Cummins, P. F., and G. Holloway (1994), On eddy-topographic stress representation, J. Phys. Oceanogr. , (3), 700–706.Dauxois, T., S. Ruffo, E. Arimondo, and M. Wilkens (Eds.)(2002), Dynamics and Thermodynamics of Systems withLong Range Interactions , Lecture Notes in Physics , vol.602, Springer, New-York.Davies, L., R. S. Plant, and S. H. Derbyshire (2009), A sim-ple model of convection with memory, Journal of Geo-physical Research: Atmospheres , (D17), n/an/a, doi:10.1029/2008JD011653.Defrasne, S. (2004), G¨ute von Projektionsverfahren in der Dy-namik relevanter kollektiver Koordinaten in ausgew¨ahltenphysikalischen Modellen.deGroot, S., and P. Mazur (1984), Non-Equilibrium Thermo-dynamics , 510 pp., Dover.Dewar, R., C. Lineweaver, R. Niven, and K. Regenauer-Lieb(2013), Beyond the second law: an overview, in Beyondthe Second Law: Entropy Production and Non-equilibriumSystems , edited by D. RC, L. CH, N. RK, and R.-L. K,Springer.DiBattista, M., and A. J. Majda (2001), Equilibrium statisti-cal predictions for baroclinic vortices: The role of angularmomentum, Theor. Comput. Fluid Dyn. , (5), 293–322.Dijkstra, H. (2013), Nonlinear Climate Dynamics , CambridgeUniversity Press, Cambridge.DiNezio, P., G. Vecchi, and A. Clement (2013), Detectabilityof changes in the walker circulation in response to globalwarming, J. Climate , .Dolaptchiev, S., U. Achatz, and I. Timofeyev (2012),Stochastic closure for local averages in the finite-difference discretization of the forced burgers equation, Theoreti-cal and Computational Fluid Dynamics , pp. 1–21, doi:10.1007/s00162-012-0270-1.Donges, J. F., Y. Zou, N. Marwan, and J. Kurths (2009),Complex networks in climate dynamics, The EuropeanPhysical Journal Special Topics , (1), 157–179, doi:10.1140/epjst/e2009-01098-2.Donohoe, A., and D. Battisti (2012), What determines merid-ional heat transport in climate models?, J. Climate , ,3832–3859.Dubinkina, S., and J. Frank (2007), Statistical me-chanics of arakawas discretizations, Journal ofComputational Physics , (2), 1286 – 1305, doi:http://dx.doi.org/10.1016/j.jcp.2007.09.002.Dutton, J. A. (1973), The global thermodynamics of atmo-spheric motions, Tellus , , 89–110.Dvorak, R. (2008), Extrasolar Planets , Wiley.Dymnikov, V. P. (2012), Mathematics of climate modeling. ,Birkhauser, [S.l.].Eady, E. T. (1949), Long waves and cyclone waves, Tellus , (3), 33–52, doi:10.1111/j.2153-3490.1949.tb01265.x.Eckhardt, B., and G. Ott (1994), Periodic orbit analysis of thelorenz attractor, Zeitschrift fr Physik B Condensed Matter , (2), 259–266, doi:10.1007/BF01316970.Eckmann, J. P., and D. Ruelle (1985), Ergodic theory of chaosand strange attractors, Rev. Mod. Phys. , , 617–656, doi:10.1103/RevModPhys.57.617.Egger, J. (1999), Numerical generation of entropies, Mon.Wea. Rev. , , 2211–2216.Emanuel, K. (1991), The theory of hurricanes, Ann. Rev.Fluid. Mech. , , 179–196.Emanuel, K., and M. Bister (1996), Moist convective velocityand buoyancy scales, J. Atmos. Sci. , , 3276–3285.Enderton, D., and J. Marshall (2009), Controls on the totaldynamical heat transport of the atmosphere and oceans, J.Atmos. Sci. , , 1593–1611.Errico, R. M. (1984), The statistical equilibrium solution of aprimitive-equation model, Tellus A , (1), 42–51.Eyink, G. L., T. W. N. Haine, and D. J. Lea (2004), Ruelle’slinear response formula, ensemble adjoint schemes and lvyflights, Nonlinearity , (5), 1867.Fasullo, J., and K. Trenberth (2008), The annual cycle of theenergy budget. part ii: Meridional structures and polewardtransports, J. Climate , , 23132325.Fofonoff, N. P. (1954), Steady flow in a frictionless homoge-neous ocean, J. Mar. Res. , , 254–262.Fox, D. G., and S. A. Orszag (1973), Inviscid dynamics oftwo-dimensional turbulence, Phys. Fluids , , 169–171.Fraedrich, K., and H. Bttger (1978), A wavenumber-frequencyanalysis of the 500 mb geopotential at 50n, Journal of theAtmospheric Sciences , (4), 745–750, doi:10.1175/1520-0469(1978)035¡0745:AWFAOT¿2.0.CO;2.Fraedrich, K., and F. Lunkeit (2008), Diagnosing the entropybudget of a climate model, Tellus A , (5), 299–304.Fraedrich, K., H. Jansen, E. Kirk, U. Luksch, and F. Lunkeit(2005), The planet simulator: Towards a user friendlymodel, Meteorologische Zeitschrift , (3), 299–304, doi:doi:10.1127/0941-2948/2005/0043.Frank, J. E., and G. A. Gottwald (2013), Stochastic homog-enization for an energy conserving multi-scale toy modelof the atmosphere, Physica D: Nonlinear Phenomena , ,46–56, doi:10.1016/j.physd.2013.03.010.Frederiksen, J. S., and T. J. O’Kane (2008), Entropy, Closuresand Subgrid Modeling, Entropy , , 635–683. Frederiksen, J. S., and B. L. Sawford (1980), Statis-tical dynamics of two-dimensional inviscid flow on asphere., J. Atmos. Sci. , , 717–732, doi:10.1175/1520-0469(1980)037¡0717:SDOTDI¿2.0.CO;2.Gallavotti, G. (1996), Chaotic hypothesis: Onsager reci-procity and fluctuation-dissipation theorem, Journal ofStatistical Physics , (5-6), 899–925.Gallavotti, G. (2006), Stationary nonequilibrium statisticalmechanics, Encyclopedia of Mathematical Physics, ed. JPFrancoise, GL Naber, TS Tsun , , 530–539.Gassmann, A. (2013), A global hexagonal C-grid non-hydrostatic dynamical core (ICON-IAP) designed for ener-getic consistency , Quart. J. Roy. Meteorol. Soc. , (B),152–175, doi:10.1002/qj.1960.Gassmann, A. (2013), A global hexagonal C-grid nonhy-drostatic dynamical core (ICON-IAP) designed for ener-getic consistency, Quart. J. Roy. Meteorol. Soc. , , doi:doi:10.1002/qj.1960.Gassmann, A., and H.-J. Herzog (2008), Towards a consistentnumerical compressible non-hydrostatic model using gen-eralized Hamiltonian tools, Quart. J. Roy. Meteorol. Soc. , (635, B), 1597–1613, doi:10.1002/qj.297.Ghil, M. (1976), Climate stability for a Sellers-type model, Journal of Atmospheric Sciences , , 3.Ghil, M., P. Yiou, S. Hallegatte, B. D. Malamud, P. Naveau,A. Soloviev, P. Friederichs, V. Keilis-Borok, D. Kon-drashov, V. Kossobokov, O. Mestre, C. Nicolis, H. W.Rust, P. Shebalin, M. Vrac, A. Witt, and I. Zaliapin(2011), Extreme events: dynamics, statistics and predic-tion, Nonlinear Processes in Geophysics , (3), 295–350,doi:10.5194/npg-18-295-2011.Gianfelice, M., F. Maimone, V. Pelino, and S. Vaient (2012),On the recurrence and robust properties of lorenz?63model, Comm Math. Phys. , , 745–779.Goodman, J. (2009), Thermodynamics of atmospheric cir-culations on hot jupiters, The Astrophysical Journal , ,1645, doi:10.1088/0004-637X/693/2/1645.Goody, R. (2000), Sources and sinks of climate entropy, Q. J.R. Meteorol. Soc. , , 1953–1970.Goody, R. (2007), Maximum entropy production in climatetheory, J. Atmos. Sci. , , 2735–2739.Gritsun, A., and G. Branstator (2007), Climate response us-ing a three-dimensional operator based on the Fluctua-tionDissipation theorem, Journal of the Atmospheric Sci-ences , (7), 2558–2575.Gritsun, A. S. (2008), Unstable periodic trajectories of abarotropic model of the atmosphere, Russ. J. Numer. Anal.Math. Modelling , , 345–367.Hald, O. H., and R. Kupferman (2001), Convergence of opti-mal prediction for nonlinear hamiltonian systems, SIAMJournal on Numerical Analysis , (3), 983–1000, doi:10.1137/S0036142900374482.Hald, O. H., and P. Stinis (2007), Optimal predictionand the rate of decay for solutions of the euler equa-tions in two and three dimensions, Proceedings of theNational Academy of Sciences , (16), 6527–6532, doi:10.1073/pnas.0700084104.Hasselmann, K. (1976), Stochastic climate models, part I.theory, Tellus , (6), 473–485.Hasselmann, K., R. Sausen, E. Maier-Reimer, and R. Voss(1993), On the cold start problem in transient simulationswith coupled atmosphere-ocean models, Climate Dynam-ics , (2), 53–61, doi:10.1007/BF00210008.Hasson, S., V. Lucarini, and S. Pascale (2013), Hydrological cycle over south and southeast asian river basins as simu-lated by pcmdi/cmip3 experiments, Earth System Dynam-ics , (2), 199–217, doi:10.5194/esd-4-199-2013.Heng, K. (2012a), The study of climates of alien worlds, American Scientist , (4), 334–341.Heng, K. (2012b), On the existence of shocks in irradiatedexoplanetary atmospheres, The Astrophysical Journal Let-ters , (L1), 6.Herbert, C. (2013), Additional invariants and statistical equi-libria for the 2D Euler equations on a spherical domain, J.Stat. Phys. , , 1084–1114.Herbert, C. (2014), Nonlinear energy transfers andphase diagrams for geostrophically balanced rotating-stratified flows, Phys. Rev. E , , 033,008, doi:10.1103/PhysRevE.89.033008.Herbert, C., B. Dubrulle, P.-H. Chavanis, and D. Paillard(2012), Statistical mechanics of quasi-geostrophic flows ona rotating sphere, J. Stat. Mech. , , P05,023, doi:10.1088/1742-5468/2012/05/P05023.Herbert, C., A. Pouquet, and R. Marino (2014), RestrictedEquilibrium and the Energy Cascade in Rotating and Strat-ified Flows, submitted to J. Fluid Mech. , arXiv:1401.2103.Hernandez-Deckers, D., and J.-S. von Storch (2012), Impactof the warming pattern on global energetics, J. Climate , , 5223–5240.Hoffman, P., and D. Schrag (2002), The snowball earth hy-pothesis: testing the limits of global change, Terra Nova , , 129.Hogg, N. G., and H. M. Stommel (1985), Hetonic explosions:the breakup and spread of warm pools as explained by baro-clinic point vortices, J. Atmos. Sci. , , 1465–1476.Holloway, G. (1986), Eddies, waves, circulation, and mixing:Statistical geofluid mechanics, Ann. Rev. Fluid Mech. , ,91–147.Holloway, G. (1992), Representing Topographic Stress forLarge-Scale Ocean Models, J. Phys. Oceanogr. , , 1033–1046.Holloway, G. (2004), From classical to statistical ocean dy-namics, Surveys in Geophysics , , 203–219.Holton, J. (2004), An introduction to Dynamic Meteorology ,531 pp., Accademic Press.Hou, T. Y. (2007), Organized structures, memory,and the decay of turbulence, Proceedings of the Na-tional Academy of Sciences , (16), 6498–6499, doi:10.1073/pnas.0700639104.IPCC (2001), IPCC Third Assessment Report: WorkingGroup I Report ”The Physical Science Basis” , CambridgeUniversity Press.IPCC (2007), IPCC Fourth Assessment Report: WorkingGroup I Report ”The Physical Science Basis” , CambridgeUniversity Press.IPCC (2013), IPCC Fifth Assessment Report: WorkingGroup I Report ”The Physical Science Basis” , CambridgeUniversity Press.Johnson, D. (1989), The forcing and maintenance of globalmonsoonal circulations: an isentropic analysis, in Advancesin Geophysics , vol. 31, edited by B. Saltzman, pp. 43–316.Johnson, D. (1997), General coldness of climate models andthe second law: Implications for modeling the earth system, J. Climate , , 2826–2846.Johnson, D. R. (2000), Entropy, the Lorenz energy cycle, andclimate, in General Circulation Model Development: Past,Present and Future , edited by D. A. Randall, pp. 659–720,Academic Press, New York. Kalnay, E. (2003), Atmospheric Modeling, Data Assimila-tion and Predictability , Cambridge University Press, Cam-bridge.Kazantsev, E., J. Sommeria, and J. Verron (1998), Subgrid-scale eddy parameterization by statistical mechanics in abarotropic ocean model, J. Phys. Oceanogr. , , 1017–1042,doi:10.1175/1520-0485(1998)028¡1017:SSEPBS¿2.0.CO;2.Kifer, Y. (2004), Some recent advances in averaging, Moderndynamical systems and applications: dedicated to AnatoleKatok on his 60th birthday , p. 385.Kifer, Y. (2008), Convergence, nonconvergence and adiabatictransitions in fully coupled averaging, Nonlinearity , (3),T27, doi:10.1088/0951-7715/21/3/T01.Kifer, Y. (2009), Large deviations and adiabatic transitionsfor dynamical systems and markov processes in fully cou-pled averaging, Memoirs of the American Mathematical So-ciety , (944), 0–0, doi:10.1090/memo/0944.Kleidon, A. (2009), Non-equilibrium thermodynamics andmaximum entropy production, Naturwissenschaften , ,653–677.Klein, R. (2010), Scale-dependent models for atmosphericflows, Annual Review of Fluid Mechanics , (1), 249–274,doi:10.1146/annurev-fluid-121108-145537.Kraichnan, R. (1967), Inertial ranges in two-dimensional tur-bulence, Phys. Fluids , , 1417, doi:10.1063/1.1762301.Kraichnan, R., and D. Montgomery (1980), Two-dimensionalturbulence, Rep. Prog. Phys. , , 547, doi:10.1088/0034-4885/43/5/001.Kubo, R. (1957), Statistical-mechanical theory of irreversibleprocesses. i. general theory and simple applications to mag-netic and conduction problems, Journal of the Physical So-ciety of Japan , (6), 570–586, doi:10.1143/JPSJ.12.570.Kuroda, Y. (1991), On the Casimir invariant of Hamiltonianfluid mechanics, J. Phys. Soc. Japan , , 727–730.Lacorata, G., and A. Vulpiani (2007), Fluctuation-responserelation and modeling in systems with fast and slow dy-namics, Nonlin. Processes Geophys. , , 681–694.Landau, L. D., and E. M. Lifshits (1996), Mechanics ,Butterworth-Heinemann, Oxford.Langen, P. L., and V. A. Alexeev (2005), Estimating 2 xCO2 warming in an aquaplanet GCM using the fluctuation-dissipation theorem, Geophysical Research Letters , (23),L23,708.Ledwell, J. R., E. T. Montgomery, K. L. Polzin, L. C. St Lau-rent, R. W. Schmitt, and J. M. Toole (2000), Evidencefor enhanced mixing over rough topography in the abyssalocean, Nature , (6766), 179–182.Lee, T. D. (1952), On some statistical properties of hydro-dynamical and magneto-hydrodynamical fields, Q. Appl.Math. , , 69–74.Legg, S., and J. Marshall (1993), A heton model of thespreading phase of open-ocean deep convection, J. Phys.Oceanogr. , , 1040–1056.Lenton, T. M., H. Held, E. Kriegler, J. W. Hall, W. Lucht,S. Rahmstorf, and H. J. Schellnhuber (2008), Tipping el-ements in the earth’s climate system, Proceedings of theNational Academy of Sciences , (6), 1786–1793, doi:10.1073/pnas.0705414105.Li, L., A. Ingersoll, X. Jiang, D. Feldmann, and Y. Yung(2008), Lorenz energy cycle of the global atmosphere basedon reanalysis datasets, Geophys. Res. Lett. , , L16,813.Liepert, B., and F. Lo (2013), CMIP5 update of inter-modelvariability and biases of the global water cycle in CMIP3coupled climate models, Environ. Res. Lett. , , 029,401. Liepert, B., and M. Previdi (2012), Inter-model variabilityand biases of the global water cycle in CMIP3 coupled cli-mate models, Environ. Res. Lett. , , 014,006.Lorenz, E. (1955), Available potential energy and the main-tenance of the general circulation, Tellus , , 157–167.Lorenz, E. (1967), The nature and theory of the general cir-culation of the atmosphere, vol. 218.TP.115, World Mete-orological Organization.Lorenz, E. (1979), Forced and free variations of weather andclimate, J. Atmos. Sci , , 1367–1376.Lorenz, E. N. (1963), Deterministic nonperiodic flow, J. At-mos. Sci. , , 130–141.Lorenz, E. N. (1996), Predictability: A problem partly solved,in GARP Publication Series , vol. 16, pp. 132–136, WMO,Geneva, Switzerland.Lovejoy, S., and D. Schertzer (2013), The Weather and Cli-mate: Emergent Laws and Multifractal Cascades , Cam-bridge University Press, Cambridge.Lucarini, V. (2008), Response theory for equilibrium andnon-equilibrium statistical mechanics: Causality and gen-eralized kramers-kronig relations, Journal of StatisticalPhysics , , 543–558, 10.1007/s10955-008-9498-y.Lucarini, V. (2009), Evidence of dispersion relations for thenonlinear response of the Lorenz 63 system, Journal of Sta-tistical Physics , , 381–400, 10.1007/s10955-008-9675-z.Lucarini, V. (2012), Stochastic perturbations to dynamicalsystems: A response theory approach, Journal of StatisticalPhysics , (4), 774–786, doi:10.1007/s10955-012-0422-0.Lucarini, V. (2013), Modeling complexity: the case of climatescience, in Models, Simulations, and the Reduction of Com-plexity , edited by U. G¨ahde, S. Hartmann, and J. Wolf, pp.229–254, De Gruyter.Lucarini, V., and M. Colangeli (2012), Beyond the lin-ear fluctuation-dissipation theorem: the role of causality, Journal of Statistical Mechanics: Theory and Experiment , (05), P05,013.Lucarini, V., and K. Fraedrich (2009), Symmetry breaking,mixing, instability, and low-frequency variability in a min-imal Lorenz-like system, Phys. Rev. E , , 026,313.Lucarini, V., and S. Pascale (2014), Entropy production andcoarse graining of the climate fields in a general circu-lation model, Climate Dynamics , (3-4), 981–1000, doi:10.1007/s00382-014-2052-5.Lucarini, V., and F. Ragone (2011), Energetics of cli-mate models: Net energy balance and meridional en-thalpy transport, Rev. Geophys. , , RG1001, doi:10.1029/2009RG000323.Lucarini, V., and S. Sarno (2011), A statistical mechanicalapproach for the computation of the climatic response togeneral forcings, Nonlin. Proc. Geophys. , , 728.Lucarini, V., J. J. Saarinen, K.-E. Peiponen, and E. M. Var-tiainen (2005), Kramers-Kronig relations in Optical Mate-rials Research , Springer, New York.Lucarini, V., R. Danihlik, I. Kriegerova, and A. Speranza(2008), Hydrological cycle in the danube basin in present-day and xxii century simulations by ipccar4 global climatemodels, Journal of Geophysical Research: Atmospheres , (D9), n/a–n/a, doi:10.1029/2007JD009167.Lucarini, V., K. Fraedrich, and F. Lunkeit (2010a), Ther-modynamics of climate change: Generalized sensitivities, Atmos. Chem. Phys , , 9729–9737.Lucarini, V., K. Fraedrich, and F. Lunkeit (2010b), Thermo-dynamic analysis of snowball earth hysteresis experiment:Efficiency, entropy production, and irreversibility, QJRMS , , 2–11.Lucarini, V., K. Fraedrich, and F. Ragone (2011), New resultson the thermodynamic properties of the climate system, J.Atmos. Sci. , , 2438–2458.Lucarini, V., D. Faranda, and M. Willeit (2012), Bistablesystems with stochastic noise: virtues and limits of effectiveone-dimensional langevin equations, Nonlinear Processes inGeophysics , (1), 9–22, doi:10.5194/npg-19-9-2012.Lucarini, V., S. Pascale, R. Boschi, E. Kirk, and N. Iro (2013),Habitability and multistablility in earth-like plantets, Astr.Nach. , (6), 576–588.Majda, A., I. Timofeyev, and E. Vanden Eijnden (2001),A mathematical framework for stochastic climate models, Communications on Pure and Applied Mathematics , (8),891–974.Marconi, U. M. B., A. Puglisi, L. Rondoni, and A. Vulpiani(2008), Fluctuation-dissipation: Response theory in statis-tical physics, Phys. Rep. , , 111.Margules, M. (1905), On the energy of storms, transl. fromgerman by c. abbe, in Smithson. Misc. Collect. , edited byD. A. Randall, pp. 533–595.Marques, C., A. Rocha, and J. Corte-Real (2011), Global di-agnostic energetics of five state-of-the-art climate models, Clim. Dyn. , , 1767–1794.Mayer, M., and L. Haimberger (2012), Poleward atmosphericenergy transports and their variability as evaluated fromECMWF reanalysis data, J. Climate , , 734–752.Merryfield, W. J. (1998), Effects of stratification on quasi-geostrophic inviscid equilibria, J. Fluid Mech. , , 345–356.Miller, J. (1990), Statistical mechanics of Euler equationsin two dimensions, Phys. Rev. Lett. , , 2137–2140, doi:10.1103/PhysRevLett.65.2137.Mitchell, M. (1976), An overview of climatic variability andits causal mechanisms, Quaternary Research , (4), 481 –493, doi:http://dx.doi.org/10.1016/0033-5894(76)90021-1.Miyazaki, T., T. Sato, H. Kimura, and N. Takahashi(2011), Influence of external flow field on the equi-librium state of quasi-geostrophic point vortices, Geo-phys. Astrophys. Fluid Dyn. , (4-5), 392–408, doi:10.1080/03091929.2010.502118.Monahan, A. H., and J. Culina (2011), Stochastic av-eraging of idealized climate models, Journal of Cli-mate , (12), 3068–3088, doi:10.1175/2011JCLI3641.1,WOS:000291585800010.Montgomery, D., and G. Joyce (1974), Statistical mechanicsof “negative temperature” states, Phys. Fluids , , 1139,doi:10.1063/1.1694856.Mori, H. (1965), Transport, collective motion, and Brownianmotion, Progress of Theoretical Physics , (3), 423–455.Mori, H., H. Fujisaka, and H. Shigematsu (1974), A newexpansion of the master equation, Progress of TheoreticalPhysics , (1), 109–122, doi:10.1143/PTP.51.109.Nambu, Y. (1973), Generalized Hamiltonian dynamics, Phys.Rev. D , , 2403–2412.Naso, A., P.-H. Chavanis, and B. Dubrulle (2011), Statisticalmechanics of Fofonoff flows in an oceanic basin, Eur. Phys.J. B , , 493–517, doi:10.1140/epjb/e2011-10440-8.Neelin, J. D., and F.-F. Jin (1993), Modes of interannualtropical oceanatmosphere interactiona unified view. partII: Analytical results in the weak-coupling limit, Jour-nal of the Atmospheric Sciences , (21), 3504–3522, doi:10.1175/1520-0469(1993)050¡3504:MOITOI¿2.0.CO;2.N´evir, P. (1998), Die Nambu-Felddarstellungen der Hydro- Thermodynamik und ihre Bedeutung f¨ur die dynamischeMeteorologie, Habilitationsschrift.N´evir, P., and R. Blender (1993), A Nambu representationof incompressible hydrodynamics using helicity and enstro-phy, J. Phys. A , , L1189–1193.N´evir, P., and M. Sommer (2009), Energy-Vorticity Theory ofIdeal Fluid Mechanics, J. Atmos. Sci. , (7), 2073–2084,doi:10.1175/2008JAS2897.1.Nikurashin, M., G. K. Vallis, and A. Adcroft (2013),Routes to energy dissipation for geostrophic flows in theSouthern Ocean, Nature Geoscience , (1), 48–51, doi:10.1038/ngeo1657.Onsager, L. (1949), Statistical hydrodynamics, Il Nuovo Ci-mento , , 279–287, doi:10.1007/BF02780991.Oort, A., L. Anderson, and J. Peixoto (1994), Estimates ofthe energy cycle of the oceans, J. Geophys. Res. , , 7665–7688.Ozawa, H., A. Ohmura, R. Lorenz, and T. Pujol (2003), Thesecond law of thermodynamics and the global climate sys-tem: A review of the maximum entropy production princi-ple, Rev. Geophys , , 1–24.Palmer, T. N., and P. Williams (2009), Stochastic Physicsand Climate Modelling , Cambridge University Press, Cam-bridge.Paltridge, G. (1975), Global dynamics and climate – a systemof minimum entropy exchange, Quarterly Journal of theRoyal Meteorological Society , , 475484.Paltridge, G. (1978), The steady state format of global cli-mate, Quarterly Journal of Royal Meteorological Society , , 927–945.Paret, J., and P. Tabeling (1997), Experimental observationof the two-dimensional inverse energy cascade, Phys. Rev.Lett. , , 4162–4165.Park, J. H., N. S. Namachchivaya, and N. Neogi (2007),Stochastic averaging and optimal prediction, Journalof Vibration and Acoustics , (6), 803–807, doi:10.1115/1.2748777.Pascale, S., J. M. Gregory, M. Ambaum, and R. Tailleux(2011), Climate entropy budget of the hadcm3 atmosphere-ocean general circulation model and of FAMOUS, its low-resolution version, Clim. Dyn. , , 1189–1206.Pascale, S., F. Ragone, V. Lucarini, Y. Wang, and R. Boschi(2013), Nonequilibrium thermodynamics of an opticallythin, dry atmosphere, Planetary and Space Science , , 48–65.Pauluis, O. (2010), Water vapor and mechanical work:A comparison of carnot and steam cycles, Jour-nal of the Atmospheric Sciences , (1), 91–102, doi:10.1175/2010JAS3530.1.Pauluis, O., and I. Dias (2012), Satellite estimated ofprecipitation-induced dissipation in the atmosphere, Sci-ence , , 953–956.Pauluis, O., and I. M. Held (2002a), Entropy budget of an at-mosphere in radiative-convective equilibrium. Part I: Max-imum work and frictional dissipation, J. Atmos. Sci. , ,125–139.Pauluis, O., and I. M. Held (2002b), Entropy budget of an at-mosphere in radiative-convective equilibrium. Part II: La-tent heat transport and moist processes, J. Atmos. Sci. , , 140–149.Pavliotis, G., and A. Stuart (2008), Multiscale methods , Textsin applied mathematics : TAM, Springer, New York, NY.Pedlosky, J. (1987), Geophysical fluid dynamics , Springerstudy edition, Springer-Verlag, New York (N. Y.). Peixoto, J. P., and A. H. Oort (1992), Physics of climate ,American Institute of Physics, New York.Pelino, V., and F. Maimone (2007), Energetics, skele-tal dynamics, and long-term predictions on kolmogorov-lorenz systems, Phys. Rev. E , , 046,214, doi:10.1103/PhysRevE.76.046214.Perna, R., K. Heng, and F. Pont (2012), The effects of irra-diation on hot jovian atmospheres: heat redistribution andenergy dissipation, Astr. Journal , , 59–76.Piriou, J.-M., J.-L. Redelsperger, J.-F. Geleyn, J.-P.Lafore, and F. Guichard (2007), An approach for con-vective parameterization with memory: Separating mi-crophysics and transport in grid-scale equations, Jour-nal of the Atmospheric Sciences , (11), 4127–4139, doi:10.1175/2007JAS2144.1.Pouquet, A., and R. Marino (2013), Geophysical turbulenceand the duality of the energy flow across scales, Phys. Rev.Lett. , , 234,501, doi:10.1103/PhysRevLett.111.234501.Prigogine, I. (1961), Thermodynamics of Irreversible Pro-cesses , Interscience, New York.Probst, P., R. Rizzi, E. Tosi, V. Lucarini, and T. Maestri(2012), Total cloud cover from satellite observations andclimate models, Atmospheric Research , (0), 161 – 170,doi:http://dx.doi.org/10.1016/j.atmosres.2012.01.005.Qi, W., and J. B. Marston (2014), Hyperviscosity and sta-tistical equilibria of Euler turbulence on the torus and thesphere, to appear in J. Stat. Mech. Ragone, F., V. Lucarini, and F. Lunkeit (2014), A new frame-work for climate sensitivity and prediction, ArXiv e-prints .Rant, Z. (1956), Energie, ein neues wort f¨ur technische ar-beitf¨higkeit, Forsch. Ing. , , 36–37.Renn`o, N. O., and A. P. Ingersoll (1996), Natural convectionas a heat engine: a theory for CAPE, J. Atmos. Sci. , ,572–585.Rhines, P. B. (1975), Waves and turbulence on a beta-plane, J. Fluid Mech. , , 417–443.Rhines, P. B. (1976), The dynamics of unsteady currents, in The Sea , vol. VI, John Wiley & Sons, New-York, New York.Rhines, P. B. (1979), Geostrophic turbulence, Ann. Rev. FluidMech. , , 401–441.Robert, R., and J. Sommeria (1991), Statistical equilibriumstates for two-dimensional flows, J. Fluid Mech. , , 291–310, doi:10.1017/S0022112091003038.Rose, B., and D. Ferreira (2013), Ocean heat transport andwater vapor greenhouse in a warm equable climate: a newlook at the low gradient paradox, J. Climate , , 2127–2136.Ruelle, D. (1989), Chaotic Evolution and Strange Attractors ,Cambridge University Press, Cambridge.Ruelle, D. (1997), Differentiation of SRB states, Communica-tions in Mathematical Physics , (1), 227–241.Ruelle, D. (1998a), General linear response formula in sta-tistical mechanics, and the fluctuation-dissipation theoremfar from equilibrium, Phys. Lett. A , , 220–224.Ruelle, D. (1998b), Nonequilibrium statistical mechanics nearequilibrium: computing higher-order terms, Nonlinearity , (1), 5–18.Ruelle, D. (2009), A review of linear response theory for gen-eral differentiable dynamical systems, Nonlinearity , (4),855–870.Salazar, R., and M. V. Kurgansky (2010), Nambu brackets influid mechanics and magnetohydrodynamics, J. Phys. A , , 305,501(1–8). Salmon, R. (1978), Two-layer quasi-geostrophic turbulence ina simple special case, Geophys. Astrophys. Fluid Dyn. , ,25–52, doi:10.1080/03091927808242628.Salmon, R. (1988), Hamiltonian fluid mechanics, Annu. Rev.of Fluid Mech. , (1), 225–256.Salmon, R. (1998), Lectures on Geophysical Fluid Dynamics ,Oxford University Press.Salmon, R. (2005), A general method for conservingquantities related to potential vorticity in numericalmodels, Nonlinearity , (5), R1–R16, doi:10.1088/0951-7715/18/5/R01.Salmon, R. (2007), A general method for conserving energyand potential enstrophy in shallow water models, J. Atmos.Sci. , , 515–531.Salmon, R., G. Holloway, and M. C. Hendershott (1976),The equilibrium statistical mechanics of simple quasi-geostrophic models, J. Fluid Mech. , , 691–703, doi:10.1017/S0022112076000463.Saltzman, B. (2001), Dynamical Paleoclimatology , AcademicPress New York.Schneider, T. (2006), The general circulationof the atmosphere, Annual Review of Earthand Planetary Sciences , (1), 655–688, doi:10.1146/annurev.earth.34.031405.125144.Schubert, G., and J. Mitchell (2013), Planetary atmospheresas heat engines, in Comparative Climatology of TerrestrialPlanets , University of Arizona Press and Lunar and Plan-etary Institute.Sch¨utte, C., J. Walter, C. Hartmann, and W. Huisinga (2004),An averaging principle for fast degrees of freedom exhibit-ing long-term correlations, Multiscale Modeling & Simula-tion , (3), 501–526, doi:10.1137/030600308.Seager, S., and D. Deming (2010), Exoplanet atmospheres, Annual Review of Astronomy and Astrophysics , , 631–672.Sellers, W. (1969), A global climatic model based on the en-ergy balance of the earthatmosphere system, J. Appl. Me-teorol. , , 392–400.Shepherd, T. (1990), Symmetries, conservation laws, andhamiltonian structure in geophysical fluid dynamics, Adv.in Geophysics , , 287–338.Sommer, M., and P. N´evir (2009), A conservative schemefor the shallow-water system on a staggered geodesic gridbased on a Nambu representation, Quart. J. Roy. Meteorol.Soc. , (639), 485–494, doi:10.1002/qj.368.Steinheimer, M., M. Hantel, and P. Bechtold (2008), Con-vection in lorenz’s global energy cycle with the ECMWFmodel, Tellus A , , 1001–1022, doi:10.1111/j.1600-0870.2008.00348.x.Stinis, P. (2006), A comparative study of two stochasticmode reduction methods, Physica D: Nonlinear Phenom-ena , (2), 197–213, doi:10.1016/j.physd.2005.11.010.Stinis, P. (2007), Higher order mori-zwanzig models for theeuler equations, Multiscale Modeling & Simulation , (3),741–760, doi:10.1137/06066504X, WOS:000252254800002.Stone, P. H. (1978), Constraints on dynamical transports ofenergy on a spherical planet, Dyn. Atmos. Oceans , , 123–139.Storch, J.-S. V., C. Eden, I. Fast, H. Haak, D. Hernndez-Deckers, E. Maier-Reimer, J. Marotzke, and D. Stammer(2012), An estimate of the Lorenz energy cycle for theWorld Ocean Based on the 1 / ◦ STORM/NCEP simu-lation, J. Phys. Oceanogr. , , 2185–2205.Strounine, K., S. Kravtsov, D. Kondrashov, and M. Ghil (2010), Reduced models of atmospheric low-frequency vari-ability: Parameter estimation and comparative perfor-mance, Physica D: Nonlinear Phenomena , (3-4), 145–166, doi:10.1016/j.physd.2009.10.013.Tailleux, R. (2013), Available potential energy and exergy instratified fluids, Ann. Rev. Fluid. Mech , , 35–58.Taylor, K. E., R. J. Stouffer, and G. A. Meehl (2012), Anoverview of CMIP5 and the experiment design, Bull. Amer.Meteor. Soc. , , 485–498.Tobias, S. M., and J. B. Marston (2013), Directstatistical simulation of out-of-equilibrium jets, Physical Review Letters , (10), 104,502, doi:10.1103/PhysRevLett.110.104502.Trenberth, K., and J. Fasullo (2012), Tracking earths energy:From el nio to global warming, Surveys in Geophysics , (3-4), 413–426, doi:10.1007/s10712-011-9150-2.Trenberth, K. E., and J. M. Caron (2001), Estimates of merid-ional atmosphere and ocean heat transports, J. Clim , ,3433–3443.Trenberth, K. E., and J. T. Fasullo (2010), Simulation ofpresent-day and twentyfirst-century energy budgets of thesouthern oceans, J. Climate , , 440–454.Turner, A., and H. Annamalai (2012), Climate change andthe south asian summer monsoon, Nature Clim. Change , , doi:10.1038/nclimate1495.Vallis, G. K. (2006), Atmospheric and Oceanic Fluid Dynam-ics: Fundamentals and Large-scale Circulation , CambridgeUniversity Press, Cambridge.Vallis, G. K., and M. E. Maltrud (1993), Generation of meanflows and jets on a beta-plane and over topography, J. Phys.Oceanogr. , , 1346–1362.Venaille, A. (2012), Bottom-trapped currents as statisticalequilibrium states above topographic anomalies, J. FluidMech. , , 500, doi:10.1017/jfm.2012.146.Venaille, A., and F. Bouchet (2011), Solvable phase diagramsand ensemble inequivalence for two-dimensional and geo-physical turbulent flows, J. Stat. Phys. , , 346–380, doi:10.1007/s10955-011-0168-0.Venaille, A., G. K. Vallis, and S. M. Griffies (2012), The cat-alytic role of beta effect in barotropization processes, J.Fluid Mech. , , 490–515.Wang, J., and G. K. Vallis (1994), Emergence of Fofonoffstates in inviscid and viscous ocean circulation models, J.Mar. Res. , , 83–127.Warn, T. (1986), Statistical mechanical equilibria of the shal-low water equations, Tellus A , (1), 1–11.Wild, M., D. Folini, C. Schr, N. Loeb, E. Dutton, andG. Knig-Langlo (2013), The global energy balance froma surface perspective, Climate Dynamics , (11-12), 3107–3134, doi:10.1007/s00382-012-1569-8.Wilks, D. S. (2005), Effects of stochastic parametrizationsin the Lorenz ’96 system, Quarterly Journal of the RoyalMeteorological Society , (606), 389–407.Woollings, T., and J. Thuburn (2006), Entropy sources ina dynamical core atmosphere model, Quarterly Journalof the Royal Meteorological Society , (614), 43–59, doi:10.1256/qj.04.189.Wouters, J., and V. Lucarini (2012), Disentangling multi-levelsystems: averaging, correlations and memory, Journal ofStatistical Mechanics: Theory and Experiment , (03),P03,003.Wouters, J., and V. Lucarini (2013), Multi-level dynami-cal systems: Connecting the ruelle response theory andthe mori-zwanzig approach, Journal of Statistical Physics , (5), 850860, doi:10.1007/s10955-013-0726-8.Wunsch, C. (2005), The total meridional heat flux and itsoceanic and atmospheric partition, J. Clim. , , 4374–4380.Wunsch, C. (2012), Discrete Inverse and State Estimation -With Geophysical Fluid Applications , Cambridge Univer-sity Press, Cambridge.Wunsch, C., and R. Ferrari (2004), Vertical Mix-ing, Energy, and the General Circulation of theOceans,