## Abstract

Systematic strategies from applied mathematics for stochastic modelling in climate are reviewed here. One of the topics discussed is the stochastic modelling of mid-latitude low-frequency variability through a few teleconnection patterns, including the central role and physical mechanisms responsible for multiplicative noise. A new low-dimensional stochastic model is developed here, which mimics key features of atmospheric general circulation models, to test the fidelity of stochastic mode reduction procedures. The second topic discussed here is the systematic design of stochastic lattice models to capture irregular and highly intermittent features that are not resolved by a deterministic parametrization. A recent applied mathematics design principle for stochastic column modelling with intermittency is illustrated in an idealized setting for deep tropical convection; the practical effect of this stochastic model in both slowing down convectively coupled waves and increasing their fluctuations is presented here.

## 1. Introduction

Stochastic modelling for climate is important for understanding the intrinsic variability of dominant low-frequency teleconnection patterns in climate, to provide cheap low-dimensional computational models for the coupled atmosphere–ocean system and to reduce model error in standard deterministic computer models for extended-range prediction through appropriate stochastic noise (Palmer 2001).

The present contribution is a research-expository paper on systematic strategies for stochastic climate modelling from the perspective of modern applied mathematics. In the modern applied mathematics, ‘modus operandi’ rigorous mathematical analysis, qualitative, asymptotic and numerical modelling are all blended together in a multidisciplinary fashion to provide systematic guidelines to address real-world problems (Majda 2000). For stochastic modelling in climate, the modern applied mathematics tool kit includes stochastic differential equations and discontinuous Markov jump processes (Gardiner 1985), systematic asymptotic reduction techniques, nonlinear dynamical systems theory and ideas from both statistical physics (Majda & Wang 2006) and mathematical statistics (Kravtsov *et al*. 2005; Majda *et al*. 2006*a*); mathematical rigour provides unambiguous guidelines in idealized models. Another facet of the modern applied mathematics philosophy is the development of qualitative models that represent a Platonic ideal for central issues simultaneously in diverse scientific disciplines such as material science, biomolecular dynamics and climate science.

In §2, we illustrate and apply this modern applied mathematics philosophy to stochastic modelling of the low-frequency variability of the atmosphere. The systematic mathematical theory (Majda *et al*. 1999, 2001, 2002, 2003, 2006*b*; collectively referred to as MTV hereafter; Franzke *et al*. 2005; Franzke & Majda 2006) for these problems is briefly reviewed including the central role and physical mechanisms responsible for multiplicative noise in the low-frequency dynamics. In this context, the Platonic ideal from applied mathematics is the truncated Burgers–Hopf model (Majda & Timofeyev 2000). A new simplified low-dimensional stochastic model that reproduces key features of atmospheric general circulation models (GCMs) is used there to test the fidelity of stochastic mode reduction techniques. A recent diagnostic statistical test with firm mathematical underpinning for understanding and interpreting the dynamical sources of the small departures from Gaussianity in low-frequency variables (Franzke *et al*. 2007) is also developed there.

While §2 deals with applied mathematical modelling through stochastic differential equations and §3 is devoted to the systematic development of stochastic lattice models to capture unresolved features that are highly intermittent in space and time such as deep convective clouds, cloud cover in subtropical boundary layers, sub-mesoscale eddies in the ocean and mesoscale sea ice cover. Here the mathematical tools involve a family of discontinuous Markov jump processes with multi-scale behaviour in space–time called stochastic spin-flip models. The key mathematical development involves systematic strategies to coarse grain such stochastic spin-flip models to achieve computational efficiency while retaining crucial features of the microscale interactions (Katsoulakis & Vlachos 2003; Katsoulakis *et al*. 2003*a*,*b*). The use of such stochastic lattice models to parametrize key features of tropical convection is briefly reviewed (Majda & Khouider 2002; Khouider *et al*. 2003). For the coupling of continuum models like a GCM to a stochastic lattice model as well as in many diverse applications, an applied mathematics Platonic ideal model has recently been introduced and analysed by Katsoulakis *et al*. (2004, 2005, 2006, 2007; hereafter KMS). This model consists of a system of ordinary differential equations (ODEs) for continuum variables ** X**,(1.1)two-way coupled to a stochastic spin-flip model written abstractly here as(1.2)where denotes the spatial coverage;

*L*is the generator;

*f*is a test function; and denotes the expected value. This idealized class of models has been used to systematically analyse the effects of various coarse-graining procedures on processes with intermittency, large-scale bifurcations and microscale phase transitions (KMS 2004, 2005, 2006, 2007). A concrete example for tropical convection in climate is given in §3. A new application of these stochastic lattice models to capture intermittent features and improve the fidelity of deterministic parametrizations of convection with clear deficiencies is developed in §3. First, the systematic design principles for (1.1) and (1.2) (KMS 2006, 2007) are used to calibrate a stochastic column model for tropical convection with intermittency and then the new results are presented on the practical effect of slowing down convectively coupled waves and increasing their fluctuations through the stochastic lattice models.

## 2. Systematic low-dimensional stochastic mode reduction and atmospheric low-frequency variability

A remarkable fact of Northern Hemisphere low-frequency variability is that it can be efficiently described by only a few teleconnection patterns that explain most of the total variance (e.g. Wallace & Gutzler 1981). These few teleconnection patterns not only exert a strong influence on regional climate and weather but are also related to climate change (Hurrell 1995). These properties of teleconnection patterns make them an attractive choice as basis functions for climate models with a highly reduced number of d.f. The development of such reduced climate models involves the solution of two major issues: (i) how to properly account for the unresolved modes that are also known as the closure problem, and (ii) how to define a small set of basis functions that optimally represent the dynamics of the major teleconnection patterns. This section addresses primarily the first issue and presents a rigorous strategy of how to systematically account for the unresolved d.f.

The simplest approach to derive highly truncated models of teleconnection patterns is to empirically fit simple stochastic models (e.g. autoregressive models and fractionally differenced models) to individual scalar teleconnection indices (Feldstein 2000; Stephenson *et al*. 2000; Percival *et al*. 2001). Statistical tests usually cannot distinguish if short- or long-memory models provide the better fit. A more complex approach, which also tries to capture deterministic interactions between different teleconnection patterns, is to linearize the equations of motion around a climatological mean state. Such models can be determined empirically from data or by using the linearized equations of motion. These models can be forced either by a random forcing (Branstator 1990; Newman *et al*. 1997; Whitaker & Sardeshmukh 1998; Zhang & Held 1999) or by an external forcing representing tropical heating (Branstator & Haupt 1998). To ensure stability of these linear models, damping is added according to various ad hoc principles. There is a recent survey of such modelling strategies (Delsole 2004).

