Get 20M+ Full-Text Papers For Less Than $1.50/day. Subscribe now for You or Your Team.

Learn More →

A statistical mechanical approach for the computation of the climatic response to general forcings

A statistical mechanical approach for the computation of the climatic response to general forcings Nonlin. Processes Geophys., 18, 7–28, 2011 Nonlinear Processes www.nonlin-processes-geophys.net/18/7/2011/ doi:10.5194/npg-18-7-2011 in Geophysics © Author(s) 2011. CC Attribution 3.0 License. A statistical mechanical approach for the computation of the climatic response to general forcings 1,2,3 3 V. Lucarini and S. Sarno Department of Meteorology, University of Reading, Reading, UK Department of Mathematics, University of Reading, Reading, UK Walker Institute for Climate System Research, University of Reading, Reading, UK Received: 4 August 2010 – Revised: 18 November 2010 – Accepted: 4 December 2010 – Published: 12 January 2011 Abstract. The climate belongs to the class of non- equations for such parameters allow to define such properties equilibrium forced and dissipative systems, for which most as an explicit function of the unperturbed forcing parameter results of quasi-equilibrium statistical mechanics, including alone for a general class of chaotic Lorenz 96 models. We the fluctuation-dissipation theorem, do not apply. In this pa- then verify the theoretical predictions from the outputs of the per we show for the first time how the Ruelle linear response simulations up to a high degree of precision. The theory is theory, developed for studying rigorously the impact of per- used to explain differences in the response of local and global turbations on general observables of non-equilibrium statisti- observables, to define the intensive properties of the system, cal mechanical systems, can be applied with great success to which do not depend on the spatial resolution of the Lorenz analyze the climatic response to general forcings. The crucial 96 model, and to generalize the concept of climate sensitivity value of the Ruelle theory lies in the fact that it allows to com- to all time scales. We also show how to reconstruct the linear pute the response of the system in terms of expectation values Green function, which maps perturbations of general time of explicit and computable functions of the phase space aver- patterns into changes in the expectation value of the consid- aged over the invariant measure of the unperturbed state. We ered observable for finite as well as infinite time. Finally, choose as test bed a classical version of the Lorenz 96 model, we propose a simple yet general methodology to study gen- which, in spite of its simplicity, has a well-recognized pro- eral Climate Change problems on virtually any time scale totypical value as it is a spatially extended one-dimensional by resorting to only well selected simulations, and by tak- model and presents the basic ingredients, such as dissipa- ing full advantage of ensemble methods. The specific case of tion, advection and the presence of an external forcing, of globally averaged surface temperature response to a general the actual atmosphere. We recapitulate the main aspects of pattern of change of the CO concentration is discussed. We the general response theory and propose some new general believe that the proposed approach may constitute a mathe- results. We then analyze the frequency dependence of the re- matically rigorous and practically very effective way to ap- sponse of both local and global observables to perturbations proach the problem of climate sensitivity, climate prediction, having localized as well as global spatial patterns. We derive and climate change from a radically new perspective. analytically several properties of the corresponding suscep- tibilities, such as asymptotic behavior, validity of Kramers- Kronig relations, and sum rules, whose main ingredient is the 1 Introduction causality principle. We show that all the coefficients of the leading asymptotic expansions as well as the integral con- A crucial goal in the study of general dynamical and statisti- straints can be written as linear function of parameters that cal mechanical systems is to understand how their statistical describe the unperturbed properties of the system, such as properties are altered when we introduce a perturbation re- its average energy. Some newly obtained empirical closure lated to changes in the external forcing or in the value of some internal parameters. The ability to compute the re- sponse of the system is of great relevance for purely math- ematical reasons as well as in many fields of science and Correspondence to: V. Lucarini technology. (v.lucarini@reading.ac.uk) Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union. 8 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings The climate system is an outstanding example of a non- lent. Whereas natural fluctuations of the system are restricted equilibrium, forced and dissipative complex system, forced to the unstable manifold, because, by definition, asymptot- in first instance by spatial differences and temporal variabil- ically there is no dynamics along the stable manifold, ex- ity in the net energy flux at the top of the atmosphere. On a ternal forcings will induce motions – of exponentially de- macroscopic level, as a result of being far from equilibrium, caying amplitude – out of the attractor with probability one. the climate system behaves as an engine, driven by the tem- The fluctuation-dissipation relation can be recovered only if perature difference between a warm and a cold thermal pool, we consider perturbations with the somewhat artificial prop- so that the atmospheric and oceanic motions are at the same erty of being everywhere tangent to the unstable manifold or, time the result of the mechanical work (then dissipated in a in a more fundamental way, if we add a stochastic forcing, turbulent cascade) produced by the engine, and are processes which has the crucial effect of smoothing the invariant mea- which re-equilibrate the energy balance of the climate sys- sure (Lacorata and Vulpiani, 2007; Marini Bettolo Marconi tem (Lorenz, 1967; Peixoto and Oort, 1992; Johnson, 2000; et al., 2008). Potential links to these issues can be found in Lucarini, 2009a). recent papers proposing new algorithms for three (Trevisan A primary goal of climate science is to understand how the and Uboldi, 2004) and four (Trevisan et al., 2010) dimen- statistical properties – mean values, fluctuations, and higher sional variational data assimilation, where it is shown that order moments – of the climate system change as a result of the quality of the procedure improves if the increment of the modulations to some crucial external (e.g. solar irradiance) variables due to the assimilation is performed only along the or internal (e.g. atmospheric composition) parameters of the unstable manifold. system occurring on various time scales. A large class of Recently, Ruelle (1998a, 2009) introduced a rigorous problems – those involving climate sensitivity, climate vari- mathematical theory allowing for computing analytically, ability, climate change, climate tipping points – fall into this ab initio, the response of a large class of non-equilibrium category. In a system as complex and as extended as the cli- systems to general external perturbations featuring arbitrary mate, where lots of feedbacks are active on a variety of spa- time modulation. The crucial result is that the changes in tial and temporal scales, this is in general a very difficult task. the expectation value of a physical observable can be ex- The need for scientific advance in this direction is outstand- pressed as a perturbative series in increasing powers of the ing as one considers that even after several decades of intense intensity of the external perturbation, where each term of the scientific efforts, the accurate evaluation of the climate sensi- series can be written as the expectation value of some well- tivity par excellence, i.e., the change of the globally averaged defined observable over the unperturbed state. In a previ- surface temperature for doubling of CO concentration with ous paper (Lucarini, 2008b) we showed that the Ruelle the- respect to pre-industrial levels (280 ppm to 560 ppm circa), is ory is, thanks to this property, formally analogous to usual a tantalizing endeavor, and large uncertainties are still present Kubo response theory (Kubo, 1957), which applies for quasi- (IPCC, 2007). equilibrium system. The crucial difference lies on the math- Such efforts have significant relevance also in the context ematical properties of the invariant measure, which is abso- of the ever-increasing attention paid by the scientific commu- lutely continuous in the quasi-equilibrium case and singular nity to the quest for reliable metrics to be used for the valida- in the non-equilibrium case. tion of climate models of various degrees of complexity and Ruelle’s analysis applies for non-equilibrium steady state for the definition of strategies aimed at the radical improve- systems (Gallavotti, 2006) possessing a Sinai-Ruelle-Bowen ment of their performance (Held, 2005; Lucarini, 2008a). (SRB) invariant measure, often referred to as Axiom A sys- The pursuit of a quantum leap in climate modelling – which tem (Eckmann and Ruelle, 1985; Ruelle, 1989). This class definitely requires new scientific ideas rather than just faster of systems, even if mathematically non-generic, includes on supercomputers – is becoming more and more of a key issue the other hand excellent models for general physical systems, in the climate community (Shukla et al., 2009). as made clear by the chaotic hypothesis (Gallavotti and Co- A serious, fundamental difficulty in the analysis of the non hen, 1995; Gallavotti, 1996), which can be interpreted as an equilibrium systems is that the fluctuation-dissipation rela- extension of the ergodic hypothesis to non-equilibrium sys- tion (Kubo, 1966), cannot be applied (Ruelle, 1998a). This tems (Gallavotti, 2006). We also remind that any time in an greatly limits the ability of understanding the response of the numerical integration we assume that the time average of a systems to external perturbations by looking at its variabil- given observable, after discarding an initial transient, basi- ity. In the specific case of climate, this can be rephrased by cally coincides with its average over the invariant measure saying that climate change signals need not project on the giving the attractor of the system, we are actually assuming natural modes of climate variability. The non-equivalence Axiom A-like hypothesis. See Penland (2003) for an original between free and forced climate fluctuations had been sug- geophysical perspective. gested by Lorenz (1979). The basic reason for this behav- The Ruelle response theory, with the support of chaotic ior is that, since the dynamics is forced and dissipative, with hypothesis, has interesting conceptual implications for cli- the asymptotic dynamics taking place in a strange attractor, mate studies. In fact, the possibility of defining a response natural fluctuations and forced motions cannot be equiva- function basically poses the problem of climate response to Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 9 forcings and of climate change in a well-defined context, and, from the use of the fluctuation-dissipation theorem in a non- when considering the procedures aimed at improving climate equilibrium context. As discussed above, this just cannot be models, justifies rigorously the procedures of tuning and ad- rigorously the case. Since they use, albeit in a convolute way, justing the free parameters. Moreover, the response theory the forcing term to compute (or at least to parameterise) the allows to compute the climate sensitivity, in the special case response of the system, what they do is actually to compute when static perturbations to the system parameters are con- (or at least estimate) the linear response. Whereas it is in the sidered. very nature of a perturbation theory the possibility of com- Previously, a response formula was proposed by Cacuci puting the response starting only from the statistical proper- for evaluating the linearized change of the solution of a time- ties of the unperturbed state (as done also in Abramov and independent generic system of nonlinear equations as a result Majda, 2007; Majda and Wang, 2010), the specific property of a change in the system’s parameters (Cacuci, 1981a,b). of the systems where the fluctuation-dissipation theorem can This can be interpreted as a special case of Ruelle’s theory, be applied is that the linear Green function can be written in where the unperturbed attractor is constituted by a fixed point the specific form of a correlation of two observables of the and a static perturbation to the system evolution equation is system. See the discussion in Ruelle (1998a, 2009). considered. Cacuci proposed to study this problem using the It is also important to note that, following the pioneer- adjoint operator to the original system, which provided an ing study by Leith (1975), several recent studies (Langen efficient way to determine the impact of small perturbations. and Alexeev, 2005; Gritsun and Branstator, 2007; Ring and Interestingly, early prominent applications of the so-called Plumb, 2008; Gritsun et al., 2008) have attempted with a adjoint method and its extension to time-dependent prob- certain (and sometimes good) degree of success to recon- lems, which allowed for evaluating all possible linear sen- struct the climate response to external perturbations start- sitivities of an evolving model in just one simulation, were ing from the internal fluctuations of the system by using a been proposed for climate related problems. In particular, it severely simplified version of the fluctuation dissipation the- was used to evaluate the sensitivities of a simple radiative- orem, based upon the assumption of a quasi-gaussian proba- convective model possessing an attractor constituted by just bility distribution for the system (see discussion in Abramov one fixed point (Hall et al., 1982; Hall and Cacuci, 1983), and Majda, 2007). The overall good results seem to suggest and later, in an inherently heuristic way, for studying the re- that, at least for the considered models, observables (typi- sponse of a (chaotic) simplified general circulation model to cally very large scale ones) and baseline climate conditions, doubling of the CO concentration (Hall, 1986). Whereas the fractal nature of the invariant measure of the system and the adjoint method did not find much space in further cli- the role of the stable direction of the flow seem not to be ex- matic studies, mostly due to early discouragement for the ceedingly relevant. This may be related to the specific choice computational burden of constructing the suitable operators of the observable or to the fact that internally generated nu- for evaluating the sensitivities, it subsequently reached great merical noise mimics the effect of stochastic perturbations, success in data assimilation problems for geophysical fluid but further studies are surely necessary. dynamics (Ghil and Malanotte-Rizzoli, 1991; Errico, 1997), In the last decade on one side a great effort has been to the point that a tangent and adjoint model compiler able directed at extending the Ruelle response theory for more to automatically generate adjoint model code was has been general classes of dynamical systems (see, e.g., Dolgopyat, introduced (Giering and Kaminski, 1998). More recently, a 2004; Baladi, 2007), and recent studies (Lucarini, 2008b) link between advanced adjoint techniques and the Ruelle the- have shown that, thanks only to the causal nature of the re- ory has been proposed (Eyink et al., 2004). sponse, it is possible to apply all the machinery of the the- Majda and collaborators (Abramov and Majda, 2007; Ma- ory of Kramers-Kronig (KK) relations (Nussenzveig, 1972; jda and Wang, 2010) have taken a different approach for Peiponen et al., 2005; Lucarini et al., 2003, 2005) for linear analysing the response of a system to external perturbation. and nonlinear processes to study accurately and rigorously Using the formalism of the Fokker-Planck equation, they the susceptibilities describing in the frequency domain the have developed a response theory which is analogous to Ru- response of a general observable to a general perturbation. elle’s. These authors have especially emphasized the pos- Moreover, the actual applicability of the theory has been sibility of splitting the response into the components rela- successfully tested in a number of simple dynamical systems tive to the projection of the perturbation vector flow onto the case for the linear (Reick, 2002; Cessac and Sepulchre, 2007) stable, unstable and neutral components of the unperturbed and nonlinear (Lucarini, 2009b) response. Such numerical flow, as already discussed in Ruelle (1998a), and have in- investigations have clarified that even in systems which are troduced algorithms to compute efficiently all of these com- not Axiom A, like the Lorenz 63 system (Lorenz, 1963), it is ponents. Their approach has very recently been extended possible to successfully use the response theory to construct to take into account stochastic and time-periodic perturba- linear (Reick, 2002) and nonlinear susceptibilities (Lucarini, tions (Majda and Wang, 2010). Somewhat confusingly, they 2009b) which obey all of the constraints imposed by the KK mostly refer to the rather interesting applications of their the- theory up to a high degree of precision. ory and algorithms to specific dynamical systems as resulting www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 10 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings These investigations definitely motivate further studies Whereas the paper aims at proposing new methods for aimed at understanding to what extent the response theory tackling classical problems of climate science, most of the is an efficient tool for analyzing complex and complicated results and of the methodologies proposed are of more gen- systems. In this paper, we take up such a challenge and con- eral interest. In this paper we limit our attention to the linear sider the Lorenz 96 (L96) system (Lorenz, 1996; Lorenz and response. We refer to Lucarini (2008b, 2009b) for a theoreti- Emanuel, 1998; Lorenz, 2004), which provides an excellent cal and numerical studies of higher-order effects of perturba- and celebrated prototypical model of a one dimensional at- tions. mosphere. The variables of the L96 model can be thought The paper is organized as follows. In Sect. 2 we briefly as generic meteorological quantities extending around a lat- analyze the general theoretical background of the linear re- itudinal circle and sampled at a regular interval. In spite of sponse theory and of the properties of the frequency depen- not being realistic in the usual sense, the L96 model presents dent susceptibility and present some new useful results. In the basic ingredients, such as dissipation, advection and the Sect. 3 we present the main features of the L96 system, in- presence of an external forcing, of the actual atmosphere. For troduce the considered perturbations to the forcing, derive this reason, L96 has quickly become the standard model to some basic properties of the response of various observables, be used for predictability studies (Orrell, 2003; Haven at al., and present the theoretical predictions. In Sect. 4 we present 2006; Hallerberg et al., 2010), when testing data assimila- the results of our numerical investigations and describe how tion techniques (Trevisan and Uboldi, 2004; Trevisan et al., they can be generalized to the entire family of L96 models. 2010; Fertig et al., 2007), and new parameterizations (Wilks, In Sect. 5 we provide a relevant example to illustrate how the 2006). The L96 model had already been taken as test-bed for results presented in this paper can be used to devise simple studying the linear response (the applicability of the fluctua- yet rigorous methods to study the climate response at all time tion dissipation theorem, in their language) in Abramov and scales on models of any degree of complexity. In Sect. 6 we Majda (2007). discuss the conclusions and present perspectives for future Although we are unable to prove that the unperturbed L96 work. is an Axiom A system, in general and for the specific choice of parameters used in our numerical simulations in particular, 2 Theoretical background: Ruelle theory and we adopt the chaotic hypothesis and present the first thor- dispersion relations ough investigation of a spatially extended system by using the rigorous statistical mechanical methodologies presented 2.1 Definition of the linear susceptibility in Ruelle (1998a, 2009) and Lucarini (2008b, 2009b). More- over, since L96 is a spatially extended system, we also ex- We consider an Axiom A dynamical system described by the plore the applicability of the response theory in all possi- evolution equation x˙ = F (x), so that the invariant probabil- ble combinations of global/local observables and global/local ity measure ρ of the associated flow is of the SRB type perturbations. We compute rigorously the corresponding lin- (Ruelle, 1998a). Let h8i be the expectation value of the ear susceptibilities, verify the KK relations and the related general observable 8 defined as ρ (dx),8(x). We perturb sum rules, and find an empirical power law. This, as in the the flow of the system by adding a on the right hand side of case of discussed in Lucarini et al. (2007), supports the valid- the evolution equation a vector field X(x)f(t), where X(x) ity of the chaotic hypothesis, allowing to extend the results defines the pattern of the perturbation, and f(t) is its time obtained for our specific choice of model’s parameters to a modulation. The resulting evolution equation results to be rather general class of L96 systems. We also show how to x˙ = F(x) + X(x)f(t). Following Ruelle (1998a), we ex- go from the frequency back to the time domain, thus deriv- press the expectation value of 8(x) in the perturbed system ing from the susceptibility the Green function, which acts as using a perturbative expansions as: time propagator of the considered perturbation for the con- sidered observable. The Green function allows to predict, X (n) h i h i h i 8 (t) = 8 + 8 (t). (1) in an ensemble mean sense, the change in the observable at 0 n=1 any time horizon as a result of a perturbation with the same spatial patter as that considered in the calculation of the sus- Each term of the perturbative series can be expressed as an ceptibility but featuring a general time modulation. th n-convolution integral of the n order causal Green function Finally, building upon the results presented here, we pro- with n delayed time perturbation functions (Ruelle, 1998b; pose a simple yet general methodology to study general Cli- Lucarini, 2008b). Limiting our attention at the linear case mate Change problems on virtually any time scale by resort- we have: ing to only well selected simulations, and by taking full ad- +∞ vantage of ensemble methods. The specific case of globally (1) (1) h8i (t) = dσ G (σ )f(t −σ ) (2) 1 1 1 averaged surface temperature response to a general pattern of −∞ change of the CO concentration is discussed. Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 11 The first order Green function can be expressed as follows: Green function is causal. Assuming, on heuristic physical ba- (1) sis, that G (t) ∈ L , we can apply the Titchmarsh theorem (1) 8 G (σ ) = ρ (dx)2(σ )35 (σ )8(x), (3) 1 0 1 0 1 8 (Nussenzveig, 1972; Peiponen et al., 2005; Lucarini et al., (1) 2003, 2005) and deduce that the linear susceptibility χ (ω) where 3(•) = X(x)·∇(•) takes into account the effects of is a holomorphic function in the upper complex ω-plane and the perturbative vector field, 2 is the usual Heaviside dis- the real and the imaginary part of χ(ω) are connected to each tribution and 5 the unperturbed time evolution operator so other by Hilbert transform. that 5 K(x) = K(x(t)) for any function K , with x(t) fol- According to a general property of Fourier transform we lowing the unperturbed flow. Note that it is possible to ex- (1) know that the short term behavior of G (t) determines the press the Green function as the expectation value of a non- (1) asymptotic properties of χ (ω). We shall obtain a more trivial but computable observable over the unperturbed SRB quantitative result by exploiting that: measure ρ . Therefore the knowledge of the unperturbed fea- tures of the flow is sufficient to define the effects of any ex- Z +∞ k d i k k ternal perturbation over any observable of our system. In the dt2(t)t exp[iωt] =(−i) P +πδ(ω) dω ω −∞ frequency domain we find that the first term of the perturba- (k+1) tive series can be written as: i Z ≈k! (6) +∞ (k+1) (1) (1) h8i (ω) = dω χ (ω )f(ω )×δ(ω−ω ) 1 1 1 1 −∞ where in the second equality we have neglected the fact that (1) the solution is a distribution and considered ω 6= 0. There- =χ (ω)f(ω), (4) fore, if the Taylor expansion of the Green function in the where the Dirac delta implies that we are analyzing the im- limit t → 0 is of the form: pact of perturbations in the frequency-domain at the fre- (1) quency ω. The linear susceptibility is defined as: β β G (t) ≈α¯ 2(t)t +o(t ) (7) +∞ (1) (1) χ (ω) = dtG (t)exp[iωt]. (5) the high frequency behavior of the linear susceptibility, i.e. 8 8 −∞ the limit ω → ∞, is: It is important to underline with a thought experiment the (1) −β−1 −β−1 computational relevance of the last equations and the impor- χ (ω) ≈αω +o(ω ) (8) tance of the susceptibility function. (β+1) Let’s suppose we introduce a time dependent perturbation where α =α¯i β!. The parameters β (which is an integer f (t) to a given pattern of forcing X(x), simulate the sys- α number) and α¯ depend on the observable 8, on the specific tem and observe the time response of an arbitrary observ- features of the unperturbed system, and on the forcing un- (1) able h8 i (t). We now compute the Fourier transform of der consideration. Taking into account Eq. (5) and assuming (1) (1) the observed signal and of the forcing modulation. Invert- that ω is real, we obtain that χ (ω) = [χ (−ω)] , so that 8 8 (1) ing Eq. (4), we can find the linear susceptibility χ (ω). Re[χ] is an even function while Im[χ] is odd function of ω. Let’s now consider a different time-modulating function of Thus α =α is real if β is odd, whereas α = iα is imaginary R I the forcing f (t) and its corresponding Fourier transform β if β is even. f (ω). Taking into account Eq. (4), if we multiply f (ω) β β Taking into account the Titchmarsch theorem, using that (1) (1) (1) times the previously computed function χ (ω) we directly χ (ω) = [χ (−ω)] , and considering the asymptotic be- ¯ 8 8 (1) havior of the susceptibility, it is possible to show that the real obtain h8 i (ω), the frequency-dependent response of the and imaginary part of the linear susceptibility obey the fol- observable 8 to the forcing X(x) modulated by the new lowing set of general KK dispersion relations (Lucarini et al., function. By applying the inverse Fourier transform we ob- (1) 2005): tain the time-dependent response h8 i (t) without needing any additional simulation. (1) 2p π ν Re[χ (ν)] (1) Moreover, the knowledge of the susceptibility function al- 2p−1 8 − ω Im[χ (ω)] =P dν (9) (1) 2 2 2 ν −ω lows us to reconstruct the G (t) by inverting Eq. (5). Oth- 0 erwise, the Green function can be obtained directly from observing the response signal by performing a simulation (1) ∞ 2p+1 π ν Im[χ (ν)] where f(t) =δ(t): in this case (see Eq. (2)) we simply have (1) 2p 8 ω Re[χ (ω)] =P dν (10) (1) 8 (1) 2 2 2 ν −ω h8i (t) =G (t). with P indicating integration in principal part and p = 2.2 Kramers-Kronig relations and sum rules 0,...,(β − 1)/2 if β is odd and p = 0,...,β/2 if β is even. As we see from Eq. (3), for an arbitrary choice of the ob- Note that the faster the asymptotic decrease of the suscepti- servable and of the perturbation the corresponding linear bility, the higher the number of independent constraints due www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 12 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings to KK relations it has to unavoidably obeys. As thoroughly discussed in Kubo (1966), in the case of quasi-equilibrium (1) system, the fluctuation-dissipation theorem ensures that the ν Re[χ (ν)]dν = α , (15) imaginary part of the susceptibility describing the response of a given observable to a perturbation is proportional to where the α constants are defined in Eq. (8). These sum rules a suitably defined power spectrum in the unperturbed sys- provide additional general constraints that must be obeyed by tem. Therefore, observing the unperturbed system and using any system and can be used to test the quality of the output of Eq. (10) it is possible to reconstruct the entire linear suscepti- any model wishing to describe it. If we are not in the perfect bility, and so know everything about the response properties model scenario (e.g., we use a simplified representation of of the system. In the case of a non-equilibrium system, as some degrees of freedom) the sum rules can in principle be discussed in the Introduction, this procedure is not possible. used to provide a fit for the parametrization. It is possible to use the KK relations to define specific self- We underline that it is possible to generalize the KK the- consistency properties of the real and imaginary part of the ory for specific classes of nonlinear susceptibilities for both susceptibility. We first consider the following application: quasi-equilibrium and non-equilibrium systems. Such re- we set p = 0 in Eq. (10) and take the limit ω → 0. We obtain sults, which are particularly suited for studying the funda- that for any observable: mental properties of harmonic generation processes, are thor- oughly discussed in Lucarini (2008b) and will not be re- (1) 2 Im[χ (ν)] (1) Re[χ (0)] = P dν , (11) ported here. π ν 2.3 A practical formula for the linear susceptibility and which says that the static susceptibility (i.e., in a more com- consistency relations between susceptibilities of mon language, the linear sensitivity of the system) is related different observables to the out-of-phase response of the system at all frequencies. In other terms, Eq. (11) is an exact formula for the linear sus- As discussed above, the definition of the linear susceptibility ceptibility of the system. Note that the static susceptibility does not depend on the function f(t) modulating the addi- is a real number because, thanks to the symmetry proper- (1) tional forcing, so that it is possible to draw general conclu- ties discussed above, Im[χ (0)] = 0. Moreover, we know (1) sions on its properties even by choosing a specific function that Re[χ (0)] is finite because the susceptibility function f(t). is analytic (and so in particular non singular). This is con- Let’s consider f(t) = 2 cos(ωt). The impact of the per- sistent with the fact that as ω → 0, the imaginary part of the turbation on the evolution of a general observable 8(x) is susceptibility goes to zero at least as fast as a linear func- defined as: tion, as only odd positive integer exponents can appear in its Taylor expansion around ω = 0), so that the integrand in δ8 (t,t ,x ) =8 (t,t ,x )−8 (t,t ,x ) (16) 0 0  0 0 0 0 0 Eq. (11) is not singular. Similarly, we obtain that for ω ∼ 0 the real part of the susceptibility is in general of the form where x and t are the initial condition and the initial 0 0 2 3 c +c ω. +o(ω ), where c and c are two constants and c 1 2 1 2 1 time, and we associate the lower index  to the strength is exactly given by Eq. (11). of the forcing. The Ruelle’s response theory ensures that By exploring the ω → ∞ limit in Eqs. (9–10) we obtain δ8 (t,t ,x ) = O(). Following Reick (2002) and Lucarini 0 0 further integral constraints. By applying the superconver- (2009b) the linear susceptibility results to be: gence theorem (Frye and Warnock, 1963), we obtain the fol- (1) (1) lowing set of vanishing sum rules (see Lucarini et al., 2003): χ (ω) ≡ lim lim χ (ω,x ,,T ) (17) 8 φ →0T →∞ 0 ≤p ≤β/2− 1, β even (1) 2p+1 dνν Im[χ (ν)] = 0 . (12) 8 where: 0 ≤p ≤(β − 3)/2, β odd 1 1 (1) ∞ χ (ω,x ,,T ) = dtδ8 (t,x )exp(iωt) (18) 0  0 0 ≤p ≤β/2− 1, β even (1) 2p ν Re[χ (ν)]dν = 0 . (13) 0 ≤p ≤(β − 1)/2, β odd is the total susceptibility, affected by the finite time and finite Note that if β = 0 no vanishing sum rules can be written for size response of the system. This quantity depends on the the susceptibility, whereas if β = 1 only Eq. (13) provides a initial condition and in principle contains information about zero-sum constraint. For each set of KK relations, an addi- the response of the system at all order of nonlinearity. tional, non-vanishing sum rule can be obtained. If β is odd, Since d(δ8 (t,x ))/dt = δ(d8 (t,x ))/dt , thanks to the 0  0 the non-vanishing sum rule is: linearity of the time derivative, by considering Eqs. (17–18) and performing an integration by parts, we obtain that (Re- (1) ν Im[χ (ν)]dν = − α , (14) ick, 2002): (1) (1) while if β is even, we have: χ (ω) = −iωχ (ω). (19) ˙ 8 Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 13 (1) Let’s now find a different expression for χ (ω). The time Moreover, since, as shown in Eq. (7–8), the susceptibility derivative of 8 in the unperturbed system is of a generic observable decreases to zero at least as fast as (1) −1 ω , for large values of ω we have that χ (ω) ≈ i/ωh4i 8(x) =0(x), (20) unless ωh4i = 0. If ωh4i 6= 0, we also have that the lead- 0 0 ing order of the short-time expansion of the Green function where 0 = F ·∇8. Similarly, the time derivative for 8 in is of the form: the case of the perturbed motion described by x˙ = F(x)+ X(x)f(t) = F(x)+ 2 cos(ωt)X(x) is: (1) G (t) =2(t)h4i +o(t ), (26) 8(x) =0(x)+ 2 cos(ωt)4(x), (21) in agreement with what can be found by direct inspection of Eq. (3). where 4 = X·∇8. From Eqs. (20–21) we obtain that δ8 (t,x ) =δ0 (t,x )+ 2 cos(ωt)4(t,x ), (22) 0  0 0 3 Application of the response theory to the Lorenz where all terms are of O(). Furthermore, we integrate each 96 model term in Eq. (22) as in Eq. (18), take the limits  → 0 and T → ∞, and obtain: 3.1 Statistical properties of the unperturbed Lorenz 96 Model 1 1 lim lim dtδ8 (t,x )exp(iωt) T →∞ →0 T 0 The Lorenz 96 model (Lorenz, 1996; Lorenz and Emanuel, 1 1 1998; Lorenz, 2004) describes the evolution of a generic at- = lim lim dtδ0 (t,x )exp(iωt) mospheric variable defined in N equally spaced grid points →0T →∞T along a latitudinal circle and provides a simple, unrealistic 1 1 + lim lim dt4 (t,x )[ exp(iωt) but conceptually satisfying representation of some basic at- →0T →∞T mospheric processes, even if such one-dimensional model it +exp(−iωt) ]exp(iωt). (23) cannot be derived ab-initio from any dynamic equation via subsequent approximations. The evolution equations can be Using the definition in Eq. (17) and the identity given in written in a scaled form as follows: Eq. (19), we derive: dx =x (x −x )−x +F (27) i−1 i+1 i−2 i (1) (1) dt −iωχ (ω) =χ (ω))+ lim lim dt4 (t,x ) 8 0 →0T →∞T where i = 1,2,.....,N , and the index i is cyclic so that x = Z i+N 1 1 x = x . The quadratic term in the equations simulates i−N i + lim lim dt4 (t,x )exp(2iωt). →0T →∞ advection, the linear one represents thermal or mechanical (24) damping and the constant one is an external forcing. The evolution equations are invariant under i →i + 1, so that the The first limit in Eq. (24) gives, by definition, h4i , whereas dynamics is the same for all variable. The time scale of the the second limit vanishes as the expression under integral is 2 system is given by the damping time, which corresponds to O( ), since it is related to second order harmonic gener- five days. The L96 system shows different features, as dif- ation nonlinear process (Lucarini, 2009b). Concluding, we ferent choices of F and N may strongly alter the topology of obtain the following general consistency relation for the lin- the attractor, alternating periodic, quasi-periodic and chaotic ear susceptibility: behavior in a non trivial way. However with a suitable choice i i (1) (1) of the parameters N and F , the system is markedly chaotic. χ (ω) = χ (ω)+ h4i . (25) 8 0 ω ω In particular, as F controls the energy input into the system, we expect that for relatively high values of this parameter Such an identity related the susceptibility of an observable the system should simulate a turbulent behavior and live on 8 to the susceptibility of the projection of its gradient along a strange attractor. As an example, setting N = 40 and F = 8 the unperturbed flow 0 and to the average value in the unper- the system possesses 13 positive Lyapunov exponents, the turbed state of the projection of its gradient along the pertur- largest corresponding to a doubling time of 2.1 days, while bation flow. Note that the two terms on the right hand side the fractal dimension of the attractor (Ruelle, 1989) is about are radically different. Whereas the first term is related to the 27.1 (Lorenz and Emanuel, 1998). projection of the dynamics along the unstable manifold, the When computing the time derivative of the total energy of second term depends on the structure of the forcing X(x), P the system, defined as E = 1/2 x , the advection terms i i which may be entirely unrelated to that of the unstable man- cancel. The evolution equation for E results to be: ifold. This is the fundamental reason why the fluctuation- dissipation theorem does not apply in the non-equilibrium E = −2E +F x . (28) case. www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 14 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings As the dynamics takes place inside a compact set, 9(x ,...,x ) is bounded for any choice of the function 9 . i N 1+γ Therefore the ensemble mean with respect the ρ (or time e(F) = F , (36) average) of the temporal derivative 9 vanishes. Therefore, defining M = x as the total momentum of the system, with λ ≈ 1.15 and γ ≈ 0.35. Such a smooth dependence of we obtain the following identity: the intensive energy and momentum with respect to the forc- D E X X ing parameter F is indeed in agreement with the hypothesis 2hEi = x =F hx i =F hMi . (29) 0 i 0 0 that the invariant measure is deformed in a very regular fash- i i ion not only locally, but over a large range of the parameter’s Similarly, we can deduce an additional consistency relation space. by investigating the expression of the time derivative of M : Note that, at the fixed point of the system corresponding to a purely zonally symmetric dynamics (x = F , ∀i) we have hMi =NF +hC i −hC i (30) 2 3 0 0 0 m =F and e =F /2. These formulas give much higher val- P P where hC i = hx x i and hC i = hx x i . 2 i i−2 3 i i−3 0 0 0 ues for both m and e than what found with our empirical i i Higher order consistence relations can be obtain in a simi- power laws for the attractor in the chaotic regime. In fact, lar fashion. at such an equilibrium, which is unstable in the parametric The equivalence of all the variables implies that over the range explored here, the energy dissipation is much weaker unperturbed flow each observable O of the whole system sat- than in the co-existing chaotic attractor, which corresponds to isfies hO(x )i = N O(x ) ∀j . Therefore, we define i j the case where breaking nonlinear waves and turbulent mo- i 0 the average energy per grid point e(N,F) and the average tions are present. Interestingly, the presence of well-defined momentum per grid point m(N,F) as: scaling laws with respect to the forcing parameters for the 2 energy and momentum of the system with different charac- hEi i 0 e = = , (31) teristic exponents in the chaotic regime and in the co-existing 2 N unstable equilibrium is in agreement with previous finding recently obtained in a simple baroclinic quasi-geostrophic hMi model (Lucarini et al., 2007). m = hx i = , (32) i 0 3.2 Asymptotic properties of the linear susceptibility where the choice of i is arbitrary and the N and F - dependence is dropped for shortness. Defining c = We perturb the L96 model by adding a small perturbation hC i /N , c = 2e, we can rewrite Eqs. (29–30) as follows: i 0 0 modulated by f(t) = 2 cos(ωt). The resulting evolution equation is: 2e =F m, (33) dx =x (x −x )−x +F + 2 cos(ωt)X (37) i−1 i+1 i+2 i i dt m =F +c −c . (34) 2 3 where X = X (x ,...,x ) is a generic function of vari- i i 1 N Expressing either the average energy e or the average mo- ables x . We adopt the chaotic hypothesis (Gallavotti, 1996; mentum m per grid point as a function of the two free pa- Gallavotti and Cohen, 1995) and we follow the theory pro- rameters N and F would allow to get a closure equation for posed by Ruelle (1998a) and discussed in Sect. 2 in order the statistical properties of the unperturbed Lorenz 96 sys- to study the linear response of suitably defined observables tem. We have computed e(N,F) and m(N,F) by performing to the perturbation. We first propose to study the high- long integrations for values of F ranging from 6 to 50 with frequency, response by analyzing in detail the asymptotic step 1 and for values of N ranging from 10 to 200 with step properties of the resulting susceptibilities. As discussed in 10. In all of these cases, chaotic motions are observed. We Sect. 2, this constitutes a crucial step for constructing the set have consistently found that, within 0.5%, e(N,F) = e(F) of applicable KK relations and for computing the value of the and m(N,F) = m(F), so that they can be considered inten- sum rules. sive quantities. Therefore, we can interpret Eqs. (33–34) as We consider two different forcing patters X . In the first equations providing a definition of the thermodynamics of case, we apply the perturbation over all the grid points and this simple one-dimensional model of atmosphere. we choose X = 1∀i. Given an observable 8, we refer to the The F -dependence of e and m can be closely approxi- linear Green function and linear susceptibility resulting from (1) (1) mated in terms of simple power laws. We obtain, within a this choice of X as G and χ , respectively, where the 8,N 8,N precision of about 1% in the considered domain, that lower index N indicates that the perturbation acts over all γ the variables x . We refer to this pattern of forcing as global m(F) =λF , (35) perturbation. In the second case, we apply the perturbation and, consistently with Eq. (33), only on the variable x of the L96 model, and we choose Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 15 X = 0∀i 6= j , X = 1. Since all the points are equivalent in Equations (40) and (41) imply that the imaginary part 1 j the unperturbed case, the choice of j is arbitrary. In this case, dominates the asymptotic behavior of the susceptibility, so when referring to the linear Green function and the linear that at high frequency the response is shifted by about π/2 susceptibility, the lower index 1 substitutes the N , indicating with respect the forcing. Observing that the leading term of −1 that the perturbation is localized to one point. We refer to asymptotic χ is of order ω just one sum rules apply for this pattern of forcing as local perturbation. either susceptibilities. Limiting our attention to the intensive quantity e, by applying Eq. (15) we obtain: 3.2.1 Global perturbation (1) Re[χ (ω)]dω = m. (42) ε,N We consider perturbations with spatial pattern given by X = i 2 1∀i and analyze the response of the observable E. Following Along similar lines, if we select as observable the total mo- (1) Eq. (3), the linear Green function G (t) can be explicitly E,a mentum M , we derive that the asymptotic behavior of its lin- written as: ear susceptibility is: (1) G (t) = ρ (dx)2(t)35 (t)E(x) 0 0 N N E,N (1) (1) −2 χ (ω) =Nχ (ω) = i + +o(ω ), (43) M,N μ,N Z 2 ω ω = ρ (dx)2(t)1·∇E(x(t)) (1) where we have defined the intensive susceptibility χ (ω), Z μ,N where μ is the intensive momentum of the system. As in = ρ (dx)2(t) ∂ (E(x(t))) (38) 0 i Eq. (41), the asymptotic behavior is determined by the imag- inary part of χ , and the real part of the susceptibility provides where x(t) satisfies the unperturbed evolution Eq. (27). Tak- the following sum rule: ing into account Eq. (7) and Eq. (8), in order to obtain the (1) asymptotic behavior of the susceptibility, we need to study Re[χ (ω)]dω = . (44) μ,N the short time behavior of the Green function. Therefore, we express E(x(t)) as a Taylor series about t = 0 considering We now wish to go back to the general consistency equation the unperturbed flow, compute the integral of each coeffi- for linear susceptibilities given in Eq. (46). Considering that cient of the t -expansion over ρ , and seek the lowest order in the perturbed system the time derivative of the total energy non-vanishing term (Lucarini, 2009b). The first two terms of of the system can be written as: the Taylor expansion of E in Eq. (38) give: E = −2E +FM + 2 cos(ωt)M (45) (1) G (t) = ρ (dx)2(t) ∂ E| +tE| +o(t) 0 i t=0 t=0 E,N i the general result given in Eq. (25) can be written as follows: h   i X X = ρ (dx)2(t) x − 2x −NF t +o(t) . (39) 0 i i F 1 i i (1) (1) χ (ω) = χ (ω)+ hMi , (46) E,N M,N 2− iω 2− iω Using Eqs. (7–8), the leading terms of the asymptotic behav- ior of linear susceptibility can be written as: since in this case 0 = −2E +FM and 4 = M . It is easy X X to check that the asymptotic behavior for the susceptibilities (1) χ (ω) =i hx i /ω+ h2x i −NF /ω i i 0 0 given in Eqs. (40–43) is in agreement with Eq. (46), which is E,N i i valid at all frequencies. m F − 2m −2 −2 +o(ω ) =iN −N +o(ω ) (40) 3.2.2 Local perturbation ω ω Since the symmetry with respect the index i is valid also in We now perturb the system in a single grid point. The sym- the perturbed case, given our choice of the forcing pattern, metry of the unperturbed system implies the equivalence of the linear susceptibility of the total energy is given the sum every point of the latitude circle. Indicating with x the grid of N identical contributions, each corresponding to the sus- point where forcing is exerted, the pattern of the perturbation ceptibility of the observable ε = 1/2x for each of the N vari- vector field is X = 0∀i 6= j , X = 1. We consider the same 1 j ables x of the system. Therefore, it is possible to define an modulating monochromatic function f(t) = 2 cos(ωt) as in (1) (1) (1) intensive linear susceptibility χ = 1/Nχ , where χ the previous case. We analyze the asymptotic behavior of ε,N E,N ε,N describes the response of the local energy to the external per- the linear susceptibility of the total energy of the system E, turbation. In particular in the limit ω → ∞ we have: of the total momentum of the system M , which are global variables, and of the local variables constituted by the energy m F − 2m (1) −2 χ (ω) =i − +o(ω ). (41) E = 1/2x and momentum M = x of the perturbed grid j j j ε,N 2 j ω ω point and of its immediate neighbors. www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 16 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings Since we are looking at the linear response and the we obtain the following asymptotic behavior for the linear global perturbation is given by N spatially shifted copies susceptibility of the local perturbation, for any observable 8 of the form P m m (1) (1) (1) N −2 8(x ,...,x ) = φ(x ) we must have χ = χ = χ =i + +o(ω ). (54) 1 n i i=1 8,1 8/N,N E ,1 j 2 ω ω (1) χ . φ,N (1) (1) In the case of the observable E, it is straightforward to ver- Since, by linearity, χ = χ , comparing this result E,1 E ,1 ify the previous identity at least in the asymptotic regimes. In with what obtained in Eq. (48), we note that the suscepti- fact, the short time behavior of the Green function describing bility of the energy at the position of the forcing E pro- the response of the E to the local perturbation results to be: vides the leading asymptotic term to the susceptibility of the total energy E. Consequently, in the high-frequency range (1) (1) (1) G (t) = ρ (dx)2(t)∂ E| +tE| +o(t) 0 j t=0 t=0 χ ≈χ , and the two susceptibilities obey the same non- E,1 E ,1 E,1 vanishing sum rule, so that: h   i X X 2 2 = ρ (dx)2(t)∂ x − 2x −x F t 0 j i i i (1) Re[χ (ω)]dω = m (55) E ,1 =2(t) 2x − 2x −F t +o(t) , (47) j j 0 0 0 Nevertheless, by comparing Eqs. (48–54), we discover that so that the asymptotic behavior of the corresponding linear −2 contributions to the second leading order (∝ω ) in the high susceptibility susceptibility is: frequency range of the susceptibility of the total energy do m F − 2m not come just from the response of the energy at the per- (1) −2 χ =i − +o(ω ), (48) E,1 turbed grid point. perturbed grid point but some other point ω ω −2 give a contribution of order ω . Therefore, the asymptotic which agrees with what found for the intensive energy re- (1) (1) behavior of the real part of χ is not captured by χ . sponse when the global perturbation is applied (see Eq. (41)). E,1 E ,1 The locality of the interaction suggests to look at the energy The sum rule for the real part of the susceptibility is exactly of the closest neighbors of x . Because of the asymmetry of the same as in what given in Eq. (42): j the nonlinear terms in the L96 evolution equations, we con- (1) sider the observable ψ = 1/2(E +E +E ). It is j+1 j+2 j−1 Re[χ (ω)]dω = . (49) E,1 2 possible to prove that: Analogously, we obtain that the asymptotic behavior of x (x −x ) (F −m) j−1 j+1 j−2 (1) 0 −2 (1) χ (ω) = − = − +o(ω ). (56) χ can be written as: ψ,1 2 2 M,,1 ω ω (1) (1) 1 1 (1) It is easy to observe that the sum of χ and χ provides ψ,1 E ,1 χ (ω) = i + , (50) j M,1 ω ω the correct leading order to the asymptotic behavior of both (1) the real and imaginary parts of χ . We shall provide an with the corresponding sum rule: E,1 argument why this strongly supports the close resemblance π (1) (1) (1) of the two functions χ and χ , where E =E +ψ = Re[χ (ω)]dω = , (51) loc j E,1 E ,1 loc M,1 0 2 2 2 2 1/2x + 1/2x + 1/2x + 1/2x is the energy of the j j+1 j+2 j−1 cluster of points centered in x . in perfect agreement with Eqs. (43) and (44), respectively. j The analysis of the asymptotic behavior of the susceptibil- It is rather interesting to look into local energy observ- ities related to the local momentum of the system provides ables. Considering the energy E of the perturbed grid point additional insights. It is possible to prove that for large fre- x we have that its short term Green function can be written quency the linear susceptibility of the momentum of the per- as: turbed grid point is: (1) G (t) = ρ (dx)2(t)∂ [ x +(x (x x 0 j j j−1 j+1 E ,1 1 1 1 (1) −2 χ (ω) = i + +o(ω ), (57) x ,1 −x x −x +F ) )t +o(t) ] =2(t)[ x j j−1 j−2 j j ω ω − x x −x x −x +F t +o(t) ]. (52) j−1 j+1 j−1 j−2 j which suggests that the response of the local momentum cap- tures the correct asymptotic behavior of both the real and the Since imaginary part of the total momentum M . Concluding, we 0 = x˙ = x x −x x −x +F (53) j j−1 j+1 j−1 j−2 j obtain that the following sum rule can be stablished: 0 0 because the ρ -average of the temporal derivative of any ob- 0 π (1) Re[χ (ω)]dω = . (58) x ,1 servable vanishes, thanks to the compactness of the attractor, j Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 17 Therefore, such a constraint is exactly the same whether we t = T = 400π/ω, which corresponds to 200 full periods of analyze the the response of the momentum of a single vari- the forcing. The length of the simulations depends on the able when the perturbation acts over all the grid points, or of corresponding period of the forcing because we are inter- the total momentum in the case of a local perturbation, or, in ested in obtaining a frequency-independent quality for the this latter case, of the momentum of the grid point where the signal. We have observed that the linear response approxi- local perturbation is applied. mation is obeyed to a good degree of approximation for up As we have seen in this section, the coefficients of the lead- to  ≈ 1, which implies that the third order nonlinear effects ing asymptotic terms and the sum rules are given by simple are relatively small. See Lucarini (2008b, 2009b) for further linear functions of m (or equivalently, thanks to Eq. (33), by clarifications on this latter point. e) and by F . As we have proposed an efficient parameteri- When considering the susceptibilities describing the re- zation of m and e as functions of F alone in Sect. 3.1, our sponse to the global perturbation, we present results obtained results can be easily applied and numerically verified for a using  = 0.25 and averaging over K = 100 random initial very large class of L96 models. conditions. When assessing the linear response to the local perturbation, a reasonably clear signal is obtained using  = 1 and averaging over K = 300 initial conditions. 4 Results Note that, since we are interested in the linear response, it is could have been possible to compute the susceptibility 4.1 Simulations and data processing using a generic modulating function f(t) (see Eq. 4) rather than having to resort to multiple monochromatic perturba- The accurate calculation of the linear susceptibility of the tions. Nevertheless, for reasons of clarity, and for empha- general observable 8 is not as easy task, since the defini- sizing that chaotic dynamical systems can be analyzed using tion given in Eq. (17) requires the evaluation of two limits, tools typical of spectroscopy, we have used a more cumber- whereas we can actually compute only the quantity given in some but probably more convincing approach. Eq. (18). Averaging the response over a long time T allows We underline that the numerical results have been obtained for improving the signal-to-noise ratio. Noise is present be- using a commercial laptop rather than resorting to HPC. This cause, due to the chaotic nature of the flow, we have a contin- comes from the motivation of showing that the methodology uous spectral background. Instead, considering small values presented is robust enough that relatively low-end means al- for the perturbation strength  degrades the signal-to-noise low us to see the physical and mathematical properties of ration, but, on the other hand, it is crucial to select a small our interest. We emphasize that, using HPC, it is rather easy in order to keep the perturbations as close as possible to to greatly increase the quality of the signal by increasing K the linear regime. As discussed in Lucarini (2009b), we can and/or T by a one or two orders of magnitude. improve the signal-to-noise ratio without needing to perform very long integrations and to consider large values for  by 4.2 Global perturbation performing an ergodic averaging of the quantity averaging (1) the quantity χ (ω,x ,,T ). Therefore, we choose the best i (1) (1) We first consider χ = 1/Nχ , where ε =E/N , and fol- (1) ε,N E,N estimator of the true susceptibility χ as: low up from the discussion in Sect. 3.2.1. The measured real and imaginary parts of the susceptibility are depicted with the black lines in Fig. 1a, b. The imaginary part has a broad spec- (1) (1) χ (ω) = lim χ (ω,x ,,T ), (59) 8 8 tral feature (with two distinct peaks) spanning from ω ≈ 2 to K→∞K i=1 ω ≈ 4, which corresponds to about twice the time scale (= 1) where the x are randomly selected initial conditions chosen of the system and to four times (see Eq. 28) the relaxation on the attractor of the unperturbed system. time of the energy. This hints at the fact that it is not ob- The numerical integrations of the Lorenz 96 system have vious to constrain the spectral features of the response an been performed using the standard configuration where N , observable just by performing a scale analysis of its evolu- the number of degrees of freedom, is set to 40, and F , the tion equation. For higher values of ω, the imaginary part intensity of the unperturbed forcing, is set to 8 (Lorenz, 2004, decreases in a very regular way, so that in the upper range 1996). Equations (27) and (37) are solved using the standard a very good agreement with the asymptotic behavior ∼m/ω fourth order Runge-Kutta numerical scheme. presented in Eq. (8) is obtained. For low frequencies, the For a given observable 8, the susceptibility at angular fre- imaginary part appears to decrease towards zero, as expected quency ω is computed by applying Eq. (59) to K outputs from symmetry reasons. Instead, the real part presents a dis- of Eq. (37), each starting with a different initial condition, persive structure in correspondence with the broad maximum where the perturbation has the same angular frequency ω. of the imaginary part, and changes sign for ω ≈ 6, so that it is The angular frequency ω ranges from ω = 0.2π to ω = 20π negative for high values of the frequency range. The asymp- l h with steps of 0.01π . Each simulation performed with a totic decrease to zero in this range is also in excellent agree- perturbation of angular frequency ω runs from t = 0 up to ment with the estimate ∼ −(F − 2m)/ω given in Eq. (8), www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 18 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 0.3 Measurement Measured a) a) Extrapolation Extrapolation KK KKmea mea KK KK ext ext 0.25 1.5 0.2 0.15 0.1 0.5 0.05 −0.05 −2 −1 0 1 2 3 −0.5 −2 −1 0 1 2 3 10 10 ω 10 10 ωh 10 10 ω l 10 10 l 10 10 ω 10 10 h ω 0.2 1 Measurement b) Measured Extrapolation b) Extrapolation KK mea KKmea KK ext KK ext 0.15 0.1 0.5 0.05 −0.05 −0.1 −2 −1 0 1 2 3 ωl h −0.5 10 10 10 10 10 10 −2 −1 0 1 2 3 10 10 ωl 10 10 10 10 Fig. 2. Linear susceptibility of intensive momentum of the system Fig. 1. Linear susceptibility of intensive energy of the system E/N M/N with respect to the global perturbation with X = 1 ∀i. The with respect to the global perturbation with X = 1 ∀i. The real and real and the imaginary parts are depicted in (a) and (b), respectively. the imaginary parts are depicted in (a) and (b), respectively. The The measured and extrapolated values are shown in red and black measured and extrapolated values are shown in red and black lines, lines, respectively. The result of the Kramers-Kronig inversion done respectively. The result of the Kramers-Kronig inversion done with with the measured and with with the extrapolated data are shown in the measured and with with the extrapolated data are shown in blue blue and magenta lines, respectively. and magenta lines, respectively. totically for high frequencies as 1/ω and 1/ω, respectively. whereas for low frequencies the real susceptibility tends to a We apply the truncated KK relations to the measured data very high value, this suggesting that the strongest response is to test the quality of the data inversion process. The estimates obtained for static perturbations. of the imaginary part (starting from the measured data of the (1) The measured real and imaginary parts of χ = real part) and of the real part (starting from the measured data μ,N (1) of the imaginary part) obtained by applying Eqs. (9–10) are 1/Nχ , where μ =M/N , are depicted in black in Fig. 2a, M,N (1) (1) shown for χ in blue in Fig. 1a, b and for χ in Fig. 2a, b. Interestingly, the spectral feature of the imaginary part is ε,N μ,N b. We observe that whereas agreement is very good for the shifted to higher frequencies than in the case of the energy real part for both susceptibilities, only a qualitative match is susceptibility, so that a well-distinct peak centered on value obtained for the imaginary part, with large discrepancies for of ω ≈ 6, which approximately corresponds to the natural ω. 2. In this latter case, moreover, the well-known problem time scale of the system. For low frequencies, the suscep- of KK divergence at the boundaries of integration (Lucarini tibility has almost exclusively a real component. As opposed et al., 2003, 2005) is very serious for ω =ω . to the previous case, the largest value for the in-phase re- l sponse is not obtained for ultralow frequencies, but rather for It is crucial to test whether the discrepancies are due to the ω ≈ 4. The asymptotic behavior of both the real and imagi- finiteness of the spectral range or are, instead, due to basic nary parts is in perfect agreement with the theoretical result problems in the applicability of the Ruelle response theory, given in Eq. (32), so that they are found to decrease asymp- related to the fact that the invariant probability measure of the Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ Im[χ] Re[χ] Im[χ] Re[χ] V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 19 unperturbed system actually features large deviations from They correspond, by definition, to the static susceptibility an SRB measure. of the observables e and m, respectively, for the global We proceed testing the first case. In order to widen the perturbation with X = 1 considered here. When evaluat- spectral range over which the susceptibility is defined, we ing the derivatives of e(F) and m(F) for F = 8 we obtain will exploit the asymptotic properties obtained in Sect. 3.2 (de(F)/dF) ≈ 1.6 and (dm(F)/dF) ≈ 0.11. These F=8 F=8 as well as the low frequency behavior of the susceptibility values are in good agreement with what found by extrapo- discussed in Sect. 2. We redefine the the imaginary part of lating the corresponding real part of the susceptibilities for (1) ω → 0 via KK relations and shown in Figs. 1a and 2a. the susceptibility of χ as follows ε,N Apart from the verification of the validity of KK relations, (1) we want to provide further support for the quality of the lin- Im[χ (ω )], 0 ≤ω ≤ω , l l ω ε,N (1) ear susceptibilities considered. (1) Im[χ (ω)] = (60) Im[χ (ω)], ω ≤ω ≤ω , ε,N l h ε,N First, we test the sum rules given in Eqs. (42) and (44) for , ω ≥ω , (1) the real part of the extrapolated susceptibilities χ (ω) and ε,N (1) where the measured data are sandwiched between the low χ (ω), respectively. Our findings are presented in Fig. 3, μ,N and high frequency limit. Whereas we have a rigorous result where it is shown that an excellent agreement (within 1%) is for the high frequency limit, the low frequency limit is com- found between the theoretical values and the numerical re- (1) puted by making the reasonable assumption that the lead- sults. Since Re[χ ](ω) is negative in the high-frequency e,,a ing order of the ω → 0 limit is linear (see discussion after range, the convergence of the integral to the theoretical value (1) Eq. 11). Similarly, the real part of the susceptibility χ can of the sum rule is from above, whereas the opposite occurs ε,N (1) be redefined as follows: for Re[χ (ω)]. Extending the integral for even larger val- m,,a ues of ω would bring the numerical results to an almost per- (1) Re[χ (ω )], 0 ≤ω ≤ω , l l ε,N fect agreement with the theory. (1) (1) Re[χ (ω)] = (61) Re[χ (ω)], ω ≤ω ≤ω , l h Following the definition given in Eq. (3), the Green func- ε,N ε,N F−2m (1) − , ω ≥ω , h tion G (τ) computed for an observable 8 and a given pat- ω 8 tern of perturbation flow X (x) (in this case X = 1∀i) can i i where we have used the fact that at low frequencies the real be used to compute the time-dependent linearized impact of part of the susceptibility is constant in ω up to a quadratic all perturbations with the same spatial pattern X (x) but with term. A corresponding procedure is used to extend the spec- arbitrary time modulation. Whereas the direct estimate of the (1) tral range of the χ (ω), where the suitable asymptotic be- Green function from the time dependent dynamics can be ob- μ,N haviors described in Sect. 3.2.1 are adopted. The red lines tained by performing an ensemble of simulations where the in Figs. 1a, b–2a, b present the results of such extrapola- time modulation of the perturbation is given by a δ(t) pattern tions, and the magenta lines show the outcome of the data (see discussion in Sect. 2, we take the indirect route by con- inversion of these functions performed via KK relations. We sidering Eq. (5). By applying the inverse Fourier Transform, (1) observe that the agreement is outstanding, with almost per- we derive the Green functions corresponding to χ (ω) and ε,a fect overlap inside the region where measurement is per- (1) χ (ω). The results are presented in Fig. 4: for both observ- μ,a formed and remarkable agreement also in the low and high ables the Green functions are clearly causal, and their short- frequency range. This is a very convincing evidence that time behavior agrees remarkably well with what be deduced the Ruelle response theory can be successfully applied for by looking at the asymptotic properties of the corresponding this system. Since the KK relations provide, first and fore- susceptibilities (compare Eqs. 41 and 43). most, consistency tests, the agreement the original and the KK-transformed susceptibility automatically confirms that 4.3 Local perturbation the extrapolation procedure we have adopted is correct. A still better agreement would be found had we taken into ac- The data obtained for the numerical simulations of the re- count value of ω larger than what considered in the extrapo- sponse to the local perturbation are, given the much weaker lation used here (up to 100 π ). overall strength of the forcing, much noisier that those pre- Furthermore, let’s consider the results presented in sented in the previous section. Nevertheless, we shall see Sect. 3.1. The slopes of the functions e(F) and m(F) are that all the theoretical predictions are verified to a surpris- given by ingly good degree of approximation. The global observables E and M are of the form dm(F) γ −1 N 8(x ,...,x ) = φ(x ), where φ(x) =x /2 and φ(x) = =λγF , (62) 1 n i i=1 dF x, respectively. We have consistently verified that the iden- (1) (1) (1) tity χ = χ = χ discussed in the previous section φ,1 φ/N,N φ,N applies in the whole spectral range explored by our simula- de(F) (1+γ) (1+γ) =λ F =m(F) . (63) tions, compatibly with the (slightly) different signal-to-noise dF 2 2 www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 20 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 1.8 (1) (1) Re[χ ] ε,N ε,N (1) (1) 1.6 Re[χ ] χ E,1 μ,N (1) π Im[χ ] ε,N 1.4 (1) Im[χ ] E,1 1.2 0.8 0.6 0.4 0.2 ω ω l h −0.2 −2 −1 0 1 2 3 −2 −1 0 1 2 3 10 10 10 10 10 10 10 10 10 10 10 10 max Fig. 5. Comparison between the linear susceptibility of the intensive Figure 5: Comparison between the linear susceptibility of the intensive energy for the global Fig. 3. Sum Rules for the real part of the linear susceptibility of Figure 3: Sum Rules for the real part of the linear susceptibility of intensive energy E/N energy for the global perturbation and of the total energy for the perturbation and of the total energy for the local perturbation. The signals are the same (black line) intensi and vmomen e enertum gy EM/N /N(black (blue line) line) of and the momentum system withMresp /Nect (blue to the line) global per- except for local theperturbation. different level of The noise. signals are the same except for the different turbation with X = 1 ∀i. The theoretical values are indicated in the figure. of the system i with respect to the global perturbation with X = 1 ∀i. level of noise. The theoretical values are indicated in the figure. 0.8 Measurement a) Extrapolation (1) KKmea G (τ) 0.7 KKext Θ(τ)(m + (F − 2m)τ) (1) 2.5 Gμ (τ) Θ(τ)(1 − τ) 0.6 0.5 1.5 0.4 0.3 0.2 0.5 0.1 −0.5 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 −0.1 ωl ω −0.2 −2 −1 0 1 2 3 10 10 10 10 10 10 Fig. 4. Green functions describing the time-dependent response of Figure 4: Green functions describing the time-dependent response of the observable of in- the observable of intensive energy E/N (black line) and momentum tensive energy E/N (black line) and momentum M/N (blue line) with respect to the global 0.8 Measurement M/N (blue line) with respect to the global perturbation with X = 1 b) perturbation with X = 1 ∀i. For both observables, the short time behaviori of the Green Extrapolation KKmea 0.7 KK function∀iestimated . For both from observ the asymptotic ables, the b short ehavior time of the beha corresp vior of onding the Green susceptibilit func-y is shown ext in figure. tion estimated from the asymptotic behavior of the corresponding 0.6 susceptibility is shown in figure. 0.5 0.4 ratios in the two sets of simulations. See, e.g., Fig. 5 for the (1) (1) 0.3 comparison between the two susceptibilities χ and χ . E,1 ε,N 0.2 We then proceed to analyze more in detail the linear susceptibilities related to local observables. In Fig. 6 we 0.1 present our results concerning the real and imaginary part (1) of χ (ω). Analogously to what observed in the previous E ,1 −0.1 40 ω h subsection, we have that once the measured susceptibility l −0.2 −2 −1 0 1 2 3 is extrapolated using the theoretical results obtained via re- 10 10 10 10 10 10 sponse theory and KK relations, we have an excellent agree- ment between the original real and imaginary parts and those Fig. 6. Linear susceptibility of the energy at the grid point where the obtained using the KK inversion. The KK algorithm, instead, local perturbation is applied. The real and the imaginary parts are provides only a partially satisfying outcome when only data depicted in (a) and (b), respectively. The measured and extrapolated from the measured range are considered. Relatively discrep- values are shown in red and black lines, respectively. The result of the Kramers-Kronig inversion done with the measured and with ancies are found near the boundaries of the data range, with with the extrapolated data are shown in blue and magenta lines, an especially serious bias near ω for the imaginary part of respectively. the susceptibility. Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ max dωRe[χ] G(τ) 0 Re[χ], Im[χ] Im[χ] Re[χ] V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 21 (1) (1) 0.25 When comparing χ (ω) and χ (ω) (see Fig. 5), we Measured E ,1 E,1 a) j Extrapolation KK mea KK observe that for low frequency the response of the energy ext 0.2 at the grid point where the perturbation is applied accounts for about half of the response of the total energy, thus im- 0.15 plying that the remaining half is redistributed among the re- maining N − 1 grid points. The relevance of the response of 0.1 grid points other than the directly perturbed one also explains (1) (1) 0.05 why the peak of Imχ (ω) (and so of Imχ (ω)) than that E,1 ε,N (1) of Imχ (ω) – see the frequency range 2. ω . 4. Slower E ,1 perturbations allow other grid points x 6= x to respond ef- k j −0.05 fectively. (1) Instead, since the leading asymptotic order of χ (ω) E ,1 −0.1 −2 −1 0 1 2 3 10 10 10 10 10 10 (1) and χ (ω) is the same, at high frequencies the local en- E,1 ergy response accounts for most of the energy response of 0.25 Measured b) Extrapolated KK the whole system. In this case, the incoming perturbation is mea KK ext 0.2 so fast that the internal time scales of the system as bypassed, and mainly a local effect is observed. Nevertheless, the sec- 0.15 ond leading order of the asymptotic expansion of the two susceptibilities has opposite sign (see Eqs. 41 and 54), which 0.1 suggests that at any large but finite frequency the local energy response is only a good approximation to the response of the 0.05 total energy. The changeover between the two regimes oc- (1) curs around the frequency of the peak of Imχ (ω), which E ,1 corresponds to a perturbation with period close to 1. −0.05 (1) Thanks to the asymptotic equivalence between χ (ω) E ,1 (1) (1) −0.1 −2 −1 0 1 2 3 and χ (ω) (and χ (ω)), they must obey the same sum 10 10 10 10 10 10 E,1 ε,N for the real part of the susceptibility (see Eqs. 42–55), even if the two real parts, as discussed above, are rather different in Fig. 7. Linear susceptibility of the momentum at the grid point value in the low frequency range and even in sign in the high where the local perturbation is applied. The real and the imaginary frequency range. Figure 8 confirms that this rather counter- parts are depicted in (a) and (b), respectively. The measured and extrapolated values are shown in red and black lines, respectively. intuitive behavior is actually observed. Note also that sum The result of the Kramers-Kronig inversion done with the measured rules, resulting from an integration, are less sensitive to noise and with with the extrapolated data are shown in blue and magenta in the data, but this occurs if and only if the underlying sig- (1) lines, respectively. nal is correct. Therefore, we understand that in χ (ω) and E,1 (1) χ (ω) the strong static and quasi-static response and the ε,N The investigation of the linear susceptibility of the vari- (rather odd) negative sign for high frequencies of the real part able x is not as insightful as that of E . We find that lin- j j of the linear susceptibility, which are crucially related to the (1) (1) ear susceptibility χ (ω) is quite similar to χ (ω) (and behavior for the grid points different from the perturbed one) x ,1 M,1 (1) compensate each other to guarantee agreement with the sum χ (ω), see Fig. 2) in both the real and imaginary parts at μ,N (1) rule obtained from the real part of χ (ω), which instead all frequency. The only notable differences are that the static E ,1 has a smaller range and more regular (monotonic) behavior response x is slightly larger than than of M , and that the with frequency. imaginary features a secondary peak at slightly larger fre- A formally similar – and analogously spectacular – spec- quencies than the main spectral feature. We have verified,as tral compensation has been observed in a physical process as in the previous cases, the results of the numerical simula- different from what we are analyzing here as the electromag- tions accurately agree with the theoretical results regarding netically induced transparency (Cataliotti et al., 1997). The the asymptotic behavior of both the real and imaginary part result obtained here supports previous findings obtained on and that KK relations map to high degree of precision the quasi-equilibrium systems suggesting that sum rules do not real and the imaginary parts into each other. See Fig. 7 for depend on many-particle interactions (Lucarini et al., 2003, details. 2005). We present as main finding of the analysis of the observ- able x that, as predicted by the theory, the real part of (1) (1) χ (ω) obeys the same sum rule as the real part of χ (ω) x ,1 M,1 www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 Im[χ] Re[χ] 22 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings The first scenario envisioned here pertains to the pair of (1) Ej,1 (1) (1) (1) (1) χ , χ ε,N E,1 linear susceptibilities χ (ω) and χ (ω), whereas the (1) π x ,1 M,1 m j x ,1 2 (1) (1) χ , χ second scenario is related to the pair of linear susceptibili- μ,N M,1 (1) (1) ties χ (ω) and χ (ω). Note that, taking into account the E ,1 E,1 asymptotic properties of the susceptibility of the observable 2 2 2 2 E = 1/2x + 1/2x + 1/2x + 1/2x , discussed in loc j j+1 j+2 j−1 (1) (1) the previous section we conclude that χ and χ should E ,1 E,1 loc be similar for all values of ω. Obviously, a similar argument applies if the leading order is real. This discussion further clarifies that the higher the −2 −1 0 1 2 3 10 10 10 10 10 10 max number of independent KK relations and related sum rules verified by a susceptibility functions, the more stringent are Fig. 8. Sum rules of the real part of the linear susceptibilities in- the constraints on its properties. Figure 8: Sum rules of the real part of the linear susceptibilities indicated in the legend. The dicated in the legend. The theoretical values are indicated in the theoretical values are indicated in the figure. 4.5 Additional properties of the linear susceptibility figure. The special mathematical properties of the linear susceptibil- (1) or of χ (ω), because the corresponding imaginary parts ities allow to investigate further properties of the response. μ,N (1) feature the same asymptotic behavior. Figure 8 shows that in In particular, we note that for m ≥ 1 the function [χ ] is (1) the case of the momentum variables the cumulative integral analytic in the upper complex ω-plane just as as χ . This is rather similar for the susceptibility of the local and of the allows, as discussed in Lucarini et al. (2005) to derive, in global variable, with small discrepancies in the region around principle, an infinite set of integral relations (KK and sum the peak of the response. rules) deriving just from the holomorphic proprieties of the susceptibility. As an example, we have considered the square 4.4 Further implications of Kramers-Kronig relations (1) of the linear susceptibility [χ (ω)] . From Eq. (41), it is ε,N and sum rules easy to prove that the following asymptotic expansion holds for large values of ω: We now show how the knowledge of the asymptotic behav- ior of the real and imaginary part and the knowledge of the 2 m m(F − 2m) (1) 2 −4 [χ (ω)] = − + +o(ω ). (64) validity of the KK relations and related sum rules allow to ε,N 2 3 ω ω draw general conclusions on the similarities and differences As shown in panel (a) of Fig. 9, KK relations are found between two given linear susceptibility functions. Let’s con- to connect up to a high degree of approximation the real sider the case that these two susceptibilities feature the same (1) and imaginary part of [χ (ω)] . Moreover, thanks to the first order asymptotic expansion in the high frequency limit. ε,N asymptotic behavior given in Eq. (64), it is possible to estab- Let’s assume that it is an odd power of ω, so that the real part lish, thanks to Eqs. (13–14), the following sum rules: is negligible for high frequencies. Therefore, the two suscep- th tibilities will obey the same sum rule for, e.g. the 0 moment (1) dνRe[χ (ν)] = 0, (65) of the real part. ε,N If they agree also in the asymptotic behavior of the real part, they cannot feature large discrepancies in the low fre- (1) 2 2 quency range of the real part of the susceptibility either, or dννIm[χ (ν)] = m , (66) ε,N otherwise the agreement of the sum rules would be broken. Therefore, the real part of the two susceptibilities are similar, Panel (b) of Fig. 9 shows that the obtained numerical results and, as a consequence of the KK relations, the two imaginary are in excellent agreement with the theoretical predictions. parts will also be similar. Note that these results do not have an obvious physical in- (1) If, instead, there is a discrepancy in the asymptotic behav- 2 terpretation, as the inverse Fourier Transform of [χ (ω)] ε,N ior of the real part of the two susceptibilities, the two real is given by the convolution product of the Green function parts will necessarily be rather different in the low frequency (1) G (t) with itself, while they depend only on the formal ε,N range, again in order to comply with the sum rule constraint. properties of the linear susceptibility. As the two real parts are different, the imaginary part of the two susceptibility will also be rather different, except, from hypothesis, in the high-frequency range. 5 Practical implications for climate change studies In this paper we have constructed and verified to a high degree of accuracy the linear response theory for a simple Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ max dωRe[χ] 0 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 23 From Eq. (2), we have that (1) 2 a) Re[(χ ) ] ε,N (1) Im[(χ ) ] ε,N (1) 2 Z 2.5 Re[(χ ) ]KK ε,N +∞ (1) (1) Im[(χ ) ]KK (1) ε,N hT i (t) = dσ G (σ )f(t −σ ). (67) S 1 1 1 −∞ 1.5 In practical terms, the left hand side of this equation is noth- ing but the ensemble average of the time series of the change between the globally averaged surface temperature of the planet at a time t after the perturbation has started. Note that 0.5 (1) the direct estimate of G (σ) is likely to be overwhelmingly difficult. Using Eq. (4), we have that: −0.5 (1) (1) hT i (ω) =χ (ω)f(ω), (68) −1 −2 −1 0 1 2 3 10 10 10 10 10 10 which implies that once we compute the Fourier Transform of the time series mentioned above and we know the modu- SR0 b) SR lating function f(t) (and so its Fourier Transform f(ω)), we (1) can reconstruct χ (ω). Let’s select a particularly simple π 2 example of modulating function f(t) = (2(t)−2(t −τ)). This is just a rectangular function of width τ , of height , and shifted from the origin by a forward time translation τ/2. In practical terms, this corresponds to changing abruptly the field CO concentration by  at time t = 0 and taking it back to its original value at t =τ . we then obtain: (1) (1) hT i (ω) hT i (ω) S S (1) χ (ω) = =ω . (69) f(ω) (sin(ωτ)+i(1− cos(ωτ))) (1) Once we know χ (ω), as widely discussed in this paper, 0 T −2 −1 0 1 2 3 S 10 10 10 10 10 10 ω (1) max we can compute G (t), and we know everything about the (1) Fig. 9. Properties of the square of the linear susceptibility χ . response of the system at all time scales, including the static ε,N (1) 2 response. Note that any choice of f(t) is equally valid to The real and imaginary parts of [χ ] with their KK transforms ε,N set up this procedure as long as f(t) is square integrable. are depicted in the panel (a), the vanishing sum rule for the real part This implies that, in a very profound way, the kind of forcing and the non-vanishing sum rule for the imaginary part are depicted in panel (b). scenarios used in the various assessment Reports of IPCC, where the CO concentration typically stabilizes at a differ- ent value from the preindustrial one (so that f(t) does not yet prototypical climate model by computing the frequency- tend to 0 as t → ∞) are not necessarily the only nor the best dependent susceptibilities of several relevant observables re- ones, in spite of what could be intuitively guessed, to study lated to localized and global patterns of forcings. These re- even the steady state response of the system. sults pave the way for devising a rigorous methodology to Obviously, a similar set of experiments could be devised be used by climate models of any degree of complexity for for studying rather thoroughly the response of the climate studying climate change at, in principle, all time scales us- system to a variety of forcings, such as changes in the O ing only a very limited set of experiments, and for exploiting concentration, aerosols, solar radiance, as well as to changes effectively the currently adopted ensemble runs methods. in the parameterizations. In the case of uncoupled models of Let’s consider, for sake of simplicity, that the observable one subdomain of the climate system (e.g. atmospheric and 8 is the time-dependent globally averaged surface temper- oceanic GCMs, land-surface models), this strategy could be ature of the planet T , that F(x) represents the whole set used to study the impact of perturbations to the boundary of climate equations in a baseline scenario (e.g., with pre- conditions provided by the other subdomains of the climate industrial CO concentration), and that the perturbation field system. X(x) is nothing but a constant field of CO concentration, which directly impacts only the radiative part of the code. The perturbation is modulated by a time-dependent function f(t) to be specified below. We assume, for simplicity, that the model does not feature daily or seasonal variations in the radiative input at the top of the atmosphere. www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 24 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 6 Summary, discussion and conclusions stochastic dynamics using the results of the mathematical theory of random dynamical systems is emerging as new, The climate can be seen as a complex, non-equilibrium promising paradigm for the investigation of the structural system, which generates entropy by irreversible processes, properties of the climate system (Chekroum et al., 2011). transforms moist static energy into mechanical energy We have proposed a different perspective. In agreement (Lorenz, 1967; Peixoto and Oort, 1992) as if it were a heat with the view given above, we have taken as mathematical engine (Johnson, 2000; Lucarini, 2009a), and, when the ex- framework for the analysis of the climate system that of non- ternal and internal parameters have fixed values, achieves equilibrium statistical mechanics, and have focused on the a steady state by balancing the input and output of energy steady state properties of ergodic dynamical systems (Eck- and entropy with the surrounding environment (Ozawa et al., mann and Ruelle, 1985) possessing the special property of 2003; Lucarini, 2009a). For such basic reasons, the tool of having an invariant measure of the SRB type (Ruelle, 1989). equilibrium and quasi-equilibrium statistical mechanics can- As proposed by the chaotic hypothesis (Gallavotti, 1996; not provide suitable tools for studying the fundamental prop- Gallavotti and Cohen, 1995), this mathematical framework erties of the climate system. In particular, the fluctuation- is well suited for analyzing general non-equilibrium physical dissipation theorem, which allows for deriving the properties systems. of the response of the system to external perturbations from In this context, the impact on the system of general per- the observations of its internal variability cannot be applied. turbation can be treated using the response theory recently It is reasonable to ask whether is possible to evaluate how introduced by Ruelle (1998a,b, 2009), which allows to com- far from equilibrium the climate system actually is. It is pos- pute the change in the expectation value of a generic observ- sible to evaluate such distance in a mathematically sound able as a perturbative series where each term is given by the way by assessing the ratio of the dimensionality of the at- average over the unperturbed invariant measure of a function tractor of the system over the total number of degrees of of the phase space which depends on the considered observ- freedom. Whereas a ratio close to one indicates that only able and on the applied perturbation. In other terms, even if small deviations from equilibrium are present, a small ratio the internal dynamics of the system is nonlinear and chaotic, suggests that strongly non-equilibrium conditions are estab- the leading order of the response is in general linear with the lished. See Posch and Hoover (1998) for a detailed treatment strength of the added perturbation. This approach overcomes of this problem in the classical case of heat conduction. In the difficulties related to the singularity of the invariant mea- the case of a quasi-geostrophic atmospheric model forced by sure discussed in Thuburn (2005). Earth-like boundary conditions, the dimensionality of the at- At each order, the propagator of the perturbation, i.e. the tractor of the model is about one order of magnitude smaller Green function, is causal. This allows for applying disper- than the total number of degrees of freedom (Vannitsem and sion theory and establish general integral constraints – KK Nicolis, 1997). While not conclusive, this seems to suggest relations – connecting the real and imaginary parts of the sus- that the best framework to interpret the climate is that of a far ceptibility, i.e. the Fourier Transform of the Green function from equilibrium system. (Lucarini, 2008b, 2009b) Following either explicitly or implicitly the programme of In this paper we have first recapitulated the main as- the Catastrophe theory (Arnold, 1992), many authors have pects of the general response theory and have propose some approached the problem of understanding the fundamental new general results, which boil down to consistency rela- properties of the climate system by looking at the detailed tions between the linear susceptibilities of different observ- structure of the bifurcations of the deterministic dynamical ables. The obtained equation provides the basic idea of why system constructed heuristically in order to represent the dy- the fluctuation-dissipation theorem does not apply in non- namics of the main climate modes using as few degrees of equilibrium cases. freedom as possible. Such an approach often hardly allows to efficiently represent the fluctuations and the statistical prop- We have showed for the first time that the Ruelle linear erties of the system. The introduction of stochastic forcing response theory can be applied with great success to ana- provides a relatively simple but conceptually rich partial so- lyze the climatic response to general forcings. We have cho- lution to some of these draw-backs, even if the Hasselmann sen as test bed the L96 model (Lorenz, 1996; Lorenz and programme (Hasselmann, 1976) suffers from the need for a Emanuel, 1998; Lorenz, 2004), which, in spite of its simplic- – usually beyond reach – closure theory for the properties ity, has a well-recognized prototypical value as it is a spa- of noise. Therefore, the stochastic component is usually in- tially extended one-dimensional model and presents the ba- troduced ad hoc, with the ensuing lack of universality and/or sic ingredients, such as dissipation, advection and the pres- robustness when various levels of truncations are considered. ence of an external forcing, of the actual atmosphere. Such These strategies have anyway brought to outstanding scien- a model features a different level of complexity with respect tific results and has been suggested the existence of generic to those adopted in previous numerical investigations of Ru- mathematical structures present in hierarchies of CMs (Saltz- elle’s theory (Reick, 2002; Cessac and Sepulchre, 2007; Lu- man, 2002). Recently, the unified treatment of chaotic and carini, 2009b), and has been already tested in terms of linear Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 25 response (albeit with a different methodology and in a differ- We believe that the proposed approach, which we may ent theoretical context) in Abramov and Majda (2007). dub as spectroscopy of the climate system, may constitute a We have analyzed the frequency dependence of the re- mathematically rigorous and practically very effective way sponse of the local and global energy and momentum of the to approach the problem of evaluating climate sensitivity system to perturbations having a global spatial pattern and to and climate change from a radically new perspective. In perturbations acting only on one grid point. We have derived this regard, we have proposed a rigorous way to compute, analytically several properties of the corresponding suscep- e.g., the surface temperature response to changes in the the tibilities, such as asymptotic behavior, validity of KK rela- CO concentration at all time scales using only a specific tions, and sum rules. We have shown that all the coefficients set of simulations, and taking advantage of the theoretical of the leading asymptotic expansions as well as the integral results presented here. Given the ever-increasing interest constraints can be written as linear functions of parameters towards decadal and multidecadal climate prediction, these that describe unperturbed properties of the system, and in tools could be of relevant practical interest and their applica- particular its average energy and average momentum. The bility could benefit from technological platform aimed at cre- theory has been used to explain differences in the response of ating ensemble simulations comprising of many members. local and global observables, in defining the intensive proper- We underline that our approach takes into account all the ties of the system and in generalizing the concept of climate (linear and nonlinear) feedbacks of the system, as they are sensitivity to all time scales. included in the definition of the Green function. This, at a We have then verified the theoretical predictions from the very practical level, is the great advantage of using Ruelle’s outputs of the simulations up to a high degree of precision, formulas. even if we have used rather modest computational resources At a more basic level, whereas considering more complex (a total of about 30 cpu days of a mid-range commercial lap- models requires heavier computational resources, the modest top). We have verified that the linear response theory holds cost of the present set of simulations suggests that, at least for for perturbations of intensity accounting to up to about 10% global or regional climatic observables, it is feasible to test of the unperturbed forcing terms. Even when local pertur- the theory discussed here for simplified yet Earth-like climate bation and local observables are considered it is possible models without resorting to top-notch computing facilities. to achieve a signal-to-noise ratio which permits rather sat- Moreover, while in this paper we have computed the sus- isfactory comparisons with the theory. We have proved that ceptibilities using, on purpose, a very cumbersome method, the combined use of KK relations and the knowledge of the more efficient strategies can be devised, at least when the asymptotic behavior of the susceptibilities allows for extrap- linear regime of the response is considered. Apart from the olating in a rigorous way the observed data. We also have practical example given for the case of the impact of the CO shown how to reconstruct the linear Green function, which concentration, these include studying the response of the sys- can be used to map perturbations of general time modula- tem to δ(t) like perturbations, which gives directly the Green tion into changes in the expectation value of the considered function of the system, and including in the forcing various observable for finite as well as infinite time. monochromatic signals. Of course, in all cases, a Monte Our numerical experiments have been performed using Carlo approach is needed in order to sample effectively the one of the standard settings of the L96 model, namely the attractor of the unperturbed system in terms of the initial con- version identified by having N = 40 degrees of freedom and ditions of the simulations. forcing F = 8. Nevertheless, some newly obtained empirical These results pave the way for future investigations aimed closure equations expressing the average energy and the av- at improving and extending the theoretical framework pre- erage momentum of the unperturbed system as simple power sented here, at finding results of general applicability in the laws of F (with no dependence on N ) have allowed to ex- context of the modelling of geophysical fluid dynamics, and, tend our results to the entire class of chaotic L96 models. finally at answering specific questions of relevance for cli- The regular scaling of the properties of the system and of its mate dynamics. In this paper we have analyzed the sim- response with F agrees with what observed in Abramov and ple case of the linear response, but, as discussed in Ruelle Majda (2007). (1998b) and Lucarini (2008b, 2009b), we have the algorithm In this paper we have only used the KK relations in the to compute higher order terms, so that the treatment of the most simplistic framework, i.e., computing the KK trans- nonlinear response in entirely feasible. forms and evaluating their agreement with the original data. In the first direction, we foresee the possibility of writing Actually, several more sophisticated analysis techniques are out explicitly the linear susceptibility of a general observable available, such as recursive self-consistent algorithms, where by projecting the perturbation onto the unstable, neutral and the measured data are taken as first guess, exploiting the fact stable manifolds and analyzing separately the contributions that multiple applications of KK relations, combined with to the total response. This will probably require the adoption sum rules, automatically filter our the noise and remove most of adjoint techniques, and will benefit from the recent algo- of the spurious signal (Lucarini et al., 2005). rithms proposed in Abramov and Majda (2007) and Majda and Wang (2010). Moreover, we will be testing the radius of www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 26 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings convergence of the Ruelle response theory in some specific tives such as climateprediction.net, and the new project examples. PCMDI/CMIP5 (http://cmip-pcmdi.llnl.gov/cmip5/), which Along the second direction, we propose to study the im- will provide a crucial input for the Fifth Assessment Report pact of stochastic forcing to deterministic chaotic models by of the IPCC. treating the (additive or multiplicative) noise as a perturba- Acknowledgements. VL acknowledges the financial support of the tion to be analyzed using the linear and nonlinear Ruelle EU-ERC research grant NAMASTE and of the Walker Institute response theory and related spectral methods. We will try Research Development Fund. The authors acknowledge the fruitful to complement the results recently obtained in Majda and interactions with the participants to the programme “Mathematical Wang (2010) from an entirely different approach, i.e. bypass- and Statistical Approaches to Climate Modelling and Prediction” ing the Fokker-Planck equation formalism and using directly held at the Newton Institute for Mathematical Sciences (Cam- the first and higher order terms of the Ruelle response for- bridge, UK). mula. Moreover, we shall look into the spectral peaks of the susceptibilities and try to understand how the amplification Edited by: H. A. Dijkstra Reviewed by: two anonymous referees of the response is related to resonances of the system and to the activation of positive feedbacks. Since at tipping points the climate response to perturbations is expected to diverge, References the formalism discussed here could be used also for under- standing the processes leading to large-scale changeovers of Abramov, R. and Majda, A. J.: Blended response algorithms for the dynamics of climate. linear fluctuation-dissipation for complex nonlinear dynamical Along the third direction, we envision the analysis of the systems, Nonlinearity, 20, 2793–2821, 2007. impact of topography on the statistical properties of the cir- Arnold, V. I.: Catastrophe Theory, Springer, Berlin, 1992. culation in a quasi-geostrophic setting, thus extending in a Baladi, V.: On the susceptibility function of piecewise expanding interval maps, Commun. Math. Phys., 275, 839–859, 2007. climatic perspective what presented in Speranza et al. (1985). Brayshaw, D. J., Hoskins, B., and Blackburn, M.: The storm track Moreover, we will tackle in an idealized setting the problem response to idealised SST perturbations in an aquaplanet GCM, of computing the response of the storm track to changes in J. Atmos. Sci., 65, 2842–2860, 2008. the surface temperature (Brayshaw et al., 2008). Positive re- Cacuci, D. G.: Sensitivity theory for nonliner systems. Part I. Non- sults in this direction, albeit using a different formalism, have linear functional analysis approach, J. Math. Phys., 22, 2794– been shown in Gritsun et al. (2008). Moreover, we will try to 2802, 1981a. compute along the way discussed here the climate response Cacuci, D. G.: Sensitivity theory for nonliner systems. Part II. Ex- to changes in CO and solar irradiance using simplified but tensions to additional classes of responses, J. Math. Phys., 22, rather valuable climate models like PLASIM (Fraedrich et 2803–2812, 1981b. al., 2005). Cataliotti, F. S., Fort, C., Hansch, ¨ T. W., Inguscio, M., and The results of the linear response theory should be com- Prevedelli, M.: Electromagnetically induced transparency in cold free atoms: Test of a sum rule for nonlinear optics, Phys. Rev. A, pared to those obtained from a simplified application of the 56, 2221–2224, 1997. fluctuation-dissipation theorem. The difference between the Cessac, B. and Sepulchre, J.-A.: Linear response, susceptibility and two results would inform us on the relevance of the ef- resonances in chaotic toy models, Physica D, 225, 13–28, 2007. fect of the external perturbation along the stable direction Chekroun, M. D., Simonnet, E., and Ghil, M.: Stochastic climate of the unperturbed flow, which cannot be precisely captured dynamics: Random attractors and time-dependent invariant mea- by the internal fluctuations of the system and, thus, can- sures, Physica D, in press, 2011. not be rigorously described in the context of the fluctuation- Dolgopyat, D.: On differentiability of SRB states for partially hy- dissipation theorem. Moreover, this could also suggest ways perbolic systems, Invent. Math., 155, 389–449, 2004. to improve the empirical evaluation and parameterisation Eckmann, J.-P. and Ruelle, D.: Ergodic theory of choas and strenge of the corresponding contribution to the response. This attractors, Rev. Mod. Phys., 57, 617–655, 1985. would also be of relevance to understand why, in spite of Errico, R.: What Is an Adjoint Model?, B. Am. Meteor. Soc., 78, 2577–2591, 1997. all the involved approximations and the theoretically unjusti- Eyink, G. L., Haine, T. W. N., and Lea, D. J.: Ruelle’s linear re- fiable assumptions, previous applications of the equilibrium sponse formula, ensemble adjoint schemes and Lvy flights, Non- fluctuation-dissipation theorem (often in a hyper-simplified linearity, 17, 1867–1889, 2004, quasi-gaussian setting) have been relatively successful in Fertig, E., Harlim, J., and Hunt, B.: A comparative study of 4D-Var predicting the climate response to external forcings. and 4D ensemble Kalman filter: perfect model simulations with Finally, we would like to remark that the theory Lorenz-96, Tellus A, 59, 96–101, 2007. and the practical recipes proposed here could be of di- Fraedrich, K., Jansen, H., Kirk, E., Luksch, U., and Lunkeit, F.: The rect interest for all projects aimed at auditing climate Planet Simulator: Towards a user friendly model, Meteorol. Z., models’ performances and at studying practical problems 14, 299–304, 2005. connected to climate change, such as PCMDI/CMIP3 Frye, G. and Warnock, R. L.: Analysis of partial wave dispersion (http://www-pcmdi.llnl.gov/), distributed computing initia- relations, Phys. Rev., 130, 478–494, 1963. Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 27 Gallavotti, G.: Chaotic hypotesis: Onsanger reciprocity and Leith, C. E.: Climate response and fluctuation dissipation, J. Atmos. fluctuation-dissipation theorem, J. Stat. Phys., 84, 899–926, Sci., 32, 2022–2026, 1975. 1996. Lorenz, E. N.: Deterministic nonperiodic flow, J. Atmos. Sci., 20, Gallavotti, G.: Nonequilibrium statistical mechanics (stationary): 130–141, 1963. overview, in: Encyclopedia of Mathematical Physics, edited by: Lorenz, E. N.: The Nature and Theory of the General Circulation of Franc ¸oise, J.-P., Naber, G. L., Tsou Sheung Tsun, Elsevier, Am- the Atmosphere, WMO, Geneva, 1967. sterdam, 530–539, 2006. Lorenz, E.: Forced and free variations of weather and climate, J. Gallavotti, G. and Cohen, E. D. G.: Dynamical ensembles in sta- Atmos. Sci., 36, 1367–1376, 1979. tionary states, J. Stat. Phys., 80, 931–970, 1995. Lorenz, E. N.: Predictability-A problem partly solved, Proc. Sem- Giering, R. and Kaminski, T.: Recipes for Adjoint Code Construc- inar on Predictability, Vol. 1, Reading, Berkshire, United King- tion, ACM Trans. On Math. Software, 24, 437–474, 1998. dom, ECMWF, 1–18, 1996. Ghil, M. and Malanotte-Rizzoli, P.: Data assimilation in meteorol- Lorenz, E. N.: Designing Chaotic Models, J. Atmos. Sci., 62, 1574– ogy and oceanography, Adv. Geophys., 33, 141–266, 1991. 1587, 2004. Gritsun, A. and Branstator G.: Climate response using a three- Lorenz, E. N. and Emanuel, K.: Optimal sites for supplementary dimensional operator based on the fluctuation-dissipation theo- weather observations, J. Atmos. Sci., 55, 399–414, 1998. rem, J. Atmos. Sci., 64, 2558–2575, 2007. Lucarini, V.: Validation of climate models, in: Encyclopedia of Gritsun A., Branstator G. and Majda A. J.: Climate response of Global Warming and Climate Change, edited by: Philander, G., linear and quadratic functionals using the fluctuation-dissipation SAGE, Thousand Oak, 1053–1057, 2008a. theorem, J. Atmos. Sci., 65, 2824–2841, 2008. Lucarini, V.: Response theory for equilibrium and non-equilibium Hall, M. C. G.: Application of adjoint sesnitivity theory to an atmo- statistical mechanics: causality and generalized Kramers- spheric general circulation model, J. Atmos. Sci., 42, 2644–2651, Kroning relations, J. Stat. Phys., 131, 543–558, 2008b. 1986. Lucarini, V.: Thermodynamic efficiency and entropy produc- Hall, M. C. G. and Cacuci, D. G.: Physical interpretation of the tion in the climate system, Phys. Rev. E, 80, 021118, adjoint functions for sensitivity analysis of atmospheric models, doi:10.1103/PhysRevE.80.021118, 2009a. J. Atmos. Sci., 40, 2537–2546, 1983. Lucarini, V.: Evidence of Dispersion Relatios for the Nonlinear Re- Hall, M. C. G., Cacuci, D. G., and Schlesinger, M. E.: Sensitivity sponse of Lorenz 63 System, J. Stat. Phys., 134, 381–400, 2009b. analysis of a radiative-convective model by the adjoint method, Lucarini, V., Bassani, F., Saarinen, J. J., and Peiponen, K.-E.: Dis- J. Atmos. Sci., 39, 2038–2050, 1982. persion theory and sum rules in linear and nonlinear optics, Riv- Hallerberg, S., Pazo, ` D., Lopez, ` J. M., and Rodr´ ıguez, ista Nuovo Cimento, 26, 1–120, 2003. M. A.: Logarithmic bred vectors in spatiotemporal Lucarini, V., Saarinen, J. J., Peiponen, K.-E., and Vartiaine, E. chaos: Structure and growth, Phys. Rev. E, 81, 066204, M.: Kramers-Kronig relations in Optical Materials Research, doi:10.1103/PhysRevE.81.066204, 2010. Springer, Heidelberg, 2005. Hasselmann, K.: Stochastic climate models, Part 1: Theory, Tellus, Lucarini, V., Speranza, A., and Vitolo, R.: Parametric smoothness 28, 473–485, 1976. and self-scaling of the statistical properties of a minimal climate Haven, K., Majda, A., and Abramov, R.: Quantifying predictability model: What beyond the mean field theories?, Physica D, 234, through information theory: small sample estimation in a non- 105–123, 2007. Gaussian framework, J. Comp. Phys., 206, 334–362, 2006. Marini Bettolo Marconi, U., Puglisi, A., Rondoni, L., and Vulpi- Held, I.: The gap between simulation and understanding in climate ani, A.: Fluctuation-dissipation: Response theory in statistical modeling, B. Am. Meteor. Soc., 86, 1609–1614, 2005. physics, Phys. Rep., 461, 111–195, 2008. Intergovernmental Panel on Climate Change 2007: The Physical Majda, A. J. and Wang, X.: Linear response theory for statistical en- Science Basis. Contribution of Working Group I to the Fourth sembles in complex systems with time-periodic forcing, Comm. Assessment Report of the Intergovernmental Panel on Climate Math. Sci., 8, 145–172, 2010. Change, Cambridge University Press, Cambridge, 2007. Nussenzveig, H. M.: Causality and Dispersion Relations, Academic Johnson, D. R.: Entropy, the Lorenz Energy Cycle and Climate, Press, New York, 1972. in: General Circulation Model Development: Past, Present and Orrell, D.: Model Error and predictability over Different Timescales Future, edited by: Randall, D. A., Academic Press, New York, in the Lorenz-96 System, J. Atmos. Sci., 60, 2219–2228, 2003. 659–720, 2000. Ozawa, H., Ohmura A., Lorenz R. and Pujol T.: The second law Kubo, R.: Statistical-mechanical theory of irreversible processes. I, of thermodynamics and the global climate system: a review of J. Phys. Soc. Jpn., 12, 570–586, 1957. the maximum entropy production principle, Rev. Geophys., 41, Kubo, R.: The fluctuation-dissipation theorem, Rep. Prog. Phys., 1018, doi:10.1029/2002RG000113, 2003. 29, 255–284, 1966. Peiponen, K.-E., Vartiainen, E. M., and Asakura, T.: Dispersion, Lacorata, G. and Vulpiani, A.: Fluctuation-Response Relation and Complex Analysis and Optical Spectroscopy, Springer, Heidel- modeling in systems with fast and slow dynamics, Nonlin. Pro- berg, 1999. cesses Geophys., 14, 681–694, doi:10.5194/npg-14-681-2007, Peixoto, J. and Oort, A.: Physics of Climate, Springer, New York, 2007. 1992. Lange, P. L. and Alexeev, V. A.: Estimating 2 × CO warming Penland, C.: Noise out of chaos and why it won’t go away, B. Am. in an aquaplanet GCM using the fluctuation-dissipation theorem, Meteor. Soc., 84, 921–925, 2003. Geophys. Res. Lett., 32, L23708, doi:10.1029/2005GL024136, Posch, H. A. and Hoover, W. G.: Heat conduction in one- 2005. dimensional chains and nonequilibrium Lyapunov spectrum, www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 28 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings Phys. Rev. E, 58, 4344–4350, 1998. Shukla, J., Hagedorn, R., Hoskins, B., Kinter, J., Marotzke, J., Reick, C. H.: Linear response of the Lorenz system, Phys. Rev. E, Miller, M., Palmer, T., and Slingo J.: Revolution in climate pre- 66, 036103, doi:10.1103/PhysRevE.66.036103, 2002. diction is both necessary and possible: A declaration at the world Ring, M. J. and Plumb, R. A.: The response of a simplified modelling summit for climate prediction, B. Am. Meteor. Soc., GCM to axisymmetric forcings: Applicability of the fluctuation- 90, 175–178, 2009. dissipation theorem, J. Atmos. Sci., 65, 3880–3898, 2008. Thuburn, J.: Climate sensitivities via a Fokker-Planck adjoint ap- Ruelle, D.: Chaotic Evolution and Strange Attractors, Cambridge proach, Q. J. Roy. Meteor. Soc. 131, 73–92, 2005. University Press, Cambridge, 1989. Trevisan, A. and Uboldi, F.: Assimilation of standard and targeted Ruelle, D.: General linear response formula in statistical mechanics, observations within the unstable subspace of the observation- and fluctuation-dissipation theorem far from equilibrium, Phys. analysis-forecast cycle, J. Atmos. Sci., 61, 103–113, 2004. Lett. A, 245, 220–224, 1998a. Trevisan, A., D’Isidoro, M., and Talagrand, O.: Four-dimensional Ruelle, D.: Nonequilibrium statistical mechanics near equilibrium: variational assimilation in the unstable subspace and the opti- computing higher order terms, Nonlinearity, 11, 5–18, 1998b. mal subspace dimension, Q. J. Roy. Meteor. Soc., 136, 487–496, Ruelle, D.: A review of linear response theory for general differen- 2010. tiable dynamical systems, Nonlinearity, 22, 855–870, 2009. Vannitsem, S. and Nicolis, C.: Lyapunov Vectors and Error Growth Saltzman, B.: Dynamic Paleoclimatology, Academic Press, New Patterns in a T21L3 Quasigeostrophic Model, J. Atmos. Sci., 54, York, 2002. 347–361, 1997. Speranza, A., Buzzi, A., Trevisan, A., and Malguzzi, P.: A The- Wilks, D. S.: Effects of stochastic parametrizations in the Lorenz ory of Deep Cyclogenesis in the Lee of the Alps. Part I: Mod- ’96 system, Q. J. Roy. Meteor. Soc., 131, 389–407, 2006. ifications of Baroclinic Instability by Localized Topography, J. Atmos. Sci, 42, 1521–1535, 1985. Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Nonlinear Processes in Geophysics Unpaywall

A statistical mechanical approach for the computation of the climatic response to general forcings

Nonlinear Processes in GeophysicsJan 12, 2011

Loading next page...
 
/lp/unpaywall/a-statistical-mechanical-approach-for-the-computation-of-the-climatic-9JrD7pNOFb

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

Publisher
Unpaywall
ISSN
1023-5809
DOI
10.5194/npg-18-7-2011
Publisher site
See Article on Publisher Site

Abstract

Nonlin. Processes Geophys., 18, 7–28, 2011 Nonlinear Processes www.nonlin-processes-geophys.net/18/7/2011/ doi:10.5194/npg-18-7-2011 in Geophysics © Author(s) 2011. CC Attribution 3.0 License. A statistical mechanical approach for the computation of the climatic response to general forcings 1,2,3 3 V. Lucarini and S. Sarno Department of Meteorology, University of Reading, Reading, UK Department of Mathematics, University of Reading, Reading, UK Walker Institute for Climate System Research, University of Reading, Reading, UK Received: 4 August 2010 – Revised: 18 November 2010 – Accepted: 4 December 2010 – Published: 12 January 2011 Abstract. The climate belongs to the class of non- equations for such parameters allow to define such properties equilibrium forced and dissipative systems, for which most as an explicit function of the unperturbed forcing parameter results of quasi-equilibrium statistical mechanics, including alone for a general class of chaotic Lorenz 96 models. We the fluctuation-dissipation theorem, do not apply. In this pa- then verify the theoretical predictions from the outputs of the per we show for the first time how the Ruelle linear response simulations up to a high degree of precision. The theory is theory, developed for studying rigorously the impact of per- used to explain differences in the response of local and global turbations on general observables of non-equilibrium statisti- observables, to define the intensive properties of the system, cal mechanical systems, can be applied with great success to which do not depend on the spatial resolution of the Lorenz analyze the climatic response to general forcings. The crucial 96 model, and to generalize the concept of climate sensitivity value of the Ruelle theory lies in the fact that it allows to com- to all time scales. We also show how to reconstruct the linear pute the response of the system in terms of expectation values Green function, which maps perturbations of general time of explicit and computable functions of the phase space aver- patterns into changes in the expectation value of the consid- aged over the invariant measure of the unperturbed state. We ered observable for finite as well as infinite time. Finally, choose as test bed a classical version of the Lorenz 96 model, we propose a simple yet general methodology to study gen- which, in spite of its simplicity, has a well-recognized pro- eral Climate Change problems on virtually any time scale totypical value as it is a spatially extended one-dimensional by resorting to only well selected simulations, and by tak- model and presents the basic ingredients, such as dissipa- ing full advantage of ensemble methods. The specific case of tion, advection and the presence of an external forcing, of globally averaged surface temperature response to a general the actual atmosphere. We recapitulate the main aspects of pattern of change of the CO concentration is discussed. We the general response theory and propose some new general believe that the proposed approach may constitute a mathe- results. We then analyze the frequency dependence of the re- matically rigorous and practically very effective way to ap- sponse of both local and global observables to perturbations proach the problem of climate sensitivity, climate prediction, having localized as well as global spatial patterns. We derive and climate change from a radically new perspective. analytically several properties of the corresponding suscep- tibilities, such as asymptotic behavior, validity of Kramers- Kronig relations, and sum rules, whose main ingredient is the 1 Introduction causality principle. We show that all the coefficients of the leading asymptotic expansions as well as the integral con- A crucial goal in the study of general dynamical and statisti- straints can be written as linear function of parameters that cal mechanical systems is to understand how their statistical describe the unperturbed properties of the system, such as properties are altered when we introduce a perturbation re- its average energy. Some newly obtained empirical closure lated to changes in the external forcing or in the value of some internal parameters. The ability to compute the re- sponse of the system is of great relevance for purely math- ematical reasons as well as in many fields of science and Correspondence to: V. Lucarini technology. (v.lucarini@reading.ac.uk) Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union. 8 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings The climate system is an outstanding example of a non- lent. Whereas natural fluctuations of the system are restricted equilibrium, forced and dissipative complex system, forced to the unstable manifold, because, by definition, asymptot- in first instance by spatial differences and temporal variabil- ically there is no dynamics along the stable manifold, ex- ity in the net energy flux at the top of the atmosphere. On a ternal forcings will induce motions – of exponentially de- macroscopic level, as a result of being far from equilibrium, caying amplitude – out of the attractor with probability one. the climate system behaves as an engine, driven by the tem- The fluctuation-dissipation relation can be recovered only if perature difference between a warm and a cold thermal pool, we consider perturbations with the somewhat artificial prop- so that the atmospheric and oceanic motions are at the same erty of being everywhere tangent to the unstable manifold or, time the result of the mechanical work (then dissipated in a in a more fundamental way, if we add a stochastic forcing, turbulent cascade) produced by the engine, and are processes which has the crucial effect of smoothing the invariant mea- which re-equilibrate the energy balance of the climate sys- sure (Lacorata and Vulpiani, 2007; Marini Bettolo Marconi tem (Lorenz, 1967; Peixoto and Oort, 1992; Johnson, 2000; et al., 2008). Potential links to these issues can be found in Lucarini, 2009a). recent papers proposing new algorithms for three (Trevisan A primary goal of climate science is to understand how the and Uboldi, 2004) and four (Trevisan et al., 2010) dimen- statistical properties – mean values, fluctuations, and higher sional variational data assimilation, where it is shown that order moments – of the climate system change as a result of the quality of the procedure improves if the increment of the modulations to some crucial external (e.g. solar irradiance) variables due to the assimilation is performed only along the or internal (e.g. atmospheric composition) parameters of the unstable manifold. system occurring on various time scales. A large class of Recently, Ruelle (1998a, 2009) introduced a rigorous problems – those involving climate sensitivity, climate vari- mathematical theory allowing for computing analytically, ability, climate change, climate tipping points – fall into this ab initio, the response of a large class of non-equilibrium category. In a system as complex and as extended as the cli- systems to general external perturbations featuring arbitrary mate, where lots of feedbacks are active on a variety of spa- time modulation. The crucial result is that the changes in tial and temporal scales, this is in general a very difficult task. the expectation value of a physical observable can be ex- The need for scientific advance in this direction is outstand- pressed as a perturbative series in increasing powers of the ing as one considers that even after several decades of intense intensity of the external perturbation, where each term of the scientific efforts, the accurate evaluation of the climate sensi- series can be written as the expectation value of some well- tivity par excellence, i.e., the change of the globally averaged defined observable over the unperturbed state. In a previ- surface temperature for doubling of CO concentration with ous paper (Lucarini, 2008b) we showed that the Ruelle the- respect to pre-industrial levels (280 ppm to 560 ppm circa), is ory is, thanks to this property, formally analogous to usual a tantalizing endeavor, and large uncertainties are still present Kubo response theory (Kubo, 1957), which applies for quasi- (IPCC, 2007). equilibrium system. The crucial difference lies on the math- Such efforts have significant relevance also in the context ematical properties of the invariant measure, which is abso- of the ever-increasing attention paid by the scientific commu- lutely continuous in the quasi-equilibrium case and singular nity to the quest for reliable metrics to be used for the valida- in the non-equilibrium case. tion of climate models of various degrees of complexity and Ruelle’s analysis applies for non-equilibrium steady state for the definition of strategies aimed at the radical improve- systems (Gallavotti, 2006) possessing a Sinai-Ruelle-Bowen ment of their performance (Held, 2005; Lucarini, 2008a). (SRB) invariant measure, often referred to as Axiom A sys- The pursuit of a quantum leap in climate modelling – which tem (Eckmann and Ruelle, 1985; Ruelle, 1989). This class definitely requires new scientific ideas rather than just faster of systems, even if mathematically non-generic, includes on supercomputers – is becoming more and more of a key issue the other hand excellent models for general physical systems, in the climate community (Shukla et al., 2009). as made clear by the chaotic hypothesis (Gallavotti and Co- A serious, fundamental difficulty in the analysis of the non hen, 1995; Gallavotti, 1996), which can be interpreted as an equilibrium systems is that the fluctuation-dissipation rela- extension of the ergodic hypothesis to non-equilibrium sys- tion (Kubo, 1966), cannot be applied (Ruelle, 1998a). This tems (Gallavotti, 2006). We also remind that any time in an greatly limits the ability of understanding the response of the numerical integration we assume that the time average of a systems to external perturbations by looking at its variabil- given observable, after discarding an initial transient, basi- ity. In the specific case of climate, this can be rephrased by cally coincides with its average over the invariant measure saying that climate change signals need not project on the giving the attractor of the system, we are actually assuming natural modes of climate variability. The non-equivalence Axiom A-like hypothesis. See Penland (2003) for an original between free and forced climate fluctuations had been sug- geophysical perspective. gested by Lorenz (1979). The basic reason for this behav- The Ruelle response theory, with the support of chaotic ior is that, since the dynamics is forced and dissipative, with hypothesis, has interesting conceptual implications for cli- the asymptotic dynamics taking place in a strange attractor, mate studies. In fact, the possibility of defining a response natural fluctuations and forced motions cannot be equiva- function basically poses the problem of climate response to Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 9 forcings and of climate change in a well-defined context, and, from the use of the fluctuation-dissipation theorem in a non- when considering the procedures aimed at improving climate equilibrium context. As discussed above, this just cannot be models, justifies rigorously the procedures of tuning and ad- rigorously the case. Since they use, albeit in a convolute way, justing the free parameters. Moreover, the response theory the forcing term to compute (or at least to parameterise) the allows to compute the climate sensitivity, in the special case response of the system, what they do is actually to compute when static perturbations to the system parameters are con- (or at least estimate) the linear response. Whereas it is in the sidered. very nature of a perturbation theory the possibility of com- Previously, a response formula was proposed by Cacuci puting the response starting only from the statistical proper- for evaluating the linearized change of the solution of a time- ties of the unperturbed state (as done also in Abramov and independent generic system of nonlinear equations as a result Majda, 2007; Majda and Wang, 2010), the specific property of a change in the system’s parameters (Cacuci, 1981a,b). of the systems where the fluctuation-dissipation theorem can This can be interpreted as a special case of Ruelle’s theory, be applied is that the linear Green function can be written in where the unperturbed attractor is constituted by a fixed point the specific form of a correlation of two observables of the and a static perturbation to the system evolution equation is system. See the discussion in Ruelle (1998a, 2009). considered. Cacuci proposed to study this problem using the It is also important to note that, following the pioneer- adjoint operator to the original system, which provided an ing study by Leith (1975), several recent studies (Langen efficient way to determine the impact of small perturbations. and Alexeev, 2005; Gritsun and Branstator, 2007; Ring and Interestingly, early prominent applications of the so-called Plumb, 2008; Gritsun et al., 2008) have attempted with a adjoint method and its extension to time-dependent prob- certain (and sometimes good) degree of success to recon- lems, which allowed for evaluating all possible linear sen- struct the climate response to external perturbations start- sitivities of an evolving model in just one simulation, were ing from the internal fluctuations of the system by using a been proposed for climate related problems. In particular, it severely simplified version of the fluctuation dissipation the- was used to evaluate the sensitivities of a simple radiative- orem, based upon the assumption of a quasi-gaussian proba- convective model possessing an attractor constituted by just bility distribution for the system (see discussion in Abramov one fixed point (Hall et al., 1982; Hall and Cacuci, 1983), and Majda, 2007). The overall good results seem to suggest and later, in an inherently heuristic way, for studying the re- that, at least for the considered models, observables (typi- sponse of a (chaotic) simplified general circulation model to cally very large scale ones) and baseline climate conditions, doubling of the CO concentration (Hall, 1986). Whereas the fractal nature of the invariant measure of the system and the adjoint method did not find much space in further cli- the role of the stable direction of the flow seem not to be ex- matic studies, mostly due to early discouragement for the ceedingly relevant. This may be related to the specific choice computational burden of constructing the suitable operators of the observable or to the fact that internally generated nu- for evaluating the sensitivities, it subsequently reached great merical noise mimics the effect of stochastic perturbations, success in data assimilation problems for geophysical fluid but further studies are surely necessary. dynamics (Ghil and Malanotte-Rizzoli, 1991; Errico, 1997), In the last decade on one side a great effort has been to the point that a tangent and adjoint model compiler able directed at extending the Ruelle response theory for more to automatically generate adjoint model code was has been general classes of dynamical systems (see, e.g., Dolgopyat, introduced (Giering and Kaminski, 1998). More recently, a 2004; Baladi, 2007), and recent studies (Lucarini, 2008b) link between advanced adjoint techniques and the Ruelle the- have shown that, thanks only to the causal nature of the re- ory has been proposed (Eyink et al., 2004). sponse, it is possible to apply all the machinery of the the- Majda and collaborators (Abramov and Majda, 2007; Ma- ory of Kramers-Kronig (KK) relations (Nussenzveig, 1972; jda and Wang, 2010) have taken a different approach for Peiponen et al., 2005; Lucarini et al., 2003, 2005) for linear analysing the response of a system to external perturbation. and nonlinear processes to study accurately and rigorously Using the formalism of the Fokker-Planck equation, they the susceptibilities describing in the frequency domain the have developed a response theory which is analogous to Ru- response of a general observable to a general perturbation. elle’s. These authors have especially emphasized the pos- Moreover, the actual applicability of the theory has been sibility of splitting the response into the components rela- successfully tested in a number of simple dynamical systems tive to the projection of the perturbation vector flow onto the case for the linear (Reick, 2002; Cessac and Sepulchre, 2007) stable, unstable and neutral components of the unperturbed and nonlinear (Lucarini, 2009b) response. Such numerical flow, as already discussed in Ruelle (1998a), and have in- investigations have clarified that even in systems which are troduced algorithms to compute efficiently all of these com- not Axiom A, like the Lorenz 63 system (Lorenz, 1963), it is ponents. Their approach has very recently been extended possible to successfully use the response theory to construct to take into account stochastic and time-periodic perturba- linear (Reick, 2002) and nonlinear susceptibilities (Lucarini, tions (Majda and Wang, 2010). Somewhat confusingly, they 2009b) which obey all of the constraints imposed by the KK mostly refer to the rather interesting applications of their the- theory up to a high degree of precision. ory and algorithms to specific dynamical systems as resulting www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 10 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings These investigations definitely motivate further studies Whereas the paper aims at proposing new methods for aimed at understanding to what extent the response theory tackling classical problems of climate science, most of the is an efficient tool for analyzing complex and complicated results and of the methodologies proposed are of more gen- systems. In this paper, we take up such a challenge and con- eral interest. In this paper we limit our attention to the linear sider the Lorenz 96 (L96) system (Lorenz, 1996; Lorenz and response. We refer to Lucarini (2008b, 2009b) for a theoreti- Emanuel, 1998; Lorenz, 2004), which provides an excellent cal and numerical studies of higher-order effects of perturba- and celebrated prototypical model of a one dimensional at- tions. mosphere. The variables of the L96 model can be thought The paper is organized as follows. In Sect. 2 we briefly as generic meteorological quantities extending around a lat- analyze the general theoretical background of the linear re- itudinal circle and sampled at a regular interval. In spite of sponse theory and of the properties of the frequency depen- not being realistic in the usual sense, the L96 model presents dent susceptibility and present some new useful results. In the basic ingredients, such as dissipation, advection and the Sect. 3 we present the main features of the L96 system, in- presence of an external forcing, of the actual atmosphere. For troduce the considered perturbations to the forcing, derive this reason, L96 has quickly become the standard model to some basic properties of the response of various observables, be used for predictability studies (Orrell, 2003; Haven at al., and present the theoretical predictions. In Sect. 4 we present 2006; Hallerberg et al., 2010), when testing data assimila- the results of our numerical investigations and describe how tion techniques (Trevisan and Uboldi, 2004; Trevisan et al., they can be generalized to the entire family of L96 models. 2010; Fertig et al., 2007), and new parameterizations (Wilks, In Sect. 5 we provide a relevant example to illustrate how the 2006). The L96 model had already been taken as test-bed for results presented in this paper can be used to devise simple studying the linear response (the applicability of the fluctua- yet rigorous methods to study the climate response at all time tion dissipation theorem, in their language) in Abramov and scales on models of any degree of complexity. In Sect. 6 we Majda (2007). discuss the conclusions and present perspectives for future Although we are unable to prove that the unperturbed L96 work. is an Axiom A system, in general and for the specific choice of parameters used in our numerical simulations in particular, 2 Theoretical background: Ruelle theory and we adopt the chaotic hypothesis and present the first thor- dispersion relations ough investigation of a spatially extended system by using the rigorous statistical mechanical methodologies presented 2.1 Definition of the linear susceptibility in Ruelle (1998a, 2009) and Lucarini (2008b, 2009b). More- over, since L96 is a spatially extended system, we also ex- We consider an Axiom A dynamical system described by the plore the applicability of the response theory in all possi- evolution equation x˙ = F (x), so that the invariant probabil- ble combinations of global/local observables and global/local ity measure ρ of the associated flow is of the SRB type perturbations. We compute rigorously the corresponding lin- (Ruelle, 1998a). Let h8i be the expectation value of the ear susceptibilities, verify the KK relations and the related general observable 8 defined as ρ (dx),8(x). We perturb sum rules, and find an empirical power law. This, as in the the flow of the system by adding a on the right hand side of case of discussed in Lucarini et al. (2007), supports the valid- the evolution equation a vector field X(x)f(t), where X(x) ity of the chaotic hypothesis, allowing to extend the results defines the pattern of the perturbation, and f(t) is its time obtained for our specific choice of model’s parameters to a modulation. The resulting evolution equation results to be rather general class of L96 systems. We also show how to x˙ = F(x) + X(x)f(t). Following Ruelle (1998a), we ex- go from the frequency back to the time domain, thus deriv- press the expectation value of 8(x) in the perturbed system ing from the susceptibility the Green function, which acts as using a perturbative expansions as: time propagator of the considered perturbation for the con- sidered observable. The Green function allows to predict, X (n) h i h i h i 8 (t) = 8 + 8 (t). (1) in an ensemble mean sense, the change in the observable at 0 n=1 any time horizon as a result of a perturbation with the same spatial patter as that considered in the calculation of the sus- Each term of the perturbative series can be expressed as an ceptibility but featuring a general time modulation. th n-convolution integral of the n order causal Green function Finally, building upon the results presented here, we pro- with n delayed time perturbation functions (Ruelle, 1998b; pose a simple yet general methodology to study general Cli- Lucarini, 2008b). Limiting our attention at the linear case mate Change problems on virtually any time scale by resort- we have: ing to only well selected simulations, and by taking full ad- +∞ vantage of ensemble methods. The specific case of globally (1) (1) h8i (t) = dσ G (σ )f(t −σ ) (2) 1 1 1 averaged surface temperature response to a general pattern of −∞ change of the CO concentration is discussed. Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 11 The first order Green function can be expressed as follows: Green function is causal. Assuming, on heuristic physical ba- (1) sis, that G (t) ∈ L , we can apply the Titchmarsh theorem (1) 8 G (σ ) = ρ (dx)2(σ )35 (σ )8(x), (3) 1 0 1 0 1 8 (Nussenzveig, 1972; Peiponen et al., 2005; Lucarini et al., (1) 2003, 2005) and deduce that the linear susceptibility χ (ω) where 3(•) = X(x)·∇(•) takes into account the effects of is a holomorphic function in the upper complex ω-plane and the perturbative vector field, 2 is the usual Heaviside dis- the real and the imaginary part of χ(ω) are connected to each tribution and 5 the unperturbed time evolution operator so other by Hilbert transform. that 5 K(x) = K(x(t)) for any function K , with x(t) fol- According to a general property of Fourier transform we lowing the unperturbed flow. Note that it is possible to ex- (1) know that the short term behavior of G (t) determines the press the Green function as the expectation value of a non- (1) asymptotic properties of χ (ω). We shall obtain a more trivial but computable observable over the unperturbed SRB quantitative result by exploiting that: measure ρ . Therefore the knowledge of the unperturbed fea- tures of the flow is sufficient to define the effects of any ex- Z +∞ k d i k k ternal perturbation over any observable of our system. In the dt2(t)t exp[iωt] =(−i) P +πδ(ω) dω ω −∞ frequency domain we find that the first term of the perturba- (k+1) tive series can be written as: i Z ≈k! (6) +∞ (k+1) (1) (1) h8i (ω) = dω χ (ω )f(ω )×δ(ω−ω ) 1 1 1 1 −∞ where in the second equality we have neglected the fact that (1) the solution is a distribution and considered ω 6= 0. There- =χ (ω)f(ω), (4) fore, if the Taylor expansion of the Green function in the where the Dirac delta implies that we are analyzing the im- limit t → 0 is of the form: pact of perturbations in the frequency-domain at the fre- (1) quency ω. The linear susceptibility is defined as: β β G (t) ≈α¯ 2(t)t +o(t ) (7) +∞ (1) (1) χ (ω) = dtG (t)exp[iωt]. (5) the high frequency behavior of the linear susceptibility, i.e. 8 8 −∞ the limit ω → ∞, is: It is important to underline with a thought experiment the (1) −β−1 −β−1 computational relevance of the last equations and the impor- χ (ω) ≈αω +o(ω ) (8) tance of the susceptibility function. (β+1) Let’s suppose we introduce a time dependent perturbation where α =α¯i β!. The parameters β (which is an integer f (t) to a given pattern of forcing X(x), simulate the sys- α number) and α¯ depend on the observable 8, on the specific tem and observe the time response of an arbitrary observ- features of the unperturbed system, and on the forcing un- (1) able h8 i (t). We now compute the Fourier transform of der consideration. Taking into account Eq. (5) and assuming (1) (1) the observed signal and of the forcing modulation. Invert- that ω is real, we obtain that χ (ω) = [χ (−ω)] , so that 8 8 (1) ing Eq. (4), we can find the linear susceptibility χ (ω). Re[χ] is an even function while Im[χ] is odd function of ω. Let’s now consider a different time-modulating function of Thus α =α is real if β is odd, whereas α = iα is imaginary R I the forcing f (t) and its corresponding Fourier transform β if β is even. f (ω). Taking into account Eq. (4), if we multiply f (ω) β β Taking into account the Titchmarsch theorem, using that (1) (1) (1) times the previously computed function χ (ω) we directly χ (ω) = [χ (−ω)] , and considering the asymptotic be- ¯ 8 8 (1) havior of the susceptibility, it is possible to show that the real obtain h8 i (ω), the frequency-dependent response of the and imaginary part of the linear susceptibility obey the fol- observable 8 to the forcing X(x) modulated by the new lowing set of general KK dispersion relations (Lucarini et al., function. By applying the inverse Fourier transform we ob- (1) 2005): tain the time-dependent response h8 i (t) without needing any additional simulation. (1) 2p π ν Re[χ (ν)] (1) Moreover, the knowledge of the susceptibility function al- 2p−1 8 − ω Im[χ (ω)] =P dν (9) (1) 2 2 2 ν −ω lows us to reconstruct the G (t) by inverting Eq. (5). Oth- 0 erwise, the Green function can be obtained directly from observing the response signal by performing a simulation (1) ∞ 2p+1 π ν Im[χ (ν)] where f(t) =δ(t): in this case (see Eq. (2)) we simply have (1) 2p 8 ω Re[χ (ω)] =P dν (10) (1) 8 (1) 2 2 2 ν −ω h8i (t) =G (t). with P indicating integration in principal part and p = 2.2 Kramers-Kronig relations and sum rules 0,...,(β − 1)/2 if β is odd and p = 0,...,β/2 if β is even. As we see from Eq. (3), for an arbitrary choice of the ob- Note that the faster the asymptotic decrease of the suscepti- servable and of the perturbation the corresponding linear bility, the higher the number of independent constraints due www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 12 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings to KK relations it has to unavoidably obeys. As thoroughly discussed in Kubo (1966), in the case of quasi-equilibrium (1) system, the fluctuation-dissipation theorem ensures that the ν Re[χ (ν)]dν = α , (15) imaginary part of the susceptibility describing the response of a given observable to a perturbation is proportional to where the α constants are defined in Eq. (8). These sum rules a suitably defined power spectrum in the unperturbed sys- provide additional general constraints that must be obeyed by tem. Therefore, observing the unperturbed system and using any system and can be used to test the quality of the output of Eq. (10) it is possible to reconstruct the entire linear suscepti- any model wishing to describe it. If we are not in the perfect bility, and so know everything about the response properties model scenario (e.g., we use a simplified representation of of the system. In the case of a non-equilibrium system, as some degrees of freedom) the sum rules can in principle be discussed in the Introduction, this procedure is not possible. used to provide a fit for the parametrization. It is possible to use the KK relations to define specific self- We underline that it is possible to generalize the KK the- consistency properties of the real and imaginary part of the ory for specific classes of nonlinear susceptibilities for both susceptibility. We first consider the following application: quasi-equilibrium and non-equilibrium systems. Such re- we set p = 0 in Eq. (10) and take the limit ω → 0. We obtain sults, which are particularly suited for studying the funda- that for any observable: mental properties of harmonic generation processes, are thor- oughly discussed in Lucarini (2008b) and will not be re- (1) 2 Im[χ (ν)] (1) Re[χ (0)] = P dν , (11) ported here. π ν 2.3 A practical formula for the linear susceptibility and which says that the static susceptibility (i.e., in a more com- consistency relations between susceptibilities of mon language, the linear sensitivity of the system) is related different observables to the out-of-phase response of the system at all frequencies. In other terms, Eq. (11) is an exact formula for the linear sus- As discussed above, the definition of the linear susceptibility ceptibility of the system. Note that the static susceptibility does not depend on the function f(t) modulating the addi- is a real number because, thanks to the symmetry proper- (1) tional forcing, so that it is possible to draw general conclu- ties discussed above, Im[χ (0)] = 0. Moreover, we know (1) sions on its properties even by choosing a specific function that Re[χ (0)] is finite because the susceptibility function f(t). is analytic (and so in particular non singular). This is con- Let’s consider f(t) = 2 cos(ωt). The impact of the per- sistent with the fact that as ω → 0, the imaginary part of the turbation on the evolution of a general observable 8(x) is susceptibility goes to zero at least as fast as a linear func- defined as: tion, as only odd positive integer exponents can appear in its Taylor expansion around ω = 0), so that the integrand in δ8 (t,t ,x ) =8 (t,t ,x )−8 (t,t ,x ) (16) 0 0  0 0 0 0 0 Eq. (11) is not singular. Similarly, we obtain that for ω ∼ 0 the real part of the susceptibility is in general of the form where x and t are the initial condition and the initial 0 0 2 3 c +c ω. +o(ω ), where c and c are two constants and c 1 2 1 2 1 time, and we associate the lower index  to the strength is exactly given by Eq. (11). of the forcing. The Ruelle’s response theory ensures that By exploring the ω → ∞ limit in Eqs. (9–10) we obtain δ8 (t,t ,x ) = O(). Following Reick (2002) and Lucarini 0 0 further integral constraints. By applying the superconver- (2009b) the linear susceptibility results to be: gence theorem (Frye and Warnock, 1963), we obtain the fol- (1) (1) lowing set of vanishing sum rules (see Lucarini et al., 2003): χ (ω) ≡ lim lim χ (ω,x ,,T ) (17) 8 φ →0T →∞ 0 ≤p ≤β/2− 1, β even (1) 2p+1 dνν Im[χ (ν)] = 0 . (12) 8 where: 0 ≤p ≤(β − 3)/2, β odd 1 1 (1) ∞ χ (ω,x ,,T ) = dtδ8 (t,x )exp(iωt) (18) 0  0 0 ≤p ≤β/2− 1, β even (1) 2p ν Re[χ (ν)]dν = 0 . (13) 0 ≤p ≤(β − 1)/2, β odd is the total susceptibility, affected by the finite time and finite Note that if β = 0 no vanishing sum rules can be written for size response of the system. This quantity depends on the the susceptibility, whereas if β = 1 only Eq. (13) provides a initial condition and in principle contains information about zero-sum constraint. For each set of KK relations, an addi- the response of the system at all order of nonlinearity. tional, non-vanishing sum rule can be obtained. If β is odd, Since d(δ8 (t,x ))/dt = δ(d8 (t,x ))/dt , thanks to the 0  0 the non-vanishing sum rule is: linearity of the time derivative, by considering Eqs. (17–18) and performing an integration by parts, we obtain that (Re- (1) ν Im[χ (ν)]dν = − α , (14) ick, 2002): (1) (1) while if β is even, we have: χ (ω) = −iωχ (ω). (19) ˙ 8 Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 13 (1) Let’s now find a different expression for χ (ω). The time Moreover, since, as shown in Eq. (7–8), the susceptibility derivative of 8 in the unperturbed system is of a generic observable decreases to zero at least as fast as (1) −1 ω , for large values of ω we have that χ (ω) ≈ i/ωh4i 8(x) =0(x), (20) unless ωh4i = 0. If ωh4i 6= 0, we also have that the lead- 0 0 ing order of the short-time expansion of the Green function where 0 = F ·∇8. Similarly, the time derivative for 8 in is of the form: the case of the perturbed motion described by x˙ = F(x)+ X(x)f(t) = F(x)+ 2 cos(ωt)X(x) is: (1) G (t) =2(t)h4i +o(t ), (26) 8(x) =0(x)+ 2 cos(ωt)4(x), (21) in agreement with what can be found by direct inspection of Eq. (3). where 4 = X·∇8. From Eqs. (20–21) we obtain that δ8 (t,x ) =δ0 (t,x )+ 2 cos(ωt)4(t,x ), (22) 0  0 0 3 Application of the response theory to the Lorenz where all terms are of O(). Furthermore, we integrate each 96 model term in Eq. (22) as in Eq. (18), take the limits  → 0 and T → ∞, and obtain: 3.1 Statistical properties of the unperturbed Lorenz 96 Model 1 1 lim lim dtδ8 (t,x )exp(iωt) T →∞ →0 T 0 The Lorenz 96 model (Lorenz, 1996; Lorenz and Emanuel, 1 1 1998; Lorenz, 2004) describes the evolution of a generic at- = lim lim dtδ0 (t,x )exp(iωt) mospheric variable defined in N equally spaced grid points →0T →∞T along a latitudinal circle and provides a simple, unrealistic 1 1 + lim lim dt4 (t,x )[ exp(iωt) but conceptually satisfying representation of some basic at- →0T →∞T mospheric processes, even if such one-dimensional model it +exp(−iωt) ]exp(iωt). (23) cannot be derived ab-initio from any dynamic equation via subsequent approximations. The evolution equations can be Using the definition in Eq. (17) and the identity given in written in a scaled form as follows: Eq. (19), we derive: dx =x (x −x )−x +F (27) i−1 i+1 i−2 i (1) (1) dt −iωχ (ω) =χ (ω))+ lim lim dt4 (t,x ) 8 0 →0T →∞T where i = 1,2,.....,N , and the index i is cyclic so that x = Z i+N 1 1 x = x . The quadratic term in the equations simulates i−N i + lim lim dt4 (t,x )exp(2iωt). →0T →∞ advection, the linear one represents thermal or mechanical (24) damping and the constant one is an external forcing. The evolution equations are invariant under i →i + 1, so that the The first limit in Eq. (24) gives, by definition, h4i , whereas dynamics is the same for all variable. The time scale of the the second limit vanishes as the expression under integral is 2 system is given by the damping time, which corresponds to O( ), since it is related to second order harmonic gener- five days. The L96 system shows different features, as dif- ation nonlinear process (Lucarini, 2009b). Concluding, we ferent choices of F and N may strongly alter the topology of obtain the following general consistency relation for the lin- the attractor, alternating periodic, quasi-periodic and chaotic ear susceptibility: behavior in a non trivial way. However with a suitable choice i i (1) (1) of the parameters N and F , the system is markedly chaotic. χ (ω) = χ (ω)+ h4i . (25) 8 0 ω ω In particular, as F controls the energy input into the system, we expect that for relatively high values of this parameter Such an identity related the susceptibility of an observable the system should simulate a turbulent behavior and live on 8 to the susceptibility of the projection of its gradient along a strange attractor. As an example, setting N = 40 and F = 8 the unperturbed flow 0 and to the average value in the unper- the system possesses 13 positive Lyapunov exponents, the turbed state of the projection of its gradient along the pertur- largest corresponding to a doubling time of 2.1 days, while bation flow. Note that the two terms on the right hand side the fractal dimension of the attractor (Ruelle, 1989) is about are radically different. Whereas the first term is related to the 27.1 (Lorenz and Emanuel, 1998). projection of the dynamics along the unstable manifold, the When computing the time derivative of the total energy of second term depends on the structure of the forcing X(x), P the system, defined as E = 1/2 x , the advection terms i i which may be entirely unrelated to that of the unstable man- cancel. The evolution equation for E results to be: ifold. This is the fundamental reason why the fluctuation- dissipation theorem does not apply in the non-equilibrium E = −2E +F x . (28) case. www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 14 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings As the dynamics takes place inside a compact set, 9(x ,...,x ) is bounded for any choice of the function 9 . i N 1+γ Therefore the ensemble mean with respect the ρ (or time e(F) = F , (36) average) of the temporal derivative 9 vanishes. Therefore, defining M = x as the total momentum of the system, with λ ≈ 1.15 and γ ≈ 0.35. Such a smooth dependence of we obtain the following identity: the intensive energy and momentum with respect to the forc- D E X X ing parameter F is indeed in agreement with the hypothesis 2hEi = x =F hx i =F hMi . (29) 0 i 0 0 that the invariant measure is deformed in a very regular fash- i i ion not only locally, but over a large range of the parameter’s Similarly, we can deduce an additional consistency relation space. by investigating the expression of the time derivative of M : Note that, at the fixed point of the system corresponding to a purely zonally symmetric dynamics (x = F , ∀i) we have hMi =NF +hC i −hC i (30) 2 3 0 0 0 m =F and e =F /2. These formulas give much higher val- P P where hC i = hx x i and hC i = hx x i . 2 i i−2 3 i i−3 0 0 0 ues for both m and e than what found with our empirical i i Higher order consistence relations can be obtain in a simi- power laws for the attractor in the chaotic regime. In fact, lar fashion. at such an equilibrium, which is unstable in the parametric The equivalence of all the variables implies that over the range explored here, the energy dissipation is much weaker unperturbed flow each observable O of the whole system sat- than in the co-existing chaotic attractor, which corresponds to isfies hO(x )i = N O(x ) ∀j . Therefore, we define i j the case where breaking nonlinear waves and turbulent mo- i 0 the average energy per grid point e(N,F) and the average tions are present. Interestingly, the presence of well-defined momentum per grid point m(N,F) as: scaling laws with respect to the forcing parameters for the 2 energy and momentum of the system with different charac- hEi i 0 e = = , (31) teristic exponents in the chaotic regime and in the co-existing 2 N unstable equilibrium is in agreement with previous finding recently obtained in a simple baroclinic quasi-geostrophic hMi model (Lucarini et al., 2007). m = hx i = , (32) i 0 3.2 Asymptotic properties of the linear susceptibility where the choice of i is arbitrary and the N and F - dependence is dropped for shortness. Defining c = We perturb the L96 model by adding a small perturbation hC i /N , c = 2e, we can rewrite Eqs. (29–30) as follows: i 0 0 modulated by f(t) = 2 cos(ωt). The resulting evolution equation is: 2e =F m, (33) dx =x (x −x )−x +F + 2 cos(ωt)X (37) i−1 i+1 i+2 i i dt m =F +c −c . (34) 2 3 where X = X (x ,...,x ) is a generic function of vari- i i 1 N Expressing either the average energy e or the average mo- ables x . We adopt the chaotic hypothesis (Gallavotti, 1996; mentum m per grid point as a function of the two free pa- Gallavotti and Cohen, 1995) and we follow the theory pro- rameters N and F would allow to get a closure equation for posed by Ruelle (1998a) and discussed in Sect. 2 in order the statistical properties of the unperturbed Lorenz 96 sys- to study the linear response of suitably defined observables tem. We have computed e(N,F) and m(N,F) by performing to the perturbation. We first propose to study the high- long integrations for values of F ranging from 6 to 50 with frequency, response by analyzing in detail the asymptotic step 1 and for values of N ranging from 10 to 200 with step properties of the resulting susceptibilities. As discussed in 10. In all of these cases, chaotic motions are observed. We Sect. 2, this constitutes a crucial step for constructing the set have consistently found that, within 0.5%, e(N,F) = e(F) of applicable KK relations and for computing the value of the and m(N,F) = m(F), so that they can be considered inten- sum rules. sive quantities. Therefore, we can interpret Eqs. (33–34) as We consider two different forcing patters X . In the first equations providing a definition of the thermodynamics of case, we apply the perturbation over all the grid points and this simple one-dimensional model of atmosphere. we choose X = 1∀i. Given an observable 8, we refer to the The F -dependence of e and m can be closely approxi- linear Green function and linear susceptibility resulting from (1) (1) mated in terms of simple power laws. We obtain, within a this choice of X as G and χ , respectively, where the 8,N 8,N precision of about 1% in the considered domain, that lower index N indicates that the perturbation acts over all γ the variables x . We refer to this pattern of forcing as global m(F) =λF , (35) perturbation. In the second case, we apply the perturbation and, consistently with Eq. (33), only on the variable x of the L96 model, and we choose Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 15 X = 0∀i 6= j , X = 1. Since all the points are equivalent in Equations (40) and (41) imply that the imaginary part 1 j the unperturbed case, the choice of j is arbitrary. In this case, dominates the asymptotic behavior of the susceptibility, so when referring to the linear Green function and the linear that at high frequency the response is shifted by about π/2 susceptibility, the lower index 1 substitutes the N , indicating with respect the forcing. Observing that the leading term of −1 that the perturbation is localized to one point. We refer to asymptotic χ is of order ω just one sum rules apply for this pattern of forcing as local perturbation. either susceptibilities. Limiting our attention to the intensive quantity e, by applying Eq. (15) we obtain: 3.2.1 Global perturbation (1) Re[χ (ω)]dω = m. (42) ε,N We consider perturbations with spatial pattern given by X = i 2 1∀i and analyze the response of the observable E. Following Along similar lines, if we select as observable the total mo- (1) Eq. (3), the linear Green function G (t) can be explicitly E,a mentum M , we derive that the asymptotic behavior of its lin- written as: ear susceptibility is: (1) G (t) = ρ (dx)2(t)35 (t)E(x) 0 0 N N E,N (1) (1) −2 χ (ω) =Nχ (ω) = i + +o(ω ), (43) M,N μ,N Z 2 ω ω = ρ (dx)2(t)1·∇E(x(t)) (1) where we have defined the intensive susceptibility χ (ω), Z μ,N where μ is the intensive momentum of the system. As in = ρ (dx)2(t) ∂ (E(x(t))) (38) 0 i Eq. (41), the asymptotic behavior is determined by the imag- inary part of χ , and the real part of the susceptibility provides where x(t) satisfies the unperturbed evolution Eq. (27). Tak- the following sum rule: ing into account Eq. (7) and Eq. (8), in order to obtain the (1) asymptotic behavior of the susceptibility, we need to study Re[χ (ω)]dω = . (44) μ,N the short time behavior of the Green function. Therefore, we express E(x(t)) as a Taylor series about t = 0 considering We now wish to go back to the general consistency equation the unperturbed flow, compute the integral of each coeffi- for linear susceptibilities given in Eq. (46). Considering that cient of the t -expansion over ρ , and seek the lowest order in the perturbed system the time derivative of the total energy non-vanishing term (Lucarini, 2009b). The first two terms of of the system can be written as: the Taylor expansion of E in Eq. (38) give: E = −2E +FM + 2 cos(ωt)M (45) (1) G (t) = ρ (dx)2(t) ∂ E| +tE| +o(t) 0 i t=0 t=0 E,N i the general result given in Eq. (25) can be written as follows: h   i X X = ρ (dx)2(t) x − 2x −NF t +o(t) . (39) 0 i i F 1 i i (1) (1) χ (ω) = χ (ω)+ hMi , (46) E,N M,N 2− iω 2− iω Using Eqs. (7–8), the leading terms of the asymptotic behav- ior of linear susceptibility can be written as: since in this case 0 = −2E +FM and 4 = M . It is easy X X to check that the asymptotic behavior for the susceptibilities (1) χ (ω) =i hx i /ω+ h2x i −NF /ω i i 0 0 given in Eqs. (40–43) is in agreement with Eq. (46), which is E,N i i valid at all frequencies. m F − 2m −2 −2 +o(ω ) =iN −N +o(ω ) (40) 3.2.2 Local perturbation ω ω Since the symmetry with respect the index i is valid also in We now perturb the system in a single grid point. The sym- the perturbed case, given our choice of the forcing pattern, metry of the unperturbed system implies the equivalence of the linear susceptibility of the total energy is given the sum every point of the latitude circle. Indicating with x the grid of N identical contributions, each corresponding to the sus- point where forcing is exerted, the pattern of the perturbation ceptibility of the observable ε = 1/2x for each of the N vari- vector field is X = 0∀i 6= j , X = 1. We consider the same 1 j ables x of the system. Therefore, it is possible to define an modulating monochromatic function f(t) = 2 cos(ωt) as in (1) (1) (1) intensive linear susceptibility χ = 1/Nχ , where χ the previous case. We analyze the asymptotic behavior of ε,N E,N ε,N describes the response of the local energy to the external per- the linear susceptibility of the total energy of the system E, turbation. In particular in the limit ω → ∞ we have: of the total momentum of the system M , which are global variables, and of the local variables constituted by the energy m F − 2m (1) −2 χ (ω) =i − +o(ω ). (41) E = 1/2x and momentum M = x of the perturbed grid j j j ε,N 2 j ω ω point and of its immediate neighbors. www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 16 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings Since we are looking at the linear response and the we obtain the following asymptotic behavior for the linear global perturbation is given by N spatially shifted copies susceptibility of the local perturbation, for any observable 8 of the form P m m (1) (1) (1) N −2 8(x ,...,x ) = φ(x ) we must have χ = χ = χ =i + +o(ω ). (54) 1 n i i=1 8,1 8/N,N E ,1 j 2 ω ω (1) χ . φ,N (1) (1) In the case of the observable E, it is straightforward to ver- Since, by linearity, χ = χ , comparing this result E,1 E ,1 ify the previous identity at least in the asymptotic regimes. In with what obtained in Eq. (48), we note that the suscepti- fact, the short time behavior of the Green function describing bility of the energy at the position of the forcing E pro- the response of the E to the local perturbation results to be: vides the leading asymptotic term to the susceptibility of the total energy E. Consequently, in the high-frequency range (1) (1) (1) G (t) = ρ (dx)2(t)∂ E| +tE| +o(t) 0 j t=0 t=0 χ ≈χ , and the two susceptibilities obey the same non- E,1 E ,1 E,1 vanishing sum rule, so that: h   i X X 2 2 = ρ (dx)2(t)∂ x − 2x −x F t 0 j i i i (1) Re[χ (ω)]dω = m (55) E ,1 =2(t) 2x − 2x −F t +o(t) , (47) j j 0 0 0 Nevertheless, by comparing Eqs. (48–54), we discover that so that the asymptotic behavior of the corresponding linear −2 contributions to the second leading order (∝ω ) in the high susceptibility susceptibility is: frequency range of the susceptibility of the total energy do m F − 2m not come just from the response of the energy at the per- (1) −2 χ =i − +o(ω ), (48) E,1 turbed grid point. perturbed grid point but some other point ω ω −2 give a contribution of order ω . Therefore, the asymptotic which agrees with what found for the intensive energy re- (1) (1) behavior of the real part of χ is not captured by χ . sponse when the global perturbation is applied (see Eq. (41)). E,1 E ,1 The locality of the interaction suggests to look at the energy The sum rule for the real part of the susceptibility is exactly of the closest neighbors of x . Because of the asymmetry of the same as in what given in Eq. (42): j the nonlinear terms in the L96 evolution equations, we con- (1) sider the observable ψ = 1/2(E +E +E ). It is j+1 j+2 j−1 Re[χ (ω)]dω = . (49) E,1 2 possible to prove that: Analogously, we obtain that the asymptotic behavior of x (x −x ) (F −m) j−1 j+1 j−2 (1) 0 −2 (1) χ (ω) = − = − +o(ω ). (56) χ can be written as: ψ,1 2 2 M,,1 ω ω (1) (1) 1 1 (1) It is easy to observe that the sum of χ and χ provides ψ,1 E ,1 χ (ω) = i + , (50) j M,1 ω ω the correct leading order to the asymptotic behavior of both (1) the real and imaginary parts of χ . We shall provide an with the corresponding sum rule: E,1 argument why this strongly supports the close resemblance π (1) (1) (1) of the two functions χ and χ , where E =E +ψ = Re[χ (ω)]dω = , (51) loc j E,1 E ,1 loc M,1 0 2 2 2 2 1/2x + 1/2x + 1/2x + 1/2x is the energy of the j j+1 j+2 j−1 cluster of points centered in x . in perfect agreement with Eqs. (43) and (44), respectively. j The analysis of the asymptotic behavior of the susceptibil- It is rather interesting to look into local energy observ- ities related to the local momentum of the system provides ables. Considering the energy E of the perturbed grid point additional insights. It is possible to prove that for large fre- x we have that its short term Green function can be written quency the linear susceptibility of the momentum of the per- as: turbed grid point is: (1) G (t) = ρ (dx)2(t)∂ [ x +(x (x x 0 j j j−1 j+1 E ,1 1 1 1 (1) −2 χ (ω) = i + +o(ω ), (57) x ,1 −x x −x +F ) )t +o(t) ] =2(t)[ x j j−1 j−2 j j ω ω − x x −x x −x +F t +o(t) ]. (52) j−1 j+1 j−1 j−2 j which suggests that the response of the local momentum cap- tures the correct asymptotic behavior of both the real and the Since imaginary part of the total momentum M . Concluding, we 0 = x˙ = x x −x x −x +F (53) j j−1 j+1 j−1 j−2 j obtain that the following sum rule can be stablished: 0 0 because the ρ -average of the temporal derivative of any ob- 0 π (1) Re[χ (ω)]dω = . (58) x ,1 servable vanishes, thanks to the compactness of the attractor, j Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 17 Therefore, such a constraint is exactly the same whether we t = T = 400π/ω, which corresponds to 200 full periods of analyze the the response of the momentum of a single vari- the forcing. The length of the simulations depends on the able when the perturbation acts over all the grid points, or of corresponding period of the forcing because we are inter- the total momentum in the case of a local perturbation, or, in ested in obtaining a frequency-independent quality for the this latter case, of the momentum of the grid point where the signal. We have observed that the linear response approxi- local perturbation is applied. mation is obeyed to a good degree of approximation for up As we have seen in this section, the coefficients of the lead- to  ≈ 1, which implies that the third order nonlinear effects ing asymptotic terms and the sum rules are given by simple are relatively small. See Lucarini (2008b, 2009b) for further linear functions of m (or equivalently, thanks to Eq. (33), by clarifications on this latter point. e) and by F . As we have proposed an efficient parameteri- When considering the susceptibilities describing the re- zation of m and e as functions of F alone in Sect. 3.1, our sponse to the global perturbation, we present results obtained results can be easily applied and numerically verified for a using  = 0.25 and averaging over K = 100 random initial very large class of L96 models. conditions. When assessing the linear response to the local perturbation, a reasonably clear signal is obtained using  = 1 and averaging over K = 300 initial conditions. 4 Results Note that, since we are interested in the linear response, it is could have been possible to compute the susceptibility 4.1 Simulations and data processing using a generic modulating function f(t) (see Eq. 4) rather than having to resort to multiple monochromatic perturba- The accurate calculation of the linear susceptibility of the tions. Nevertheless, for reasons of clarity, and for empha- general observable 8 is not as easy task, since the defini- sizing that chaotic dynamical systems can be analyzed using tion given in Eq. (17) requires the evaluation of two limits, tools typical of spectroscopy, we have used a more cumber- whereas we can actually compute only the quantity given in some but probably more convincing approach. Eq. (18). Averaging the response over a long time T allows We underline that the numerical results have been obtained for improving the signal-to-noise ratio. Noise is present be- using a commercial laptop rather than resorting to HPC. This cause, due to the chaotic nature of the flow, we have a contin- comes from the motivation of showing that the methodology uous spectral background. Instead, considering small values presented is robust enough that relatively low-end means al- for the perturbation strength  degrades the signal-to-noise low us to see the physical and mathematical properties of ration, but, on the other hand, it is crucial to select a small our interest. We emphasize that, using HPC, it is rather easy in order to keep the perturbations as close as possible to to greatly increase the quality of the signal by increasing K the linear regime. As discussed in Lucarini (2009b), we can and/or T by a one or two orders of magnitude. improve the signal-to-noise ratio without needing to perform very long integrations and to consider large values for  by 4.2 Global perturbation performing an ergodic averaging of the quantity averaging (1) the quantity χ (ω,x ,,T ). Therefore, we choose the best i (1) (1) We first consider χ = 1/Nχ , where ε =E/N , and fol- (1) ε,N E,N estimator of the true susceptibility χ as: low up from the discussion in Sect. 3.2.1. The measured real and imaginary parts of the susceptibility are depicted with the black lines in Fig. 1a, b. The imaginary part has a broad spec- (1) (1) χ (ω) = lim χ (ω,x ,,T ), (59) 8 8 tral feature (with two distinct peaks) spanning from ω ≈ 2 to K→∞K i=1 ω ≈ 4, which corresponds to about twice the time scale (= 1) where the x are randomly selected initial conditions chosen of the system and to four times (see Eq. 28) the relaxation on the attractor of the unperturbed system. time of the energy. This hints at the fact that it is not ob- The numerical integrations of the Lorenz 96 system have vious to constrain the spectral features of the response an been performed using the standard configuration where N , observable just by performing a scale analysis of its evolu- the number of degrees of freedom, is set to 40, and F , the tion equation. For higher values of ω, the imaginary part intensity of the unperturbed forcing, is set to 8 (Lorenz, 2004, decreases in a very regular way, so that in the upper range 1996). Equations (27) and (37) are solved using the standard a very good agreement with the asymptotic behavior ∼m/ω fourth order Runge-Kutta numerical scheme. presented in Eq. (8) is obtained. For low frequencies, the For a given observable 8, the susceptibility at angular fre- imaginary part appears to decrease towards zero, as expected quency ω is computed by applying Eq. (59) to K outputs from symmetry reasons. Instead, the real part presents a dis- of Eq. (37), each starting with a different initial condition, persive structure in correspondence with the broad maximum where the perturbation has the same angular frequency ω. of the imaginary part, and changes sign for ω ≈ 6, so that it is The angular frequency ω ranges from ω = 0.2π to ω = 20π negative for high values of the frequency range. The asymp- l h with steps of 0.01π . Each simulation performed with a totic decrease to zero in this range is also in excellent agree- perturbation of angular frequency ω runs from t = 0 up to ment with the estimate ∼ −(F − 2m)/ω given in Eq. (8), www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 18 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 0.3 Measurement Measured a) a) Extrapolation Extrapolation KK KKmea mea KK KK ext ext 0.25 1.5 0.2 0.15 0.1 0.5 0.05 −0.05 −2 −1 0 1 2 3 −0.5 −2 −1 0 1 2 3 10 10 ω 10 10 ωh 10 10 ω l 10 10 l 10 10 ω 10 10 h ω 0.2 1 Measurement b) Measured Extrapolation b) Extrapolation KK mea KKmea KK ext KK ext 0.15 0.1 0.5 0.05 −0.05 −0.1 −2 −1 0 1 2 3 ωl h −0.5 10 10 10 10 10 10 −2 −1 0 1 2 3 10 10 ωl 10 10 10 10 Fig. 2. Linear susceptibility of intensive momentum of the system Fig. 1. Linear susceptibility of intensive energy of the system E/N M/N with respect to the global perturbation with X = 1 ∀i. The with respect to the global perturbation with X = 1 ∀i. The real and real and the imaginary parts are depicted in (a) and (b), respectively. the imaginary parts are depicted in (a) and (b), respectively. The The measured and extrapolated values are shown in red and black measured and extrapolated values are shown in red and black lines, lines, respectively. The result of the Kramers-Kronig inversion done respectively. The result of the Kramers-Kronig inversion done with with the measured and with with the extrapolated data are shown in the measured and with with the extrapolated data are shown in blue blue and magenta lines, respectively. and magenta lines, respectively. totically for high frequencies as 1/ω and 1/ω, respectively. whereas for low frequencies the real susceptibility tends to a We apply the truncated KK relations to the measured data very high value, this suggesting that the strongest response is to test the quality of the data inversion process. The estimates obtained for static perturbations. of the imaginary part (starting from the measured data of the (1) The measured real and imaginary parts of χ = real part) and of the real part (starting from the measured data μ,N (1) of the imaginary part) obtained by applying Eqs. (9–10) are 1/Nχ , where μ =M/N , are depicted in black in Fig. 2a, M,N (1) (1) shown for χ in blue in Fig. 1a, b and for χ in Fig. 2a, b. Interestingly, the spectral feature of the imaginary part is ε,N μ,N b. We observe that whereas agreement is very good for the shifted to higher frequencies than in the case of the energy real part for both susceptibilities, only a qualitative match is susceptibility, so that a well-distinct peak centered on value obtained for the imaginary part, with large discrepancies for of ω ≈ 6, which approximately corresponds to the natural ω. 2. In this latter case, moreover, the well-known problem time scale of the system. For low frequencies, the suscep- of KK divergence at the boundaries of integration (Lucarini tibility has almost exclusively a real component. As opposed et al., 2003, 2005) is very serious for ω =ω . to the previous case, the largest value for the in-phase re- l sponse is not obtained for ultralow frequencies, but rather for It is crucial to test whether the discrepancies are due to the ω ≈ 4. The asymptotic behavior of both the real and imagi- finiteness of the spectral range or are, instead, due to basic nary parts is in perfect agreement with the theoretical result problems in the applicability of the Ruelle response theory, given in Eq. (32), so that they are found to decrease asymp- related to the fact that the invariant probability measure of the Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ Im[χ] Re[χ] Im[χ] Re[χ] V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 19 unperturbed system actually features large deviations from They correspond, by definition, to the static susceptibility an SRB measure. of the observables e and m, respectively, for the global We proceed testing the first case. In order to widen the perturbation with X = 1 considered here. When evaluat- spectral range over which the susceptibility is defined, we ing the derivatives of e(F) and m(F) for F = 8 we obtain will exploit the asymptotic properties obtained in Sect. 3.2 (de(F)/dF) ≈ 1.6 and (dm(F)/dF) ≈ 0.11. These F=8 F=8 as well as the low frequency behavior of the susceptibility values are in good agreement with what found by extrapo- discussed in Sect. 2. We redefine the the imaginary part of lating the corresponding real part of the susceptibilities for (1) ω → 0 via KK relations and shown in Figs. 1a and 2a. the susceptibility of χ as follows ε,N Apart from the verification of the validity of KK relations, (1) we want to provide further support for the quality of the lin- Im[χ (ω )], 0 ≤ω ≤ω , l l ω ε,N (1) ear susceptibilities considered. (1) Im[χ (ω)] = (60) Im[χ (ω)], ω ≤ω ≤ω , ε,N l h ε,N First, we test the sum rules given in Eqs. (42) and (44) for , ω ≥ω , (1) the real part of the extrapolated susceptibilities χ (ω) and ε,N (1) where the measured data are sandwiched between the low χ (ω), respectively. Our findings are presented in Fig. 3, μ,N and high frequency limit. Whereas we have a rigorous result where it is shown that an excellent agreement (within 1%) is for the high frequency limit, the low frequency limit is com- found between the theoretical values and the numerical re- (1) puted by making the reasonable assumption that the lead- sults. Since Re[χ ](ω) is negative in the high-frequency e,,a ing order of the ω → 0 limit is linear (see discussion after range, the convergence of the integral to the theoretical value (1) Eq. 11). Similarly, the real part of the susceptibility χ can of the sum rule is from above, whereas the opposite occurs ε,N (1) be redefined as follows: for Re[χ (ω)]. Extending the integral for even larger val- m,,a ues of ω would bring the numerical results to an almost per- (1) Re[χ (ω )], 0 ≤ω ≤ω , l l ε,N fect agreement with the theory. (1) (1) Re[χ (ω)] = (61) Re[χ (ω)], ω ≤ω ≤ω , l h Following the definition given in Eq. (3), the Green func- ε,N ε,N F−2m (1) − , ω ≥ω , h tion G (τ) computed for an observable 8 and a given pat- ω 8 tern of perturbation flow X (x) (in this case X = 1∀i) can i i where we have used the fact that at low frequencies the real be used to compute the time-dependent linearized impact of part of the susceptibility is constant in ω up to a quadratic all perturbations with the same spatial pattern X (x) but with term. A corresponding procedure is used to extend the spec- arbitrary time modulation. Whereas the direct estimate of the (1) tral range of the χ (ω), where the suitable asymptotic be- Green function from the time dependent dynamics can be ob- μ,N haviors described in Sect. 3.2.1 are adopted. The red lines tained by performing an ensemble of simulations where the in Figs. 1a, b–2a, b present the results of such extrapola- time modulation of the perturbation is given by a δ(t) pattern tions, and the magenta lines show the outcome of the data (see discussion in Sect. 2, we take the indirect route by con- inversion of these functions performed via KK relations. We sidering Eq. (5). By applying the inverse Fourier Transform, (1) observe that the agreement is outstanding, with almost per- we derive the Green functions corresponding to χ (ω) and ε,a fect overlap inside the region where measurement is per- (1) χ (ω). The results are presented in Fig. 4: for both observ- μ,a formed and remarkable agreement also in the low and high ables the Green functions are clearly causal, and their short- frequency range. This is a very convincing evidence that time behavior agrees remarkably well with what be deduced the Ruelle response theory can be successfully applied for by looking at the asymptotic properties of the corresponding this system. Since the KK relations provide, first and fore- susceptibilities (compare Eqs. 41 and 43). most, consistency tests, the agreement the original and the KK-transformed susceptibility automatically confirms that 4.3 Local perturbation the extrapolation procedure we have adopted is correct. A still better agreement would be found had we taken into ac- The data obtained for the numerical simulations of the re- count value of ω larger than what considered in the extrapo- sponse to the local perturbation are, given the much weaker lation used here (up to 100 π ). overall strength of the forcing, much noisier that those pre- Furthermore, let’s consider the results presented in sented in the previous section. Nevertheless, we shall see Sect. 3.1. The slopes of the functions e(F) and m(F) are that all the theoretical predictions are verified to a surpris- given by ingly good degree of approximation. The global observables E and M are of the form dm(F) γ −1 N 8(x ,...,x ) = φ(x ), where φ(x) =x /2 and φ(x) = =λγF , (62) 1 n i i=1 dF x, respectively. We have consistently verified that the iden- (1) (1) (1) tity χ = χ = χ discussed in the previous section φ,1 φ/N,N φ,N applies in the whole spectral range explored by our simula- de(F) (1+γ) (1+γ) =λ F =m(F) . (63) tions, compatibly with the (slightly) different signal-to-noise dF 2 2 www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 20 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 1.8 (1) (1) Re[χ ] ε,N ε,N (1) (1) 1.6 Re[χ ] χ E,1 μ,N (1) π Im[χ ] ε,N 1.4 (1) Im[χ ] E,1 1.2 0.8 0.6 0.4 0.2 ω ω l h −0.2 −2 −1 0 1 2 3 −2 −1 0 1 2 3 10 10 10 10 10 10 10 10 10 10 10 10 max Fig. 5. Comparison between the linear susceptibility of the intensive Figure 5: Comparison between the linear susceptibility of the intensive energy for the global Fig. 3. Sum Rules for the real part of the linear susceptibility of Figure 3: Sum Rules for the real part of the linear susceptibility of intensive energy E/N energy for the global perturbation and of the total energy for the perturbation and of the total energy for the local perturbation. The signals are the same (black line) intensi and vmomen e enertum gy EM/N /N(black (blue line) line) of and the momentum system withMresp /Nect (blue to the line) global per- except for local theperturbation. different level of The noise. signals are the same except for the different turbation with X = 1 ∀i. The theoretical values are indicated in the figure. of the system i with respect to the global perturbation with X = 1 ∀i. level of noise. The theoretical values are indicated in the figure. 0.8 Measurement a) Extrapolation (1) KKmea G (τ) 0.7 KKext Θ(τ)(m + (F − 2m)τ) (1) 2.5 Gμ (τ) Θ(τ)(1 − τ) 0.6 0.5 1.5 0.4 0.3 0.2 0.5 0.1 −0.5 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 −0.1 ωl ω −0.2 −2 −1 0 1 2 3 10 10 10 10 10 10 Fig. 4. Green functions describing the time-dependent response of Figure 4: Green functions describing the time-dependent response of the observable of in- the observable of intensive energy E/N (black line) and momentum tensive energy E/N (black line) and momentum M/N (blue line) with respect to the global 0.8 Measurement M/N (blue line) with respect to the global perturbation with X = 1 b) perturbation with X = 1 ∀i. For both observables, the short time behaviori of the Green Extrapolation KKmea 0.7 KK function∀iestimated . For both from observ the asymptotic ables, the b short ehavior time of the beha corresp vior of onding the Green susceptibilit func-y is shown ext in figure. tion estimated from the asymptotic behavior of the corresponding 0.6 susceptibility is shown in figure. 0.5 0.4 ratios in the two sets of simulations. See, e.g., Fig. 5 for the (1) (1) 0.3 comparison between the two susceptibilities χ and χ . E,1 ε,N 0.2 We then proceed to analyze more in detail the linear susceptibilities related to local observables. In Fig. 6 we 0.1 present our results concerning the real and imaginary part (1) of χ (ω). Analogously to what observed in the previous E ,1 −0.1 40 ω h subsection, we have that once the measured susceptibility l −0.2 −2 −1 0 1 2 3 is extrapolated using the theoretical results obtained via re- 10 10 10 10 10 10 sponse theory and KK relations, we have an excellent agree- ment between the original real and imaginary parts and those Fig. 6. Linear susceptibility of the energy at the grid point where the obtained using the KK inversion. The KK algorithm, instead, local perturbation is applied. The real and the imaginary parts are provides only a partially satisfying outcome when only data depicted in (a) and (b), respectively. The measured and extrapolated from the measured range are considered. Relatively discrep- values are shown in red and black lines, respectively. The result of the Kramers-Kronig inversion done with the measured and with ancies are found near the boundaries of the data range, with with the extrapolated data are shown in blue and magenta lines, an especially serious bias near ω for the imaginary part of respectively. the susceptibility. Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ max dωRe[χ] G(τ) 0 Re[χ], Im[χ] Im[χ] Re[χ] V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 21 (1) (1) 0.25 When comparing χ (ω) and χ (ω) (see Fig. 5), we Measured E ,1 E,1 a) j Extrapolation KK mea KK observe that for low frequency the response of the energy ext 0.2 at the grid point where the perturbation is applied accounts for about half of the response of the total energy, thus im- 0.15 plying that the remaining half is redistributed among the re- maining N − 1 grid points. The relevance of the response of 0.1 grid points other than the directly perturbed one also explains (1) (1) 0.05 why the peak of Imχ (ω) (and so of Imχ (ω)) than that E,1 ε,N (1) of Imχ (ω) – see the frequency range 2. ω . 4. Slower E ,1 perturbations allow other grid points x 6= x to respond ef- k j −0.05 fectively. (1) Instead, since the leading asymptotic order of χ (ω) E ,1 −0.1 −2 −1 0 1 2 3 10 10 10 10 10 10 (1) and χ (ω) is the same, at high frequencies the local en- E,1 ergy response accounts for most of the energy response of 0.25 Measured b) Extrapolated KK the whole system. In this case, the incoming perturbation is mea KK ext 0.2 so fast that the internal time scales of the system as bypassed, and mainly a local effect is observed. Nevertheless, the sec- 0.15 ond leading order of the asymptotic expansion of the two susceptibilities has opposite sign (see Eqs. 41 and 54), which 0.1 suggests that at any large but finite frequency the local energy response is only a good approximation to the response of the 0.05 total energy. The changeover between the two regimes oc- (1) curs around the frequency of the peak of Imχ (ω), which E ,1 corresponds to a perturbation with period close to 1. −0.05 (1) Thanks to the asymptotic equivalence between χ (ω) E ,1 (1) (1) −0.1 −2 −1 0 1 2 3 and χ (ω) (and χ (ω)), they must obey the same sum 10 10 10 10 10 10 E,1 ε,N for the real part of the susceptibility (see Eqs. 42–55), even if the two real parts, as discussed above, are rather different in Fig. 7. Linear susceptibility of the momentum at the grid point value in the low frequency range and even in sign in the high where the local perturbation is applied. The real and the imaginary frequency range. Figure 8 confirms that this rather counter- parts are depicted in (a) and (b), respectively. The measured and extrapolated values are shown in red and black lines, respectively. intuitive behavior is actually observed. Note also that sum The result of the Kramers-Kronig inversion done with the measured rules, resulting from an integration, are less sensitive to noise and with with the extrapolated data are shown in blue and magenta in the data, but this occurs if and only if the underlying sig- (1) lines, respectively. nal is correct. Therefore, we understand that in χ (ω) and E,1 (1) χ (ω) the strong static and quasi-static response and the ε,N The investigation of the linear susceptibility of the vari- (rather odd) negative sign for high frequencies of the real part able x is not as insightful as that of E . We find that lin- j j of the linear susceptibility, which are crucially related to the (1) (1) ear susceptibility χ (ω) is quite similar to χ (ω) (and behavior for the grid points different from the perturbed one) x ,1 M,1 (1) compensate each other to guarantee agreement with the sum χ (ω), see Fig. 2) in both the real and imaginary parts at μ,N (1) rule obtained from the real part of χ (ω), which instead all frequency. The only notable differences are that the static E ,1 has a smaller range and more regular (monotonic) behavior response x is slightly larger than than of M , and that the with frequency. imaginary features a secondary peak at slightly larger fre- A formally similar – and analogously spectacular – spec- quencies than the main spectral feature. We have verified,as tral compensation has been observed in a physical process as in the previous cases, the results of the numerical simula- different from what we are analyzing here as the electromag- tions accurately agree with the theoretical results regarding netically induced transparency (Cataliotti et al., 1997). The the asymptotic behavior of both the real and imaginary part result obtained here supports previous findings obtained on and that KK relations map to high degree of precision the quasi-equilibrium systems suggesting that sum rules do not real and the imaginary parts into each other. See Fig. 7 for depend on many-particle interactions (Lucarini et al., 2003, details. 2005). We present as main finding of the analysis of the observ- able x that, as predicted by the theory, the real part of (1) (1) χ (ω) obeys the same sum rule as the real part of χ (ω) x ,1 M,1 www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 Im[χ] Re[χ] 22 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings The first scenario envisioned here pertains to the pair of (1) Ej,1 (1) (1) (1) (1) χ , χ ε,N E,1 linear susceptibilities χ (ω) and χ (ω), whereas the (1) π x ,1 M,1 m j x ,1 2 (1) (1) χ , χ second scenario is related to the pair of linear susceptibili- μ,N M,1 (1) (1) ties χ (ω) and χ (ω). Note that, taking into account the E ,1 E,1 asymptotic properties of the susceptibility of the observable 2 2 2 2 E = 1/2x + 1/2x + 1/2x + 1/2x , discussed in loc j j+1 j+2 j−1 (1) (1) the previous section we conclude that χ and χ should E ,1 E,1 loc be similar for all values of ω. Obviously, a similar argument applies if the leading order is real. This discussion further clarifies that the higher the −2 −1 0 1 2 3 10 10 10 10 10 10 max number of independent KK relations and related sum rules verified by a susceptibility functions, the more stringent are Fig. 8. Sum rules of the real part of the linear susceptibilities in- the constraints on its properties. Figure 8: Sum rules of the real part of the linear susceptibilities indicated in the legend. The dicated in the legend. The theoretical values are indicated in the theoretical values are indicated in the figure. 4.5 Additional properties of the linear susceptibility figure. The special mathematical properties of the linear susceptibil- (1) or of χ (ω), because the corresponding imaginary parts ities allow to investigate further properties of the response. μ,N (1) feature the same asymptotic behavior. Figure 8 shows that in In particular, we note that for m ≥ 1 the function [χ ] is (1) the case of the momentum variables the cumulative integral analytic in the upper complex ω-plane just as as χ . This is rather similar for the susceptibility of the local and of the allows, as discussed in Lucarini et al. (2005) to derive, in global variable, with small discrepancies in the region around principle, an infinite set of integral relations (KK and sum the peak of the response. rules) deriving just from the holomorphic proprieties of the susceptibility. As an example, we have considered the square 4.4 Further implications of Kramers-Kronig relations (1) of the linear susceptibility [χ (ω)] . From Eq. (41), it is ε,N and sum rules easy to prove that the following asymptotic expansion holds for large values of ω: We now show how the knowledge of the asymptotic behav- ior of the real and imaginary part and the knowledge of the 2 m m(F − 2m) (1) 2 −4 [χ (ω)] = − + +o(ω ). (64) validity of the KK relations and related sum rules allow to ε,N 2 3 ω ω draw general conclusions on the similarities and differences As shown in panel (a) of Fig. 9, KK relations are found between two given linear susceptibility functions. Let’s con- to connect up to a high degree of approximation the real sider the case that these two susceptibilities feature the same (1) and imaginary part of [χ (ω)] . Moreover, thanks to the first order asymptotic expansion in the high frequency limit. ε,N asymptotic behavior given in Eq. (64), it is possible to estab- Let’s assume that it is an odd power of ω, so that the real part lish, thanks to Eqs. (13–14), the following sum rules: is negligible for high frequencies. Therefore, the two suscep- th tibilities will obey the same sum rule for, e.g. the 0 moment (1) dνRe[χ (ν)] = 0, (65) of the real part. ε,N If they agree also in the asymptotic behavior of the real part, they cannot feature large discrepancies in the low fre- (1) 2 2 quency range of the real part of the susceptibility either, or dννIm[χ (ν)] = m , (66) ε,N otherwise the agreement of the sum rules would be broken. Therefore, the real part of the two susceptibilities are similar, Panel (b) of Fig. 9 shows that the obtained numerical results and, as a consequence of the KK relations, the two imaginary are in excellent agreement with the theoretical predictions. parts will also be similar. Note that these results do not have an obvious physical in- (1) If, instead, there is a discrepancy in the asymptotic behav- 2 terpretation, as the inverse Fourier Transform of [χ (ω)] ε,N ior of the real part of the two susceptibilities, the two real is given by the convolution product of the Green function parts will necessarily be rather different in the low frequency (1) G (t) with itself, while they depend only on the formal ε,N range, again in order to comply with the sum rule constraint. properties of the linear susceptibility. As the two real parts are different, the imaginary part of the two susceptibility will also be rather different, except, from hypothesis, in the high-frequency range. 5 Practical implications for climate change studies In this paper we have constructed and verified to a high degree of accuracy the linear response theory for a simple Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ max dωRe[χ] 0 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 23 From Eq. (2), we have that (1) 2 a) Re[(χ ) ] ε,N (1) Im[(χ ) ] ε,N (1) 2 Z 2.5 Re[(χ ) ]KK ε,N +∞ (1) (1) Im[(χ ) ]KK (1) ε,N hT i (t) = dσ G (σ )f(t −σ ). (67) S 1 1 1 −∞ 1.5 In practical terms, the left hand side of this equation is noth- ing but the ensemble average of the time series of the change between the globally averaged surface temperature of the planet at a time t after the perturbation has started. Note that 0.5 (1) the direct estimate of G (σ) is likely to be overwhelmingly difficult. Using Eq. (4), we have that: −0.5 (1) (1) hT i (ω) =χ (ω)f(ω), (68) −1 −2 −1 0 1 2 3 10 10 10 10 10 10 which implies that once we compute the Fourier Transform of the time series mentioned above and we know the modu- SR0 b) SR lating function f(t) (and so its Fourier Transform f(ω)), we (1) can reconstruct χ (ω). Let’s select a particularly simple π 2 example of modulating function f(t) = (2(t)−2(t −τ)). This is just a rectangular function of width τ , of height , and shifted from the origin by a forward time translation τ/2. In practical terms, this corresponds to changing abruptly the field CO concentration by  at time t = 0 and taking it back to its original value at t =τ . we then obtain: (1) (1) hT i (ω) hT i (ω) S S (1) χ (ω) = =ω . (69) f(ω) (sin(ωτ)+i(1− cos(ωτ))) (1) Once we know χ (ω), as widely discussed in this paper, 0 T −2 −1 0 1 2 3 S 10 10 10 10 10 10 ω (1) max we can compute G (t), and we know everything about the (1) Fig. 9. Properties of the square of the linear susceptibility χ . response of the system at all time scales, including the static ε,N (1) 2 response. Note that any choice of f(t) is equally valid to The real and imaginary parts of [χ ] with their KK transforms ε,N set up this procedure as long as f(t) is square integrable. are depicted in the panel (a), the vanishing sum rule for the real part This implies that, in a very profound way, the kind of forcing and the non-vanishing sum rule for the imaginary part are depicted in panel (b). scenarios used in the various assessment Reports of IPCC, where the CO concentration typically stabilizes at a differ- ent value from the preindustrial one (so that f(t) does not yet prototypical climate model by computing the frequency- tend to 0 as t → ∞) are not necessarily the only nor the best dependent susceptibilities of several relevant observables re- ones, in spite of what could be intuitively guessed, to study lated to localized and global patterns of forcings. These re- even the steady state response of the system. sults pave the way for devising a rigorous methodology to Obviously, a similar set of experiments could be devised be used by climate models of any degree of complexity for for studying rather thoroughly the response of the climate studying climate change at, in principle, all time scales us- system to a variety of forcings, such as changes in the O ing only a very limited set of experiments, and for exploiting concentration, aerosols, solar radiance, as well as to changes effectively the currently adopted ensemble runs methods. in the parameterizations. In the case of uncoupled models of Let’s consider, for sake of simplicity, that the observable one subdomain of the climate system (e.g. atmospheric and 8 is the time-dependent globally averaged surface temper- oceanic GCMs, land-surface models), this strategy could be ature of the planet T , that F(x) represents the whole set used to study the impact of perturbations to the boundary of climate equations in a baseline scenario (e.g., with pre- conditions provided by the other subdomains of the climate industrial CO concentration), and that the perturbation field system. X(x) is nothing but a constant field of CO concentration, which directly impacts only the radiative part of the code. The perturbation is modulated by a time-dependent function f(t) to be specified below. We assume, for simplicity, that the model does not feature daily or seasonal variations in the radiative input at the top of the atmosphere. www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 24 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 6 Summary, discussion and conclusions stochastic dynamics using the results of the mathematical theory of random dynamical systems is emerging as new, The climate can be seen as a complex, non-equilibrium promising paradigm for the investigation of the structural system, which generates entropy by irreversible processes, properties of the climate system (Chekroum et al., 2011). transforms moist static energy into mechanical energy We have proposed a different perspective. In agreement (Lorenz, 1967; Peixoto and Oort, 1992) as if it were a heat with the view given above, we have taken as mathematical engine (Johnson, 2000; Lucarini, 2009a), and, when the ex- framework for the analysis of the climate system that of non- ternal and internal parameters have fixed values, achieves equilibrium statistical mechanics, and have focused on the a steady state by balancing the input and output of energy steady state properties of ergodic dynamical systems (Eck- and entropy with the surrounding environment (Ozawa et al., mann and Ruelle, 1985) possessing the special property of 2003; Lucarini, 2009a). For such basic reasons, the tool of having an invariant measure of the SRB type (Ruelle, 1989). equilibrium and quasi-equilibrium statistical mechanics can- As proposed by the chaotic hypothesis (Gallavotti, 1996; not provide suitable tools for studying the fundamental prop- Gallavotti and Cohen, 1995), this mathematical framework erties of the climate system. In particular, the fluctuation- is well suited for analyzing general non-equilibrium physical dissipation theorem, which allows for deriving the properties systems. of the response of the system to external perturbations from In this context, the impact on the system of general per- the observations of its internal variability cannot be applied. turbation can be treated using the response theory recently It is reasonable to ask whether is possible to evaluate how introduced by Ruelle (1998a,b, 2009), which allows to com- far from equilibrium the climate system actually is. It is pos- pute the change in the expectation value of a generic observ- sible to evaluate such distance in a mathematically sound able as a perturbative series where each term is given by the way by assessing the ratio of the dimensionality of the at- average over the unperturbed invariant measure of a function tractor of the system over the total number of degrees of of the phase space which depends on the considered observ- freedom. Whereas a ratio close to one indicates that only able and on the applied perturbation. In other terms, even if small deviations from equilibrium are present, a small ratio the internal dynamics of the system is nonlinear and chaotic, suggests that strongly non-equilibrium conditions are estab- the leading order of the response is in general linear with the lished. See Posch and Hoover (1998) for a detailed treatment strength of the added perturbation. This approach overcomes of this problem in the classical case of heat conduction. In the difficulties related to the singularity of the invariant mea- the case of a quasi-geostrophic atmospheric model forced by sure discussed in Thuburn (2005). Earth-like boundary conditions, the dimensionality of the at- At each order, the propagator of the perturbation, i.e. the tractor of the model is about one order of magnitude smaller Green function, is causal. This allows for applying disper- than the total number of degrees of freedom (Vannitsem and sion theory and establish general integral constraints – KK Nicolis, 1997). While not conclusive, this seems to suggest relations – connecting the real and imaginary parts of the sus- that the best framework to interpret the climate is that of a far ceptibility, i.e. the Fourier Transform of the Green function from equilibrium system. (Lucarini, 2008b, 2009b) Following either explicitly or implicitly the programme of In this paper we have first recapitulated the main as- the Catastrophe theory (Arnold, 1992), many authors have pects of the general response theory and have propose some approached the problem of understanding the fundamental new general results, which boil down to consistency rela- properties of the climate system by looking at the detailed tions between the linear susceptibilities of different observ- structure of the bifurcations of the deterministic dynamical ables. The obtained equation provides the basic idea of why system constructed heuristically in order to represent the dy- the fluctuation-dissipation theorem does not apply in non- namics of the main climate modes using as few degrees of equilibrium cases. freedom as possible. Such an approach often hardly allows to efficiently represent the fluctuations and the statistical prop- We have showed for the first time that the Ruelle linear erties of the system. The introduction of stochastic forcing response theory can be applied with great success to ana- provides a relatively simple but conceptually rich partial so- lyze the climatic response to general forcings. We have cho- lution to some of these draw-backs, even if the Hasselmann sen as test bed the L96 model (Lorenz, 1996; Lorenz and programme (Hasselmann, 1976) suffers from the need for a Emanuel, 1998; Lorenz, 2004), which, in spite of its simplic- – usually beyond reach – closure theory for the properties ity, has a well-recognized prototypical value as it is a spa- of noise. Therefore, the stochastic component is usually in- tially extended one-dimensional model and presents the ba- troduced ad hoc, with the ensuing lack of universality and/or sic ingredients, such as dissipation, advection and the pres- robustness when various levels of truncations are considered. ence of an external forcing, of the actual atmosphere. Such These strategies have anyway brought to outstanding scien- a model features a different level of complexity with respect tific results and has been suggested the existence of generic to those adopted in previous numerical investigations of Ru- mathematical structures present in hierarchies of CMs (Saltz- elle’s theory (Reick, 2002; Cessac and Sepulchre, 2007; Lu- man, 2002). Recently, the unified treatment of chaotic and carini, 2009b), and has been already tested in terms of linear Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 25 response (albeit with a different methodology and in a differ- We believe that the proposed approach, which we may ent theoretical context) in Abramov and Majda (2007). dub as spectroscopy of the climate system, may constitute a We have analyzed the frequency dependence of the re- mathematically rigorous and practically very effective way sponse of the local and global energy and momentum of the to approach the problem of evaluating climate sensitivity system to perturbations having a global spatial pattern and to and climate change from a radically new perspective. In perturbations acting only on one grid point. We have derived this regard, we have proposed a rigorous way to compute, analytically several properties of the corresponding suscep- e.g., the surface temperature response to changes in the the tibilities, such as asymptotic behavior, validity of KK rela- CO concentration at all time scales using only a specific tions, and sum rules. We have shown that all the coefficients set of simulations, and taking advantage of the theoretical of the leading asymptotic expansions as well as the integral results presented here. Given the ever-increasing interest constraints can be written as linear functions of parameters towards decadal and multidecadal climate prediction, these that describe unperturbed properties of the system, and in tools could be of relevant practical interest and their applica- particular its average energy and average momentum. The bility could benefit from technological platform aimed at cre- theory has been used to explain differences in the response of ating ensemble simulations comprising of many members. local and global observables, in defining the intensive proper- We underline that our approach takes into account all the ties of the system and in generalizing the concept of climate (linear and nonlinear) feedbacks of the system, as they are sensitivity to all time scales. included in the definition of the Green function. This, at a We have then verified the theoretical predictions from the very practical level, is the great advantage of using Ruelle’s outputs of the simulations up to a high degree of precision, formulas. even if we have used rather modest computational resources At a more basic level, whereas considering more complex (a total of about 30 cpu days of a mid-range commercial lap- models requires heavier computational resources, the modest top). We have verified that the linear response theory holds cost of the present set of simulations suggests that, at least for for perturbations of intensity accounting to up to about 10% global or regional climatic observables, it is feasible to test of the unperturbed forcing terms. Even when local pertur- the theory discussed here for simplified yet Earth-like climate bation and local observables are considered it is possible models without resorting to top-notch computing facilities. to achieve a signal-to-noise ratio which permits rather sat- Moreover, while in this paper we have computed the sus- isfactory comparisons with the theory. We have proved that ceptibilities using, on purpose, a very cumbersome method, the combined use of KK relations and the knowledge of the more efficient strategies can be devised, at least when the asymptotic behavior of the susceptibilities allows for extrap- linear regime of the response is considered. Apart from the olating in a rigorous way the observed data. We also have practical example given for the case of the impact of the CO shown how to reconstruct the linear Green function, which concentration, these include studying the response of the sys- can be used to map perturbations of general time modula- tem to δ(t) like perturbations, which gives directly the Green tion into changes in the expectation value of the considered function of the system, and including in the forcing various observable for finite as well as infinite time. monochromatic signals. Of course, in all cases, a Monte Our numerical experiments have been performed using Carlo approach is needed in order to sample effectively the one of the standard settings of the L96 model, namely the attractor of the unperturbed system in terms of the initial con- version identified by having N = 40 degrees of freedom and ditions of the simulations. forcing F = 8. Nevertheless, some newly obtained empirical These results pave the way for future investigations aimed closure equations expressing the average energy and the av- at improving and extending the theoretical framework pre- erage momentum of the unperturbed system as simple power sented here, at finding results of general applicability in the laws of F (with no dependence on N ) have allowed to ex- context of the modelling of geophysical fluid dynamics, and, tend our results to the entire class of chaotic L96 models. finally at answering specific questions of relevance for cli- The regular scaling of the properties of the system and of its mate dynamics. In this paper we have analyzed the sim- response with F agrees with what observed in Abramov and ple case of the linear response, but, as discussed in Ruelle Majda (2007). (1998b) and Lucarini (2008b, 2009b), we have the algorithm In this paper we have only used the KK relations in the to compute higher order terms, so that the treatment of the most simplistic framework, i.e., computing the KK trans- nonlinear response in entirely feasible. forms and evaluating their agreement with the original data. In the first direction, we foresee the possibility of writing Actually, several more sophisticated analysis techniques are out explicitly the linear susceptibility of a general observable available, such as recursive self-consistent algorithms, where by projecting the perturbation onto the unstable, neutral and the measured data are taken as first guess, exploiting the fact stable manifolds and analyzing separately the contributions that multiple applications of KK relations, combined with to the total response. This will probably require the adoption sum rules, automatically filter our the noise and remove most of adjoint techniques, and will benefit from the recent algo- of the spurious signal (Lucarini et al., 2005). rithms proposed in Abramov and Majda (2007) and Majda and Wang (2010). Moreover, we will be testing the radius of www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 26 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings convergence of the Ruelle response theory in some specific tives such as climateprediction.net, and the new project examples. PCMDI/CMIP5 (http://cmip-pcmdi.llnl.gov/cmip5/), which Along the second direction, we propose to study the im- will provide a crucial input for the Fifth Assessment Report pact of stochastic forcing to deterministic chaotic models by of the IPCC. treating the (additive or multiplicative) noise as a perturba- Acknowledgements. VL acknowledges the financial support of the tion to be analyzed using the linear and nonlinear Ruelle EU-ERC research grant NAMASTE and of the Walker Institute response theory and related spectral methods. We will try Research Development Fund. The authors acknowledge the fruitful to complement the results recently obtained in Majda and interactions with the participants to the programme “Mathematical Wang (2010) from an entirely different approach, i.e. bypass- and Statistical Approaches to Climate Modelling and Prediction” ing the Fokker-Planck equation formalism and using directly held at the Newton Institute for Mathematical Sciences (Cam- the first and higher order terms of the Ruelle response for- bridge, UK). mula. Moreover, we shall look into the spectral peaks of the susceptibilities and try to understand how the amplification Edited by: H. A. Dijkstra Reviewed by: two anonymous referees of the response is related to resonances of the system and to the activation of positive feedbacks. Since at tipping points the climate response to perturbations is expected to diverge, References the formalism discussed here could be used also for under- standing the processes leading to large-scale changeovers of Abramov, R. and Majda, A. J.: Blended response algorithms for the dynamics of climate. linear fluctuation-dissipation for complex nonlinear dynamical Along the third direction, we envision the analysis of the systems, Nonlinearity, 20, 2793–2821, 2007. impact of topography on the statistical properties of the cir- Arnold, V. I.: Catastrophe Theory, Springer, Berlin, 1992. culation in a quasi-geostrophic setting, thus extending in a Baladi, V.: On the susceptibility function of piecewise expanding interval maps, Commun. Math. Phys., 275, 839–859, 2007. climatic perspective what presented in Speranza et al. (1985). Brayshaw, D. J., Hoskins, B., and Blackburn, M.: The storm track Moreover, we will tackle in an idealized setting the problem response to idealised SST perturbations in an aquaplanet GCM, of computing the response of the storm track to changes in J. Atmos. Sci., 65, 2842–2860, 2008. the surface temperature (Brayshaw et al., 2008). Positive re- Cacuci, D. G.: Sensitivity theory for nonliner systems. Part I. Non- sults in this direction, albeit using a different formalism, have linear functional analysis approach, J. Math. Phys., 22, 2794– been shown in Gritsun et al. (2008). Moreover, we will try to 2802, 1981a. compute along the way discussed here the climate response Cacuci, D. G.: Sensitivity theory for nonliner systems. Part II. Ex- to changes in CO and solar irradiance using simplified but tensions to additional classes of responses, J. Math. Phys., 22, rather valuable climate models like PLASIM (Fraedrich et 2803–2812, 1981b. al., 2005). Cataliotti, F. S., Fort, C., Hansch, ¨ T. W., Inguscio, M., and The results of the linear response theory should be com- Prevedelli, M.: Electromagnetically induced transparency in cold free atoms: Test of a sum rule for nonlinear optics, Phys. Rev. A, pared to those obtained from a simplified application of the 56, 2221–2224, 1997. fluctuation-dissipation theorem. The difference between the Cessac, B. and Sepulchre, J.-A.: Linear response, susceptibility and two results would inform us on the relevance of the ef- resonances in chaotic toy models, Physica D, 225, 13–28, 2007. fect of the external perturbation along the stable direction Chekroun, M. D., Simonnet, E., and Ghil, M.: Stochastic climate of the unperturbed flow, which cannot be precisely captured dynamics: Random attractors and time-dependent invariant mea- by the internal fluctuations of the system and, thus, can- sures, Physica D, in press, 2011. not be rigorously described in the context of the fluctuation- Dolgopyat, D.: On differentiability of SRB states for partially hy- dissipation theorem. Moreover, this could also suggest ways perbolic systems, Invent. Math., 155, 389–449, 2004. to improve the empirical evaluation and parameterisation Eckmann, J.-P. and Ruelle, D.: Ergodic theory of choas and strenge of the corresponding contribution to the response. This attractors, Rev. Mod. Phys., 57, 617–655, 1985. would also be of relevance to understand why, in spite of Errico, R.: What Is an Adjoint Model?, B. Am. Meteor. Soc., 78, 2577–2591, 1997. all the involved approximations and the theoretically unjusti- Eyink, G. L., Haine, T. W. N., and Lea, D. J.: Ruelle’s linear re- fiable assumptions, previous applications of the equilibrium sponse formula, ensemble adjoint schemes and Lvy flights, Non- fluctuation-dissipation theorem (often in a hyper-simplified linearity, 17, 1867–1889, 2004, quasi-gaussian setting) have been relatively successful in Fertig, E., Harlim, J., and Hunt, B.: A comparative study of 4D-Var predicting the climate response to external forcings. and 4D ensemble Kalman filter: perfect model simulations with Finally, we would like to remark that the theory Lorenz-96, Tellus A, 59, 96–101, 2007. and the practical recipes proposed here could be of di- Fraedrich, K., Jansen, H., Kirk, E., Luksch, U., and Lunkeit, F.: The rect interest for all projects aimed at auditing climate Planet Simulator: Towards a user friendly model, Meteorol. Z., models’ performances and at studying practical problems 14, 299–304, 2005. connected to climate change, such as PCMDI/CMIP3 Frye, G. and Warnock, R. L.: Analysis of partial wave dispersion (http://www-pcmdi.llnl.gov/), distributed computing initia- relations, Phys. Rev., 130, 478–494, 1963. Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/ V. Lucarini and S. Sarno: Computation of the climatic response to general forcings 27 Gallavotti, G.: Chaotic hypotesis: Onsanger reciprocity and Leith, C. E.: Climate response and fluctuation dissipation, J. Atmos. fluctuation-dissipation theorem, J. Stat. Phys., 84, 899–926, Sci., 32, 2022–2026, 1975. 1996. Lorenz, E. N.: Deterministic nonperiodic flow, J. Atmos. Sci., 20, Gallavotti, G.: Nonequilibrium statistical mechanics (stationary): 130–141, 1963. overview, in: Encyclopedia of Mathematical Physics, edited by: Lorenz, E. N.: The Nature and Theory of the General Circulation of Franc ¸oise, J.-P., Naber, G. L., Tsou Sheung Tsun, Elsevier, Am- the Atmosphere, WMO, Geneva, 1967. sterdam, 530–539, 2006. Lorenz, E.: Forced and free variations of weather and climate, J. Gallavotti, G. and Cohen, E. D. G.: Dynamical ensembles in sta- Atmos. Sci., 36, 1367–1376, 1979. tionary states, J. Stat. Phys., 80, 931–970, 1995. Lorenz, E. N.: Predictability-A problem partly solved, Proc. Sem- Giering, R. and Kaminski, T.: Recipes for Adjoint Code Construc- inar on Predictability, Vol. 1, Reading, Berkshire, United King- tion, ACM Trans. On Math. Software, 24, 437–474, 1998. dom, ECMWF, 1–18, 1996. Ghil, M. and Malanotte-Rizzoli, P.: Data assimilation in meteorol- Lorenz, E. N.: Designing Chaotic Models, J. Atmos. Sci., 62, 1574– ogy and oceanography, Adv. Geophys., 33, 141–266, 1991. 1587, 2004. Gritsun, A. and Branstator G.: Climate response using a three- Lorenz, E. N. and Emanuel, K.: Optimal sites for supplementary dimensional operator based on the fluctuation-dissipation theo- weather observations, J. Atmos. Sci., 55, 399–414, 1998. rem, J. Atmos. Sci., 64, 2558–2575, 2007. Lucarini, V.: Validation of climate models, in: Encyclopedia of Gritsun A., Branstator G. and Majda A. J.: Climate response of Global Warming and Climate Change, edited by: Philander, G., linear and quadratic functionals using the fluctuation-dissipation SAGE, Thousand Oak, 1053–1057, 2008a. theorem, J. Atmos. Sci., 65, 2824–2841, 2008. Lucarini, V.: Response theory for equilibrium and non-equilibium Hall, M. C. G.: Application of adjoint sesnitivity theory to an atmo- statistical mechanics: causality and generalized Kramers- spheric general circulation model, J. Atmos. Sci., 42, 2644–2651, Kroning relations, J. Stat. Phys., 131, 543–558, 2008b. 1986. Lucarini, V.: Thermodynamic efficiency and entropy produc- Hall, M. C. G. and Cacuci, D. G.: Physical interpretation of the tion in the climate system, Phys. Rev. E, 80, 021118, adjoint functions for sensitivity analysis of atmospheric models, doi:10.1103/PhysRevE.80.021118, 2009a. J. Atmos. Sci., 40, 2537–2546, 1983. Lucarini, V.: Evidence of Dispersion Relatios for the Nonlinear Re- Hall, M. C. G., Cacuci, D. G., and Schlesinger, M. E.: Sensitivity sponse of Lorenz 63 System, J. Stat. Phys., 134, 381–400, 2009b. analysis of a radiative-convective model by the adjoint method, Lucarini, V., Bassani, F., Saarinen, J. J., and Peiponen, K.-E.: Dis- J. Atmos. Sci., 39, 2038–2050, 1982. persion theory and sum rules in linear and nonlinear optics, Riv- Hallerberg, S., Pazo, ` D., Lopez, ` J. M., and Rodr´ ıguez, ista Nuovo Cimento, 26, 1–120, 2003. M. A.: Logarithmic bred vectors in spatiotemporal Lucarini, V., Saarinen, J. J., Peiponen, K.-E., and Vartiaine, E. chaos: Structure and growth, Phys. Rev. E, 81, 066204, M.: Kramers-Kronig relations in Optical Materials Research, doi:10.1103/PhysRevE.81.066204, 2010. Springer, Heidelberg, 2005. Hasselmann, K.: Stochastic climate models, Part 1: Theory, Tellus, Lucarini, V., Speranza, A., and Vitolo, R.: Parametric smoothness 28, 473–485, 1976. and self-scaling of the statistical properties of a minimal climate Haven, K., Majda, A., and Abramov, R.: Quantifying predictability model: What beyond the mean field theories?, Physica D, 234, through information theory: small sample estimation in a non- 105–123, 2007. Gaussian framework, J. Comp. Phys., 206, 334–362, 2006. Marini Bettolo Marconi, U., Puglisi, A., Rondoni, L., and Vulpi- Held, I.: The gap between simulation and understanding in climate ani, A.: Fluctuation-dissipation: Response theory in statistical modeling, B. Am. Meteor. Soc., 86, 1609–1614, 2005. physics, Phys. Rep., 461, 111–195, 2008. Intergovernmental Panel on Climate Change 2007: The Physical Majda, A. J. and Wang, X.: Linear response theory for statistical en- Science Basis. Contribution of Working Group I to the Fourth sembles in complex systems with time-periodic forcing, Comm. Assessment Report of the Intergovernmental Panel on Climate Math. Sci., 8, 145–172, 2010. Change, Cambridge University Press, Cambridge, 2007. Nussenzveig, H. M.: Causality and Dispersion Relations, Academic Johnson, D. R.: Entropy, the Lorenz Energy Cycle and Climate, Press, New York, 1972. in: General Circulation Model Development: Past, Present and Orrell, D.: Model Error and predictability over Different Timescales Future, edited by: Randall, D. A., Academic Press, New York, in the Lorenz-96 System, J. Atmos. Sci., 60, 2219–2228, 2003. 659–720, 2000. Ozawa, H., Ohmura A., Lorenz R. and Pujol T.: The second law Kubo, R.: Statistical-mechanical theory of irreversible processes. I, of thermodynamics and the global climate system: a review of J. Phys. Soc. Jpn., 12, 570–586, 1957. the maximum entropy production principle, Rev. Geophys., 41, Kubo, R.: The fluctuation-dissipation theorem, Rep. Prog. Phys., 1018, doi:10.1029/2002RG000113, 2003. 29, 255–284, 1966. Peiponen, K.-E., Vartiainen, E. M., and Asakura, T.: Dispersion, Lacorata, G. and Vulpiani, A.: Fluctuation-Response Relation and Complex Analysis and Optical Spectroscopy, Springer, Heidel- modeling in systems with fast and slow dynamics, Nonlin. Pro- berg, 1999. cesses Geophys., 14, 681–694, doi:10.5194/npg-14-681-2007, Peixoto, J. and Oort, A.: Physics of Climate, Springer, New York, 2007. 1992. Lange, P. L. and Alexeev, V. A.: Estimating 2 × CO warming Penland, C.: Noise out of chaos and why it won’t go away, B. Am. in an aquaplanet GCM using the fluctuation-dissipation theorem, Meteor. Soc., 84, 921–925, 2003. Geophys. Res. Lett., 32, L23708, doi:10.1029/2005GL024136, Posch, H. A. and Hoover, W. G.: Heat conduction in one- 2005. dimensional chains and nonequilibrium Lyapunov spectrum, www.nonlin-processes-geophys.net/18/7/2011/ Nonlin. Processes Geophys., 18, 7–28, 2011 28 V. Lucarini and S. Sarno: Computation of the climatic response to general forcings Phys. Rev. E, 58, 4344–4350, 1998. Shukla, J., Hagedorn, R., Hoskins, B., Kinter, J., Marotzke, J., Reick, C. H.: Linear response of the Lorenz system, Phys. Rev. E, Miller, M., Palmer, T., and Slingo J.: Revolution in climate pre- 66, 036103, doi:10.1103/PhysRevE.66.036103, 2002. diction is both necessary and possible: A declaration at the world Ring, M. J. and Plumb, R. A.: The response of a simplified modelling summit for climate prediction, B. Am. Meteor. Soc., GCM to axisymmetric forcings: Applicability of the fluctuation- 90, 175–178, 2009. dissipation theorem, J. Atmos. Sci., 65, 3880–3898, 2008. Thuburn, J.: Climate sensitivities via a Fokker-Planck adjoint ap- Ruelle, D.: Chaotic Evolution and Strange Attractors, Cambridge proach, Q. J. Roy. Meteor. Soc. 131, 73–92, 2005. University Press, Cambridge, 1989. Trevisan, A. and Uboldi, F.: Assimilation of standard and targeted Ruelle, D.: General linear response formula in statistical mechanics, observations within the unstable subspace of the observation- and fluctuation-dissipation theorem far from equilibrium, Phys. analysis-forecast cycle, J. Atmos. Sci., 61, 103–113, 2004. Lett. A, 245, 220–224, 1998a. Trevisan, A., D’Isidoro, M., and Talagrand, O.: Four-dimensional Ruelle, D.: Nonequilibrium statistical mechanics near equilibrium: variational assimilation in the unstable subspace and the opti- computing higher order terms, Nonlinearity, 11, 5–18, 1998b. mal subspace dimension, Q. J. Roy. Meteor. Soc., 136, 487–496, Ruelle, D.: A review of linear response theory for general differen- 2010. tiable dynamical systems, Nonlinearity, 22, 855–870, 2009. Vannitsem, S. and Nicolis, C.: Lyapunov Vectors and Error Growth Saltzman, B.: Dynamic Paleoclimatology, Academic Press, New Patterns in a T21L3 Quasigeostrophic Model, J. Atmos. Sci., 54, York, 2002. 347–361, 1997. Speranza, A., Buzzi, A., Trevisan, A., and Malguzzi, P.: A The- Wilks, D. S.: Effects of stochastic parametrizations in the Lorenz ory of Deep Cyclogenesis in the Lee of the Alps. Part I: Mod- ’96 system, Q. J. Roy. Meteor. Soc., 131, 389–407, 2006. ifications of Baroclinic Instability by Localized Topography, J. Atmos. Sci, 42, 1521–1535, 1985. Nonlin. Processes Geophys., 18, 7–28, 2011 www.nonlin-processes-geophys.net/18/7/2011/

Journal

Nonlinear Processes in GeophysicsUnpaywall

Published: Jan 12, 2011

There are no references for this article.