A more powerful method is to empirically fit nonlinear stochastic models with possibly multiplicative (state dependent) noise by using the Fokker–Planck equation (Gardiner 1985; Sura 2003; Berner 2005). To reliably estimate the drift and diffusion coefficients in the Fokker–Planck equation is a subtle inverse problem that requires very long time series, and is further complicated by the need to retain the leading-order eigenvalue structure of the Fokker–Planck operator in order to keep the autocorrelation time scales of the original model (Crommelin & Vanden-Eijnden 2006); the fitting procedure in most of the recent work is the most attractive current regression strategy for low-frequency behaviour. Recently, Kravtsov *et al*. (2005) have developed a simplified nonlinear regression strategy that produces very good results for a three-layer quasi-geostrophic model with a realistic climate. However, order 2000 regression coefficients need to be fitted in a model with order 1000 state variables to achieve these results. Some inherent limitations of this approach in describing the correct physics are discussed briefly below in a simplified model.

All the work presented above derives reduced models by regression fitting of the resolved modes. Another approach is to take advantage of the basis function property of teleconnection patterns. Schubert (1985), Selten (1995), Achatz & Branstator (1999) and Achatz & Opsteegh (2003*a*,*b*) developed low-order models with empirical orthogonal functions (EOFs) as basis functions. Truncated EOF models experience climate drift due to the neglected interactions with the unresolved modes. Selten (1995) and Achatz & Branstator (1999) parametrize these neglected interactions by a linear damping, whose strength is determined empirically. A possibly more powerful tool to represent the dynamics of a system is principal interaction patterns (PIPs; Hasselmann 1988; Kwasniok 1996, 2004). The calculation of PIPs takes into account the dynamics of the model for which one tries to find an optimal basis and also often involves ad hoc closure through linear damping and an ansatz for nonlinear interactions. Crommelin & Majda (2004) compare different optimal bases. They find that the models based on PIPs are superior to models based on EOFs. On the other hand, they also point out that the determination of PIPs can show sensitivities regarding the calculation procedure, at least for some low-order atmospheric dynamical systems with regime transitions. Furthermore, PIPs have two more disadvantages: (i) much higher computational cost than EOFs, and (ii) one needs not only data as for EOFs but also the dynamical equations to calculate PIPs. These features can make PIPs possibly a less attractive basis.

Majda *et al*. (1999, 2001, 2002, 2003, 2005, 2006*a*,*b*) provide a systematic framework for how to account for the effect of the fast d.f. on the slow modes in combination with using the dominant teleconnection patterns as basis functions. In contrast to the empirical fitting procedures applied in the studies discussed above, the stochastic mode reduction strategy put forward in MTV *predicts* the functional form of all deterministic and stochastic correction terms and provides a *minimal* regression fitting procedure of only the *fast modes* (Franzke *et al*. 2005; Franzke & Majda 2006). In general, only an estimate for the variance and eddy turnover time for each fast mode is needed. It has been applied and tested on a wide variety of simplified models and examples.

### (a) Overview of the MTV strategy

We illustrate the ideas for stochastic climate modelling by considering the following prototype equation for geophysical flow:(2.1)The above functional form (2.1) is typical of dry dynamical cores of climate models, but the MTV strategy is easily extended to include non-quadratic nonlinearities such as associated with boundary fluxes and moist processes. In stochastic climate modelling, the variable *u* is decomposed into an orthogonal decomposition through the variables and *u*′, which are characterized by strongly differing time scales (MTV 1999; Majda *et al*. 2005). The variable denotes a slow low-frequency mode (also referred to as climate mode) of the system, which evolves slowly in time compared with the *u*′ variables (also referred to as fast mode). By decomposing in terms of some optimal energy norm basis, we can write them as(2.2)with and , where *R* is the number of climate modes; *a*_{i} denotes the expansion coefficients; and *α*_{i} and *β*_{j} are the expansion coefficients of the slow (fast) modes. The use of the energy norm ensures the conservation of energy by the nonlinear operator (Selten 1995). By properly projecting the energy norm basis, derived from the geophysical model, onto equation (2.1), we get two sets of equations for slow *α*_{i} and fast *β*_{i} modes(2.3)(2.4)where the nonlinear operators have been symmetrized, that is, *B*_{ijk}=*B*_{ikj} in (2.3) and (2.4). The upper indices *α* and *β* indicate the respective subsets of the full operators in (2.1). Here *ϵ* is a small positive parameter that controls the separation of time scale between the slow and fast modes and measures the ratio of the correlation time of the slowest non-climate mode *u*′ to the fastest climate mode . In placing the parameter in front of particular terms, we tacitly assume that they evolve on a faster time scale than the terms involving the climate modes alone (see MTV (2001) and Franzke *et al*. (2005) for more details). Ultimately, *ϵ* is set to the value *ϵ*=1 in developing all the final results (MTV 2002, 2003), that is, introducing *ϵ* is only a technical step in order to carry out the MTV mode reduction strategy. Such a use of *ϵ* has been checked on a wide variety of idealized examples where the actual value of *ϵ* ranges from quite small to order one (MTV 2002, 2003, 2006*b*; Majda & Timofeyev 2004). Following MTV (1999, 2001, 2002, 2003, 2006*b*) and Franzke *et al*. (2005), the mode elimination procedure is based on the assumption that the dynamics of the fast modes alone in equation (2.4), that is, the dynamical system(2.5)is ergodic and mixing with integrable decay of correlation. In other words, we assume that for almost all initial conditions, and suitable functions *f* and *g*, we have(2.6)where 〈.〉 denotes expectation with respect to some appropriate invariant distribution and(2.7)is an integrable function of *s*, that is, . Furthermore, we assume that the low-order statistics for the fast modes in (2.5) are Gaussian. Under the above assumptions, it can be shown in the limit *ϵ*→0 (Kurtz 1973; MTV 2001) that the dynamics of the slow modes *α*_{i} in (2.3) can be written as the following Itô stochastic equation for the slow modes alone:(2.8)where the nonlinear noise matrix *σ*^{(1)} satisfies(2.9)It is guaranteed (MTV 2001) that the operator on the l.h.s. of (2.9) is always positive definite ensuring the existence of the nonlinear noise matrix on the r.h.s. All coefficients are defined explicitly in MTV (2001) and Franzke *et al*. (2005). A comprehensive mathematical theory of the stochastic mode reduction strategy for geophysical applications is developed in MTV (2001) with many new mathematical phenomena in the resulting equations explored there.

To see which of these correction terms play a vital role in the integrations of the low-order stochastic model (2.8), we grouped the interaction terms between the slow and fast modes according to their physical origin and set a parameter *λ*_{i} in front of the corresponding interaction coefficient (see Franzke *et al*. 2005; Franzke & Majda 2006 for more details). The bare truncation is indicated by a *λ*_{B} and describes the interaction between the slow modes. The interaction between the triads *B*^{αββ} and *B*^{βαβ} gives rise to additive noise and a linear correction term and arises from the advection of the fast modes by the slow ones; we name these triads ‘additive’ triads and set a *λ*_{A} in front of them (MTV 1999, 2001, 2002, 2003). The other type of triad interaction is between *B*^{ααβ} and *B*^{βαα}. These interactions create multiplicative noises and cubic nonlinear correction terms (MTV 1999, 2001, 2002, 2003); we call them ‘multiplicative’ triads hereafter and indicate them by a *λ*_{M}. These triad interactions describe the advection of the slow modes by the fast modes that induce tendencies in the slow modes. The linear coupling between the slow and fast modes *L*^{αβ} and *L*^{βα} gives rise to additive noise and a linear correction term (MTV 2001; Franzke *et al*. 2005), which is called the augmented linearity here and is indicated by a *λ*_{L}. The augmented linearity describes the effect of the linear interaction between the fast (slow) modes and the climatological mean state onto the slow (fast) modes and is the main interaction captured in the linear stochastic modelling strategy (Delsole 2004). We set a *λ*_{F} in front of the last remaining interaction term *L*^{ββ}, the linear coupling of the fast modes. The quadratic nonlinear corrections, a forcing term and a further multiplicative noise contribution are caused by the interaction between the linear coupling terms and the multiplicative triads. Another forcing correction term comes from the interaction between additive triads and the linear coupling of the fast modes.

In ch. 3 of Majda *et al*. (2005), a simplified three-mode elementary ‘toy climate model’ is discussed and the MTV procedure is applied explicitly to that example. The origin of all the terms in (2.8) is developed in a transparent fashion in these examples. Once the low-order stochastic model has been developed from the above procedure, one can assess the importance of the various deterministic and stochastic processes systematically by varying the coefficients *λ*_{B}, *λ*_{A}, *λ*_{L}, *λ*_{M} and *λ*_{F} systematically in (2.8) (Franzke *et al*. 2005) and even develop simple physically motivated regression fitting strategies (Franzke & Majda 2006). In interesting recent work, Sura & Sardeshmukh (2008) have used scalar linear stochastic models with multiplicative and additive noise to explain non-Gaussian sea surface temperature (SST) variability. If the reduced stochastic models in (2.10) from the MTV procedure are linearized at the climate mean state, they automatically produce vector systems of linear stochastic equations with both multiplicative and additive noise with the same structure, with clear sources for the underlying physical contributions to this equation.

### (b) Idealized models for stochastic mode reduction

The idealized models, where the procedure has been tested, have order 100 d.f. and include those with trivial climates (MTV 2002), periodic orbits or multiple equilibria (MTV 2003) and heteroclinic chaotic orbits coupled to a deterministic bath of modes satisfying the truncated Burgers–Hopf equation (Majda & Timofeyev 2000, 2004). The truncated Burgers–Hopf equation is a toy model with some remarkable features mimicking behaviour in the real atmosphere; it has a well-defined equipartition spectrum and a simple scaling theory for correlations with the large scales decorrelating more slowly than the small scales, that is, low-frequency variability. Furthermore, these predictions are confirmed with very high precision by numerical simulations (Majda & Timofeyev 2000; Majda & Wang 2006). The MTV procedure has been validated in these examples even when there is little separation of time scales between the slow and fast modes. In the example of a four-dimensional resonant system with chaotic dynamics coupled to the truncated Burgers system, only one empirical regression fitting coefficient is used and complex bifurcation diagrams and probability distribution functions (pdfs) in a climate change scenario are reproduced by the four-dimensional stochastic mode reduction resulting from the MTV procedure applied to the 104 d.f. deterministic system. An especially stringent recent test is the application of this procedure to the first few large-scale modes of the truncated Burgers equations in the turbulent cascade (MTV 2006*b*).

### (c) A priori stochastic modelling for mountain torque

The ideal barotropic quasi-geostrophic equations with a large-scale zonal mean flow *U* on a 2*π*×2*π* periodic domain (Carnevale & Frederiksen 1987) are given by(2.10)where *q* is the potential vorticity; *U* is the large-scale zonal mean flow; *ψ* is the stream function; and *h* is the topography. In (2.10), the mean flow changes in time through the topographic stress; this effect is the direct analogue for periodic geometry of the change in time of angular momentum due to mountain torque in spherical geometry (Frederiksen *et al*. 1996; Majda & Wang 2006). Here the *a priori* stochastic modelling strategy (MTV 1999, 2001) is applied to the stochastic modelling of the topographic stress terms in (2.10) as an analogue for mountain torque; thus, the variable *U* is the slow variable while all the modes *ψ*_{k} are fast variables for the MTV procedure.

In this example, the systematic stochastic modelling procedure (MTV 2003) results in the predicted nonlinear reduced equation for *U*,(2.11)where and(2.12)Here is the Rossby wave frequency Doppler shifted by the mean flow. Under the additional assumption that , a standard predicted linear stochastic model for *U* emerges from (2.11) with *γ*=*γ*(0) from (2.12) and *γ*′(0)=0 (MTV 1999, 2001; see Majda *et al*. (2003) for definitions of *μ* and *γ*_{k}). It is shown in MTV (2003) that this nonlinear stochastic equation is superior to the linear one for large amplitude topography, where *ϵ*≅0.7. This is the simplest example with multiplicative noise. Egger (2005) has used the systematic strategies from MTV (1999, 2003) to improve regression strategies for analysing observational data for angular momentum. In Majda *et al*. (2006*a*), the recent systematic regression fitting strategy mentioned earlier (Crommelin & Vanden-Eijnden 2006) is applied to (2.10) and independently confirms the predictions in (2.11) and (2.12).

### (d) Geophysical and climate models

Franzke *et al*. (2005) put the above systematic stochastic mode reduction strategy in a form that makes the practical implementation of the MTV procedure in the complex geophysical models simpler with the same reduced stochastic equations for the fast modes. EOFs in an appropriate metric are used to distinguish between the slow and fast modes in high-dimensional geophysical systems, such as climate models. Even though EOFs do not strictly decompose the modes according to their autocorrelation time scale, the leading EOFs usually evolve on a slower time scale than the higher EOFs. In the study by Franzke *et al*. (2005), a T21-truncated barotropic model on the sphere with a realistic climate was used to derive low-order stochastic models based on kinetic energy norm EOFs by the MTV strategy. Low-order models with as little as two slow modes succeed in capturing the geographical distributions of the climatological mean field, the variance and the eddy forcing. Furthermore, the envelope of the autocorrelation functions is captured reasonably well.

Recently, Franzke & Majda (2006) applied the systematic stochastic mode reduction strategy to a baroclinic three-layer quasi-geostrophic model on the sphere (Marshall & Molteni 1993), which mimics the climatology of the European Centre for Medium-Range Weather Forecasts reanalysis data. The low-order stochastic climate model consists of climate modes as slow modes defined as the leading total energy norm EOFs and the stochastic mode reduction procedure predicts all forcing, linear, quadratic and cubic correction terms as well as additive and multiplicative noises; these correction terms and noises account for the interaction of the climate modes with the neglected non-climate modes and the self-interaction among the non-climate modes. For the three-layer quasi-geostrophic model low-order stochastic models with 10 or less climate modes reproduce the geographical distributions of the standard deviation and eddy forcing well. They underestimate the standard deviations by at most a factor of approximately 1.5. Furthermore, they reproduce the autocorrelation functions reasonably well. A budget analysis shows that both linear and nonlinear correction terms as well as both additive and multiplicative noises are important. The physical intuition behind the noises as derived from the MTV procedure is as follows: the additive noise stems from the linear interaction between the fast modes and the climatological mean state and the multiplicative noise comes from the advection of the slow modes by the fast modes. All these deterministic correction terms and noises (both additive and multiplicative) are *predicted* by the systematic stochastic mode reduction strategy, whereas previous studies *a priori* approximate the nonlinear part of the equations by a linear operator and additive noise. This noise is typically white in time but may be spatially correlated. In other words, these studies truncate the dynamics on both the slow and fast modes and add ad hoc damping in order to stabilize the linear model (Whitaker & Sardeshmukh 1998; Zhang & Held 1999). The systematic MTV approach summarized briefly above truncates the dynamics only on the fast modes and predicts the functional form of all necessary nonlinear correction terms and noises; therefore, it also predicts the necessary damping.

The MTV stochastic climate models for this application experience some climate drift. A minimal empirical MTV model without climate drift can be constructed through three-parameter regression fitting by downscaling the bare truncation terms and upscaling the two important MTV processes (augmented linearity and multiplicative triads). These empirical MTV stochastic climate models with minimal regression fitting still capture the geographical distribution of the standard deviation and eddy forcing and the autocorrelation functions reasonably well, while not experiencing climate drift. This surprising result can be interpreted as the fact that the climate modes are predominantly driven by the fast modes and the self-interactions among the slow modes are less important, as can already be seen from the bare truncation models, which do not capture any feature of the actual dynamics. Furthermore, these empirical MTV stochastic climate models suggest that the bare truncation is probably the cause of the climate drift. Integrations of bare truncation models (without any MTV correction terms) already produce a big climate drift (Franzke & Majda 2006). The MTV mode reduction procedure is able to reduce the climate drift in most of the slow modes, but is not able to overcome it completely. Previous results with a variety of simplified models show no climate drift in an MTV framework (MTV 1999, 2001, 2002, 2003; Majda & Timofeyev 2004). This is probably because these simplified models are constructed in such a way that they have an optimal basis that captures the dynamics of the climate modes. This gives evidence that total energy norm EOFs are not an adequate dynamical basis in capturing the dynamics of the slow modes. Further details of this application can be found in Franzke & Majda (2006).

### (e) A simple stochastic model with key features of atmospheric low-frequency variability

In this section, we present a four-mode stochastic climate model of the kind put forward in Majda *et al*. (2005, ch. 3). This simple stochastic climate model is set-up in such a way that it features many of the important dynamical features of comprehensive GCMs but with many fewer d.f. Such simple toy models allow the efficient exploration of the whole parameter space that is impossible to conduct with GCMs. Thus, we are able to test the predictions of the above framework with direct model experiments by switching on and off certain terms rather than relying on diagnostic methods.

While this model is not rigorously derived from a geophysical flow model (e.g. barotropic vorticity equation), it has the same functional form one would end up with when deriving a reduced stochastic model from a geophysical model. Thus, consistent with geophysical flow models, the toy model has a quadratically nonlinear part that conserves energy, a linear operator and a constant forcing, which in a geophysical model would represent the effects of external forcing such as solar insulation and sea surface temperature. The linear operator has two contributions: one is a skew-symmetric part formally similar to the Coriolis effect and topographic Rossby wave propagation. The other is a negative definite symmetric part formally similar to dissipative processes such as surface drag and radiative damping.

The model is constructed in such a way that there are two modes, denoted by ** x**, that evolve more slowly than the other two modes,

**. In the realistic models, there would be very many additional fast modes representing, for example, synoptic weather systems or convection. To mimic their combined effect, we include damping and stochastic forcing in the equations for**

*y***, where**

*y***denotes a Wiener process. The motivation for this approximation is that these fast modes are associated with turbulent energy transfers and strong mixing and that we do not require a more detailed description since we are only interested in their effect on the slow modes. The two fast modes carry most of the variance in this model but as noted earlier, these two modes are surrogates in the model for the entire bath of fast modes so this is very natural. Note that this stochastic Ansatz is different from most of the previous studies (e.g. Newman**

*W**et al*. 1997; Whitaker & Sardeshmukh 1998) that truncated the model on

*both*the resolved and unresolved modes and replaced all nonlinear terms by a stochastic process. In this study, only the nonlinear interactions among the fast modes are modelled by this stochastic Ansatz. The parameter

*ϵ*controls the time-scale separation between the slow and fast variables. For testing the predictions of the general framework that is derived in the previous section, we will treat the two slow modes

**as the climate modes and the two fast modes as the non-climate modes**

*x***. Therefore, our toy model has the following form:(2.13a)(2.13b)(2.13c)(2.13d)Note that in this case we use a different scaling in powers of**

*y**ϵ*in order to derive easily a solution by working directly on the stochastic differential equation instead of the Fokker–Planck equation as for equation (2.8) (MTV 1999, 2001). To ensure energy conservation of the nonlinear operator, the coefficients have to satisfy the following relation:

*b*

_{123}+

*b*

_{213}+

*b*

_{312}=0, while the nonlinear bare truncation terms also conserve energy. In this particular set-up, the slow and the fast modes are coupled through two mechanisms: one is a skew-symmetric linear coupling and the other is nonlinear triad interaction. The nonlinear coupling involving

*b*

_{ijk}produces multiplicative noise in the MTV framework (1999, 2003; see also equation (2.15

*a*)–(2.15

*d*)) so we refer to it as a multiplicative triad. One advantage of the model in (2.13

*a*)–(2.13

*d*) is that the stochastic mode reduction as

*ϵ*→0 is done explicitly through elementary manipulations (MTV 1999, 2001; Majda

*et al*. 2005). The corresponding reduced Itô stochastic differential equation (SDE) for the climate variables alone is given by(2.14a)(2.14b)Note that coarse-graining time as

*t*→

*t*/

*ϵ*amounts to setting

*ϵ*=1, because by coarse-graining time the factors of

*ϵ*disappear (see Majda

*et al*. (2005) for more details).

To evaluate the performance of the reduced dynamics, we calculate the autocorrelation function and the third-order two-time moment , which is a measure of deviations from Gaussianity (MTV 2002; Majda *et al*. 2005). The comparison of the reduced model (2.14*a*) and (2.14*b*) with the full model (2.13*a*)–(2.13*d*) results shows excellent agreement for moderate values of *ϵ*=0.1, 0.5 and still good agreement for *ϵ*=1.0 (figures 1 and 2). Especially, the non-Gaussian features are reproduced with high accuracy for *ϵ*=0.1 and also *ϵ*=0.5, even though slightly less well. For large values of *ϵ*, the reduced dynamics get the sign of the non-Gaussianity right. All of the reduced models in (2.11)–(2.14*a*) and (2.14*b*) cannot be approximated by the interesting regression strategy of Kravtsov *et al*. (2005) because there is nonlinear multiplicative noise in (2.13*a*)–(2.13*d*), (2.14*a*) and (2.14*b*), nonlinear triad interaction in (2.13*a*)–(2.13*d*) and augmented cubic nonlinearity in (2.14*a*) and (2.14*b*). Thus, these regression strategies necessarily have large model error in this example.

### (f) A mathematical framework for the mean tendency equation as a dynamical diagnostic

In this section, we provide a general framework to estimate the origin of nonlinear signatures of planetary wave dynamics and how important the observed deviations from Gaussianity are for the planetary waves (Franzke *et al*. 2007). This general framework diagnoses contributions to the mean values of state-dependent tendencies in a low-dimensional subspace of complex geophysical systems. The mean phase-space tendencies in some GCMs show distinct nonlinear signatures in certain planes, which are spanned by its leading EOFs (Selten & Branstator 2004; Branstator & Berner 2005; Franzke *et al*. 2007) that also show only weak deviations from Gaussianity; mostly in the form of weak skewness and kurtosis and in the case of joint pdfs multiple radial ridges of enhanced density (Franzke & Majda 2006; Berner & Branstator 2007). This mathematical framework can be applied to complex geophysical systems and reveals how much the self-interaction among the modes spanning those planes (climate modes) and how much the unresolved modes contribute to these mean phase-space tendencies and also the effect of the observed small deviations from Gaussianity.

To derive the mean tendency equation for the resolved modes, we split the state vector of a quadratic dynamical system into resolved modes *α* and unresolved modes *β* and use conditional mean probability density relations. Note that while there is a similarity in motivation in decomposing the system into the resolved and unresolved modes compared with the MTV strategy, where the system is decomposed into the slow and fast modes, the mean tendency equation is more general such that it applies also to systems without time-scale separation. The general formula for the *conditional mean tendencies* for the dynamics in the resolved variables is derived in Franzke *et al*. (2007) and is given by(2.15a)(2.15b)(2.15c)(2.15d)Note that the r.h.s. of (2.15*a*) is the *bare truncation* restricted to interactions among the resolved modes, while (2.15*b*) and (2.15*c*) involve all of the conditional mean statistics ; the terms of (2.15*c*) are associated with *multiplicative triad interactions* (leading to multiplicative noise in an MTV (1999, 2003) framework of *β*_{j}; the terms in (2.15*d*) involve all of the conditional interaction statistics of second moments and include all of the *additive triad interactions* (leading to additive noise in an MTV framework) of *β*_{j} and *β*_{k}. Thus, equation (2.15*a*)–(2.15*d*) provides a general framework to investigate the mean phase-space tendencies, which allows the decomposition into contributions from interactions among the resolved planetary waves themselves and various contributions involving unresolved d.f.

Now we simplify the above conditional mean tendency equation for purely Gaussian EOF modes. The pdfs of EOFs of GCM data are nearly Gaussian (Hsu & Zwiers 2001; Franzke *et al*. 2005; Franzke & Majda 2006; Majda *et al*. 2006*a*; Berner & Branstator 2007) and there are many geophysical models without damping and forcing that exactly satisfy these assumptions such as barotropic flow on the sphere with topography (Carnevale & Frederiksen 1987; Salmon 1998; Majda & Wang 2006); thus, it is reasonable as a starting point to *assume* that the pdf is *exactly Gaussian*. Thus, in the EOF basis, the pdf factors like(2.16)where and are Gaussian distributions with mean zero. Thus, in the *Gaussian* case, the conditional mean tendency equation simplifies to (see Franzke *et al*. (2007) for more details)(2.17a)(2.17b)Note that contributions from (2.15*b*) and (2.15*c*), that is, linear coupling and multiplicative triad interactions, are identically zero. These terms vanish because the modes are uncorrelated in the EOF basis. Thus, in this Gaussian case, the conditional mean tendency equation recognizes bare truncation (2.15*a*) and a constant forcing from additive triad interactions (2.15*d*). Since the leading EOFs of GCM data have pdfs with only small, but significant departures from Gaussianity the behaviour in (2.16), (2.17*a*) and (2.17*b*) serves as a ‘null hypothesis’ for these deviations from Gaussianity in the climate models.

This general framework predicts that in the case of purely Gaussian modes the nonlinear signatures are stemming from the bare truncation (i.e. the self-interaction among the planetary waves resolved in the low-dimensional plane). In Franzke *et al*. (2007), these diagnostics were applied to a plane of two low-frequency EOFs in a three-layer climate model with a nonlinear double swirl (see figure 5a in Franzke *et al*. 2007); the origin of this double swirl is primarily the three contributions in (2.15*b*)–(2.15*d*) from the unresolved modes and not nonlinear effects from the bare truncation. The EOFs of the three-layer model have small but significant departures from Gaussianity (see figure 9 in Franzke *et al*. 2007). Therefore, these deviations account for the contribution of the unresolved onto the resolved modes in the mean phase-space tendencies. Thus, these results show that small deviations from Gaussianity have a substantial impact on the dynamics of the leading EOFs.

## 3. Coarse-grained stochastic lattice models for climate: tropical convection

The current practical models for prediction of both weather and climate involve GCMs where the physical equations for these extremely complex flows are discretized in space and time, and the effects of unresolved processes are parametrized according to various recipes. With the current generation of supercomputers, the smallest possible mesh spacings are approximately 10–50 km for short-term weather simulations and of order 100 km for short-term climate simulations. There are many important physical processes that are unresolved in such simulations such as the mesoscale sea ice cover, the cloud cover in subtropical boundary layers and deep convective clouds in the tropics. Most of these features are highly intermittent in space and time. An appealing way to represent these unresolved features is through a suitable coarse-grained stochastic model that simultaneously retains crucial physical features of the interaction between the unresolved and resolved scales in a GCM. In work from 2002 and 2003, two of the authors have developed a new systematic stochastic strategy (Majda & Khouider 2002; Khouider *et al*. 2003) to parametrize key features of deep convection in the tropics involving suitable stochastic spin-flip models and also a systematic mathematical strategy to coarse grain such microscopic stochastic models to practical mesoscopic meshes in a computationally efficient manner while retaining crucial physical properties of the interaction.

As regards tropical convection, crucial scientific issues involve the fashion in which a stochastic model effects the climate mean state and the strength and nature of fluctuations about the climate mean. Here the strategy to develop a new family of coarse-grained stochastic models for tropical deep convection is briefly reviewed (Majda & Khouider 2002; Khouider *et al*. 2003) as an illustrative example of the potential use of stochastic lattice models. In Khouider *et al*. (2003), it has been established that in suitable regimes of parameters, the coarse-grained stochastic parametrizations can significantly alter the climatology as well as increase wave fluctuations about the climatology. This was established in Khouider *et al*. (2003) in the simplest scenario for tropical climate involving the Walker circulation, the east–west climatological state that arises from local region of enhanced surface heat flux, mimicking the Indonesian marine continent. Convectively coupled waves in the tropics such as the Madden–Julian oscillation play an important role in medium-range forecasts, yet the current generation of computer models fails to represent such waves adequately (Lin *et al*. 2006). Palmer (2001) has emphasized the potential of stochastic parametrization to reduce the model error in a deterministic computer model. Here in an idealized setting, we show how to develop a stochastic parametrization to modify and improve the behaviour of convectively coupled waves in a reasonable prototype GCM; this is achieved by following a path guided by the systematic design principle for the idealized model in (1.1) (KMS 2006, 2007) to build in suitable intermittency effects.

### (a) The microscopic stochastic model for convective inhibition

In a typical GCM, the fluid dynamical and thermodynamical variables denoted here by the generic vector ** u** are regarded as known only over a discrete horizontal mesh with denoting these discrete values. Throughout the discussion, one horizontal spatial dimension along the equator in the east–west direction is assumed for simplicity in notation and explanation. As mentioned above, the typical mesh spacing in a GCM is coarse with Δ

*x*ranging from 50 to 250 km depending on the time duration of the simulation. The stochastic variable used to illustrate the approach is convective inhibition (CIN). Observationally, CIN is known to have significant fluctuations on a horizontal spatial scale on the order of a kilometre, the microscopic scale here, with changes in CIN attributed to different mechanisms in the turbulent boundary layer, such as gust fronts, gravity waves and turbulent fluctuations in equivalent potential temperature (Mapes 2000). In Majda & Khouider (2002) and Khouider

*et al*. (2003), it was proposed that all of these different microscopic physical mechanisms and their interaction that increase and decrease CIN are too complex to model in detail in a coarse-mesh GCM parametrization and instead, as in statistical mechanics, should be modelled by a simple order parameter

*σ*

_{I}, taking only two discrete values,(3.1)

The value of CIN at a given coarse-mesh point is determined by the averaging of CIN over the microscopic states in the vicinity of the given mesh point, i.e.(3.2)

Note that the mesh size Δ*x* is mesoscopic, that is, between the microscale *O*(1 km) and the macroscale *O*(10 000 km), and that can have any value in the range . Discrete sums over microscopic mesh values (of order 1 km) and continuous integrals are used interchangeably for notational convenience (Majda & Khouider 2002).

### (b) The simplest coarse-grained stochastic model

In practical parametrization, it is desirable for computational feasibility to replace the microscopic dynamics by a process on the coarse mesh, which retains critical dynamical features of the interaction. Following the general procedure developed and tested in Katsoulakis & Vlachos (2003) and Katsoulakis *et al*. (2003*a*,*b*), the simplest local version of the systematic coarse-grained stochastic process is developed in Khouider *et al*. (2003) and summarized here.

Each coarse cell Δ*x*_{k}, *k*=1, …, *m*, of the coarse-grained lattice is divided into *q* microscopic cells, such that Δ*x*_{k}≡1/*q* {1, 2, …, *q*}, *k*=1, …, *m*. In the coarse-grained procedure, given the coarse-grained sequence of random variables(3.3)so that the average in (3.2) verifies , for *j*=*k* in some sense, the microscopic dynamics is replaced by a birth/death Markov process defined on the variables {0, 1, …, *q*} for each *k*, such that *η*_{t}(*k*) evolves according to the following probability law:(3.4)Here *C*_{a} and *C*_{d} are the coarse-grained adsorption and desorption rates representing the rates at which CIN sites are being created and destroyed, respectively, within one coarse-grained cell, that is, the birth/death rates for the Markov process is(3.5)where(3.6)with the coarse-grained interaction potential within the coarse cell given by , where *U*_{0} is the mean strength of the microscopic interaction potential *J* between neighbouring sites (Katsoulakis *et al*. 2003*a*,*b*). The coarse-grained energy content for CIN is given by the coarse-grained Hamiltonian(3.7)where *h*_{ext} is the external potential that modifies the energy for CIN depending on the large-scale flow variables. The canonical invariant Gibbs measure for the coarse-grained stochastic process is a product measure given by(3.8)where is a given explicit/prior distribution (Katsoulakis *et al*. 2003*b*). As shown in Katsoulakis *et al*. (2003*b*), the coarse-grained birth/death process above satisfies detailed balance with respect to the Gibbs measure in (3.8) as well as a number of other attractive theoretical features. The simplest coarse-grained approximation given above assumes that the effect of the microscopic interactions on the mesoscopic scales occurs within the mesoscopic coarse-mesh scale Δ*x*, otherwise systematic non-local couplings are needed (Katsoulakis *et al*. 2003*b*). The accuracy of these approximations is tested for diverse examples from material science elsewhere (Katsoulakis & Vlachos 2003; Katsoulakis *et al*. 2003*a*,*b*) and for the instructive idealized coupled models in (1.1) (KMS 2004, 2005, 2006, 2007).

The practical implementation of the coarse-grained birth/death process in (3.3)–(3.6) requires specification of the parameters, *τ*_{I}, *U*_{0}, *q* and the external potential as well as the statistical parameter *β*. The advantages of such a stochastic lattice model are as follows:

retains systematically the energetics of unresolved features through the coarse-grained Gibb's measure;

has minimal computational overhead since there are rapid algorithms for updating birth–death processes;

incorporates feedbacks of the resolved modes on the unresolved modes and their energetics through an external field; and

includes dynamical coupling through not only sampling the probability distributions of unresolved variables but also their evolving behaviour in time is constrained by the large-scale dynamics.

### (c) The model deterministic convective parametrization

A prototype mass flux parametrization with crude vertical resolution (Majda & Shefter 2001; Khouider & Majda 2006*b*) is used to illustrate the fashion in which the coarse-grained stochastic model for CIN can be coupled to a deterministic convective mass flux parametrization. The prognostic variables are the *x*-component of the fluid velocity, *u*, the potential temperature in the middle troposphere, *θ*, the equivalent potential temperatures, *θ*_{eb} and *θ*_{em}, measuring the potential temperatures plus moisture content of the boundary layer and middle troposphere, respectively. The vertical structure is determined by projection on a first baroclinic heating mode (Majda & Shefter 2001; Khouider & Majda 2006*b*). The dynamic equations for these variables in the parametrization are given by(3.9)while the constants and are externally imposed and represent the radiative cooling at equilibrium in the upper troposphere and saturation equivalent potential temperature in the boundary layer. The constants *h* and *H* measure the depths of the boundary layer and the troposphere above the boundary layer, respectively. The typical values used here are *h*=500 m and *H*=16 km, while *u*_{0}=2 m s^{−1} and , corresponding to a momentum spin up time of approximately 3 days. The explicit values for the other constants used in (3.9) and elsewhere in this section can be found in Majda & Shefter (2001) and Khouider *et al*. (2003).

The crucial quantities in the prototype mass flux parametrization are the terms and *D*, where represents the middle troposphere heating due to deep convection, while *D* represents the downward mass flux on the boundary layer. The heating term is given by(3.10)where *M* is a fixed constant; *σ*_{c} is the area fraction for deep convective mass flux; and is the convectively available potential energy. Here *R* is a dimensional constant (Majda & Shefter 2001; Khouider *et al*. 2003). The downward mass flux on the boundary layer *D* includes the environmental downdrafts *m*_{e} and the downward mass flux due to convection *m*_{−}, which are non-negative quantities with explicit formulae described elsewhere (Majda & Shefter 2001; Khouider *et al*. 2003). The notation (*X*)^{+} denotes the positive part of *X*. This parametrization respects conservation of vertically integrated moist static energy.

The equations (3.9) and (3.10) represent an idealized GCM with crude vertical resolution based on a reasonable design principle for deep convection including basic conservation principles. However, like the current generation of GCMs, the model has major deficiencies as regards convectively coupled waves. First, instabilities for the full model in (3.9) arise only through the nonlinearity from and are called the wind-induced surface heat exchange (WISHE) instability (Majda & Shefter 2001, and references therein). Second, if the bulk aerodynamic flux coefficients in (3.9), are replaced by the constant value |*u*_{0}|, that is no WISHE, all waves in the system are stable and no fluctuations emerge in an aquaplanet simulation (Majda & Shefter 2001). There is no observational evidence supporting the role of WISHE in driving the large-scale convectively coupled waves in nature (Lin *et al*. 2006 and references therein); here the WISHE term is regarded as a deterministic fix in a GCM parametrization to generate instabilities. The simulation results reported in figure 5*a* show that in an aquaplanet model above the equator, two regular periodic convectively coupled waves moving eastward and westward at roughly equal strength are generated by the idealized GCM. Since GCMs often have convectively coupled waves that move too fast and are too regular (Lin *et al*. 2006), the goal here is to see whether the stochastic lattice coupling will slow down the waves in the deterministic parametrization and simultaneously increase the spatio-temporal fluctuations in those waves. It is important to note here that there are recent deterministic multicloud models (Khouider & Majda 2006*a*, 2007, 2008) for convectively coupled waves involving the three cloud types in observations above the boundary layer, congestus, deep and stratiform, and their heating structure that reproduces key features of the observational record for convectively coupled waves (Lin *et al*. 2006). The mechanism of instability in these models (Khouider & Majda 2006*a*) is completely different from WISHE that is not active in the multicloud models. For the idealized setting of flow above the equator, the multicloud models can produce packets of convectively coupled waves moving in one direction at 15–20 m s^{−1} with their low-frequency envelopes moving at 4–7 m s^{−1} in the opposite direction across the warm pool in a fashion like the Madden–Julian oscillation (Khouider & Majda 2007, 2008).

### (d) Coupling of the stochastic CIN model into the parametrization

The equations (3.9) and (3.10) are regarded here as the prototype deterministic GCM parametrization when discretized in a standard fashion using central differences on a coarse mesh Δ*x* with Δ*x* ranging from 50 to 250 km. In the simulations from Khouider *et al*. (2003), and those presented below, Δ*x*=80 km. The coarse-grained stochastic CIN model is coupled to this basic parametrization. The area fraction for deep convection *σ*_{c} governing the upward mass flux strength is allowed to vary on the coarse mesh and is given by(3.11)where is the average in (3.2) and is a threshold constant and equal to 0.002 (Majda & Shefter 2001; Khouider *et al*. 2003). When the order parameter *σ*_{I} signifies strong CIN locally so that the flux of deep convection is diminished to zero while with PAC locally active , this flux increases to the maximum allowed by the value . To complete the coupling of the stochastic CIN model into the parametrization, the coarse-mesh external potential , from (3.6) and (3.7), needs to be specified from the coarse-mesh values *u*_{j}. There is no unique choice of the external potential but its form can be dictated by simple physical reasoning. In Khouider *et al*. (2003), the plausible physical assumption is made that when the convective downward mass flux *m*_{−} decreases the energy for CIN decreases. Since the convective downward mass flux results from the evaporative cooling induced by precipitation falling into dry air, it constitutes a mechanism that carries negatively buoyant cool and dry air from the middle troposphere onto the boundary layer, hence tending to reduce CAPE and deep convection.

Another natural external potential is the boundary-layer equivalent potential, since the flux at the boundary is crucial physically. As such, it yields somewhat stronger and more intermittent fluctuations. Thus, the choice(3.12)is used here with a calibration factor. The other parameters in the stochastic lattice model are chosen as *τ*_{I}=2 hours, *β*=1, *U*_{0}=1 so that CIN sites are favoured in the equilibrium Gibbs measure.

### (e) The stochastic single-column model and intermittency

A central issue is how to calibrate the stochastic lattice model to generate intermittent fluctuations with plausible magnitudes as observed in tropical convection. A natural design framework is first to achieve such behaviour in the stochastic single-column model given by the following equations:(3.13)Here *σ*_{I} is computed from the birth–death stochastic model in (3.3)–(3.6) with the interaction potential defined through *θ*_{eb} from (3.7) and the definition of *h*_{ext} in (3.12) above. The different time scales were computed for a radiative convective equilibrium consistent with the Jordan tropical sounding and the values *τ*_{eb}≈6 hours, *τ*_{em}≈8 days, *τ*_{e}≈8 hours, K d^{−1}, *τ*_{R}=50 days and *τ*_{D}=75 days are used here (Khouider *et al*. 2003). This is a specific example of the prototype models from (1.1) studied extensively in idealized settings (KMS 2004, 2005, 2006). As mentioned earlier, the choices *β*=1, *U*_{0}=1 guarantee that the Gibbs measure in (3.8) has a high probability for CIN states; this is a natural requirement for deep convection, where the area fraction of deep convection is much smaller than the area with positive CAPE. In order to generate intermittency and fluctuations in the stochastic column model, it is very natural to require that the ratio of the creation rate of PAC sites to the creation rate of CIN sites, *C*_{d}/*C*_{a}, becomes large as *σ*_{I}↗1 (KMS 2006, 2007). In figure 3, it is shown that this property of the stochastic model is satisfied for the coefficient in the external potential in (3.12) with similar behaviour for *q*=10, 100 and 1000. Recall that *q* is the number of microscopic states in the birth–death process. A time series of the stochastic column model implemented with these parameters is displayed in figure 4 with *q*=10. The strongly intermittent fluctuations in *θ*_{eb} over several degrees Kelvin as well as similar intermittency in the mass flux is evident. They are consistent with SST variations observed in the Indian Ocean warm pool (e.g. Fasullo & Webster 1999). Also note that the fluctuations in the mid-troposphere potential temperature are much weaker in magnitude as occurs in the actual tropics, where middle tropospheric temperature anomalies are on the order of 0.1–0.5°C (e.g. Wheeler *et al*. 2000).

### (f) The effects of the stochastic parametrization on convectively coupled waves

All parameters in the stochastic lattice model have been determined through systematic design principles in the previous section. Here the results of numerical simulations of the stochastic model in statistical steady state are reported for flow above the equator in a standard aquaplanet set-up with uniform SST (Majda & Khouider 2002; Khouider *et al*. 2003). The results of various simulations are reported in figure 5 for comparison with the deterministic parametrization with WISHE reported in figure 5*a*. As shown in figure 5*b*, the effect of the fluctuations of the stochastic lattice model is simultaneously to create more realistic intermittency in the convectively coupled waves and to slow down their phase speed from 15 to 11 m s^{−1}; as mentioned earlier, these are desirable qualitative features of a stochastic parametrization. In figure 5*c*, we report the result of running the model parametrization in (3.9) without WISHE but coupled to the stochastic lattice model; recall that this is the situation where the deterministic model without WISHE produces no waves in the statistical steady state. In figure 5*c*, there is a clear evidence for stochastic-generated convectively coupled waves with the reduced phase speed of roughly 8 m s^{−1} and roughly half the amplitude; these waves have been created by coupling alone to the intermittent stochastic lattice model and without deterministic instability. This is another attractive feature of the present stochastic lattice models in changing the character of model error for convection (Palmer 2001).

## 4. Concluding remarks

This paper both reviewed and provided new illustrations and examples of the fashion in which modern applied mathematics can provide new perspectives and systematic design principles for stochastic modelling for climate. Section 2 was devoted to systematic stochastic modelling of low-frequency variability including the quantitative sources for multiplicative noise. A new simplified low-dimensional stochastic model with key features of atmospheric low-frequency variability was introduced in §2*e* in order to test stochastic mode reduction strategies. A recent diagnostic test with firm mathematical underpinning for exploring the subtle departures from Gaussianity and their sources was discussed in §2*f*. Section 3 was devoted to developing stochastic lattice models to capture intermittent features and improve the fidelity of deterministic parametrizations. Recent systematic design principles (KMS 2006, 2007) were used to calibrate a stochastic column model for tropical convection in §3*e*; the practical effect of slowing down convectively coupled waves and increasing their fluctuations through these stochastic lattice models was presented in §3*f*.

## Acknowledgments

The authors thank their collaborators, Eric Vanden-Eijnden, Ilya Timofeyev, Markos Katsoulakis and Alex Sopasakis, for their explicit and implicit contributions to the work presented here. The research of A.J.M. is partially funded by the National Science Foundation grant no. DMS-0456713, the Office of Naval Research grant no. N00014-05-1-0164 and the Defense Advanced Research Project Agency grant no. N00014-07-1-0750. The research of B.K. is partly supported by a grant from the National Sciences and Engineering Research Council of Canada.

## Footnotes

One contribution of 12 to a Theme Issue ‘Stochastic physics and climate modelling’.

- © 2008 The Royal Society