## Abstract

Mean phase space tendencies are investigated to systematically identify the origin of nonlinear signatures and the dynamical significance of small deviations from Gaussianity of planetary low-frequency waves. A general framework for the systematic investigation of mean phase space tendencies in complex geophysical systems is derived. In the special case of purely Gaussian statistics, this theory predicts that the interactions among the planetary waves themselves are the source of the nonlinear signatures in phase space, whereas the unresolved waves contribute only an amplitude-independent forcing, and cannot contribute to any nonlinear signature.

The predictions of the general framework are studied for a simple stochastic climate model. This toy model has statistics that are very close to being Gaussian and a strong nonlinear signature in the form of a double swirl in the mean phase space tendencies of its low-frequency variables, much like recently identified signatures of nonlinear planetary wave dynamics in prototype and comprehensive atmospheric general circulation models (GCMs). As predicted by the general framework for the Gaussian case, the double swirl results from nonlinear interactions of the low-frequency variables. Mean phase space tendencies in a reduced space of a prototype atmospheric GCM are also investigated. Analysis of the dynamics producing nonlinear signatures in these mean tendencies shows a complex interplay between waves resolved in the subspace and unresolved waves. The interactions among the resolved planetary waves themselves do not produce the nonlinear signature. It is the interaction with the unresolved waves that is responsible for the nonlinear dynamics. Comparing this result with the predictions of the general framework for the Gaussian case shows that the impact of the unresolved waves is due to their small deviations from Gaussianity. This suggests that the observed deviations from Gaussianity, even though small, are dynamically relevant.

## 1. Introduction

The driving forces and dominant functional form of planetary wave dynamics is an important topic in recent literature. Depending on the school of thought, planetary wave dynamics are either dominated by deterministic linear dynamics that are driven by additive white noise (Penland and Sardeshmukh 1995; Newman et al. 1997; Branstator and Haupt 1998; Branstator and Frederiksen 2003) or are intrinsically nonlinear. The latter view is motivated by studies of highly truncated models (Charney and DeVore 1979; Wiin-Nielsen 1979; Charney and Strauss 1980; Schubert 1985; Selten 1995; Ghil and Robertson 2002) and by the behavior of low-order dynamical systems that mimic certain aspects of the climate system in an idealized fashion (Palmer 1999).

Since these low-frequency planetary waves have a significant impact on long-range predictability, regional climate, and climate change, a better understanding of planetary wave dynamics is not only of academic interest but will be crucial in improving long-range weather and climate predictions. The better understanding will also guide the development of systematically reduced stochastic low-order models of low-frequency waves for predictability and climate change response studies and for efficient modeling of the coupled atmosphere–ocean system.

The climate system is governed by fundamentally nonlinear equations of motion. However, it is not obvious if the evolution of the dominant planetary waves is strongly nonlinear, because some of their key characteristics have been explained largely in linear terms (Rossby 1939; Hoskins and Karoly 1981; Simmons et al. 1983; Branstator 1990). Furthermore, there are also nonlinear processes that can be represented in a linear fashion. One example for such a nonlinear process, which can be explained by a linear relationship with planetary waves, is momentum fluxes from synoptic-scale waves (Branstator and Haupt 1998).

What is the observable imprint of linear planetary wave dynamics? In this linear view the atmospheric flow would consist of preferred flow patterns; any linear combination of these preferred patterns would again be a preferred flow pattern, and the probability density functions (PDFs) of these flow patterns would not have multiple extrema. If driven by additive white noise, these PDFs would be Gaussian. Another indication of linear dynamics is the form of mean tendencies in phase space. For a stable linear system these tendencies must be either pure exponential damping or a damped oscillation around the origin (Branstator and Berner 2005).

What are observable signs of nonlinear planetary wave dynamics? One possible sign of nonlinearity is deviation from Gaussianity in PDFs (Hansen and Sutera 1986; Kimoto and Ghil 1993a; Cheng and Wallace 1993; Corti et al. 1999; Christiansen 2005; Berner and Branstator 2007). Another possible sign of nonlinearity is the structure of the mean values of tendencies in a low-dimensional phase space as studied by Berner (2003), Selten and Branstator (2004), and Branstator and Berner (2005). Additional suggestions of the existence of such nonlinear tendency signatures have been presented by Haines and Hannachi (1995) and Hannachi (1997a, b), who found multiple local minima in the tendencies of quasigeostrophic models after truncating to only a few EOF dimensions, and by Kimoto and Ghil (1993b) and Crommelin (2003), who found evidence of preferred regime transition routes in nature and a barotropic model, respectively.

The huge body of work on planetary wave dynamics hints to the fact that both linear and nonlinear processes are contributing to planetary wave behavior. Berner (2003, 2005), Branstator and Berner (2005), and Berner and Branstator (2007) found non-Gaussian features and signatures of both linear and nonlinear planetary wave dynamics in a study of a very long GCM integration. They detected dominant linear signatures in certain planes of phase space, and distinct nonlinear signatures in others. Consistent with the above finding is that in some studies certain teleconnection patterns like the Pacific–North America (PNA) pattern are mainly driven by linear processes (Cash and Lee 2001; Feldstein 2002; Franzke and Feldstein 2005), whereas other patterns like the North Atlantic Oscillation (NAO) are mainly driven by nonlinear processes (Feldstein 2003; Vallis et al. 2004; Franzke and Feldstein 2005). These teleconnection patterns are prominent manifestations of planetary wave dynamics.

To further a unified view of planetary wave dynamics, in this paper we identify and quantify which dynamical processes are responsible for driving planetary waves. In the next section we will derive a general theoretical framework for doing so. This framework will enable us to diagnose in complex GCMs how important the contributions of the interactions among a highly truncated representation of the planetary waves are (referred to as resolved modes in the following). Some past studies have raised the possibility that the interaction among the planetary waves themselves produces a significant imprint on planetary wave behavior. The framework will also reveal how important the remaining degrees of freedom are (also referred to as unresolved waves). This general framework not only quantifies the contributions from the unresolved waves, it also reveals the functional form of their primary interactions with the resolved modes. An understanding of this interaction is crucial in the development of systematic low-order models of the planetary waves (Majda et al. 1999, 2003, hereafter MTV99, MTV03, respectively; Franzke et al. 2005; Franzke and Majda 2006; Crommelin and Vanden-Eijnden 2006a, b; Selten 1995; Kwasniok 1996) and will be fundamental in every unified theory of planetary wave dynamics.

Studies of very long integrations of various GCMs show that the PDFs of the dominant planetary waves (measured by the leading EOFs) are nearly Gaussian and do not have multiple peaks (Hsu and Zwiers 2001; Berner and Branstator 2007; Franzke and Majda 2006). On the other hand, these PDFs do have detectable departures from Gaussianity, including weak skewness, and in the case of bivariate PDFs multiple radial ridges of enhanced density (Berner and Branstator 2007). How important are these (small) departures from Gaussianity? Are these small departures able to produce climatically significant signatures in the planetary waves? The theoretical framework, derived below, will help us in answering these questions by allowing us to derive for arbitrary dynamical systems exact mean phase space tendency equations for the case of purely Gaussian modes. Any departure from the behavior that is indicated by this equation when applied to a particular system can then be attributed to non-Gaussianity.

In this study we try to further a unified view of planetary wave dynamics and do this with a hierarchy of approaches. First, in section 2, we develop the general theoretical framework that is alluded to above. Next we investigate the low-frequency dynamics of a simple stochastic toy model that has a strong nonlinear signature in the form of a double swirl in the mean phase space tendencies of its low-frequency variables in section 3. This toy model allows us to test the predictions of the framework in a direct way by changing model parameters. In section 4 we use a prototype GCM to investigate its mean phase space dynamics and diagnose its interactions with unresolved degrees of freedom. This model is closely related to the model in which Selten and Branstator (2004) detected nonlinear planetary wave mean phase space tendencies. Furthermore, the mean tendencies of this model are very similar to those in the comprehensive National Center for Atmospheric Research (NCAR) Community Climate Model, version 0 (CCM0) GCM, where nonlinear phase space tendencies were found to be important (Berner 2003; Branstator and Berner 2005). We end this paper with a discussion of our results in section 5.

## 2. Mean phase space tendencies and their information

In this section we develop a general framework for diagnosing contributions to the mean values of state-dependent tendencies in a low-dimensional subspace of complex geophysical systems. As has been already alluded to in the introduction, the mean phase space tendencies in some GCMs show distinct nonlinear signatures in certain planes, which are spanned by its leading EOFs (Berner 2003; Branstator and Berner 2005; Selten and Branstator 2004). Those leading EOFs also show only weak deviations from Gaussianity, mostly in the form of weak skewness, kurtosis, and, in the case of joint PDFs, multiple radial ridges of enhanced density (Berner and Branstator 2007; Franzke and Majda 2006). We are especially interested in revealing how much the self-interaction among the modes spanning those planes (resolved modes) and how much the unresolved modes contribute to these mean phase space tendencies and what effect the observed small deviations from Gaussianity have. To highlight the dynamical impact of the deviations from Gaussianity we consider the case of a system with purely Gaussian statistics. For this case the mean phase space tendency equation simplifies considerably.

### a. Total mean tendency equation

We derive a general form of the mean phase space tendencies for dynamical equations with quadratic nonlinearity of the same form as the equations of motion. We write these equations symbolically as

where **x** denotes the state vector, *F̃* a constant forcing, *L̃* a linear, and *B̃* a quadratic nonlinear operator. Because we are mainly interested in the dynamics of a low-order subspace spanned by only a few high-amplitude, large-scale patterns we decompose **x** into EOFs of model data in some appropriate metric [·, ·] [e.g., the traditionally used squared norm of geopotential height or streamfunction, or a more dynamically based total energy norm (Selten 1995)],

where *a _{i}* denotes the

*i*th principal component,

**e**

*the*

_{i}*i*th EOF, 〈

**x**〉 the time mean state, and

the variance, where *λ _{i}* denotes the

*i*th eigenvalue of the EOF decomposition. The decomposition of model output into empirical basis functions is standard practice in low-order modeling (Schubert 1985; Selten 1995; Kwasniok 1996; Branstator and Haupt 1998). By inserting (2) into (1) we obtain the equivalent equations of motion in EOF space:

where *F _{i}*,

*L*, and

_{ij}*B*are the constant forcing, linear, and nonlinear interaction tensors, respectively.

_{ijk}### b. Conditional mean tendency equations

To derive the mean tendency equation for the resolved modes we split the state vector **a** into resolved modes *α* and unresolved modes *β*. The variables *α* and *β* are vectors where the subscript *j* denotes the *j*th vector element. Now the part of (4) for the tendencies of the resolved waves *α* can be written as

where the indices either belong to the set of resolved (*I _{R}*) or unresolved (

*I*) modes. For simplicity of notation we symmetrize

_{U}*B*about the last two indices, so that

_{ijk}*B*=

_{ijk}*B*. Equation (5) gives a full description of the dynamics of the resolved modes, where (5a) are the tendencies stemming from the self-interaction among the resolved modes (bare truncation) and (5b) are the tendencies from interactions involving unresolved modes. The remainder of this section is devoted to estimating the relative importance of those two contributions to the mean tendencies of the resolved modes, as in the mean phase space tendency studies listed in the introduction.

_{ikj}We start with the complete PDF of the system (4), which we denote by *p*(*α*, *β*). By assuming ergodicity the expectation or average of an arbitrary function *f* is given by

for **a** = (*α*, *β*). The conditional mean tendency equation involves utilizing the conditional mean probability densities

Thus, the conditional expectation of an arbitrary function *f* for a given value of *α* is given by

Such a conditional probability density can be gotten from *p*(*x*, *y*) by averaging so that

The conditional mean tendencies for the resolved modes *α* involve evaluating linear and quadratic nonlinear terms like

With (5) and (10) it is easy to compute a general formula for the *conditional mean tendencies* for the dynamics [(4)] in the resolved variables,

Note that the right-hand side of (11a) is the *bare truncation* restricted to interactions among the resolved modes, while (11b) and (11c) involve all of the conditional mean statistics 〈*β _{j}*/

*α*〉; the terms of (11c) are associated with

*multiplicative triad interactions*(leading to multiplicative noise in a MTV99 and MTV03 framework) of

*β*, and the terms in (11d) involve all of the conditional interaction statistics 〈

_{j}*β*/

_{j}β_{k}*α*〉 of second moments and include all of the

*additive triad interactions*(leading to additive noise in a MTV99 and MTV03 framework) of

*β*and

_{j}*β*. Thus, Eq. (11) 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 degrees of freedom.

_{k}### c. Explicit evaluation for a Gaussian PDF

Now we simplify the above conditional mean tendency equation for purely Gaussian EOF modes. The PDFs in GCMs are nearly Gaussian (Hsu and Zwiers 2001; Berner and Branstator 2007; Franzke et al. 2005; Franzke and Majda 2006; Majda et al. 2006), and there are many geophysical models without damping and forcing that exactly satisfy these assumptions such as barotropic flow on the sphere with topography (Salmon 1998; Carnevale and Frederiksen 1987; Majda and 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

where *p*^{G}_{i}(*α*_{i}), *p*^{G}_{j}(*β*_{j}) are Gaussian distributions with mean zero. This is true because the EOFs are uncorrelated. With the Gaussian assumption and the factorization of the PDF in (12) we have

Thus, in the *Gaussian* case, the conditional mean tendency equation simplifies to

Note that contributions from (11b) and (11c), that is, linear coupling and multiplicative triad interactions, are identically zero. Thus, in this Gaussian case the conditional mean tendency equation recognizes bare truncation [(11a)] and a constant forcing from additive triad interactions [(11d)]. Note that the conditional averages in (13) only require that the PDF *p*(*α*, *β*) = *p*(*α*)*p*(*β*) factorizes into the PDFs of the resolved and unresolved modes, and that the unresolved modes factor into Gaussians.

## 3. Mean phase space tendencies of a simple stochastic climate model

### a. Stochastic climate model

To test the predictions of the framework that are derived in the previous section we utilize a four-mode stochastic climate model of the kind put forward in Majda et al. (2005, their chapter 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 degrees of freedom. Such simple toy models allow for the efficient exploration of the whole parameter space, which is impossible to conduct with GCMs. Thus, we are able to test the predictions of the above framework with direct model experiments by switching certain terms on and off rather than relying on diagnostic methods.

While this model is not rigorously derived from a geophysical flow model (e.g., a 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 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, **y**. In realistic models there would be many fast modes representing, for example, synoptic weather systems or convection. To mimic their combined effect, we include damping and stochastic forcing −(*γ*/ɛ)**y** + (*σ*/ɛ)**Ẇ** in the equations for **y** where **W** denotes a Wiener process. The motivation for this approximation is that these fast modes are associated with turbulent energy transfers and strong mixing, and we do not require a more detailed description because we are only interested in their effect on the slowly resolved modes. The two fast modes carry most of the variance in this model. 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 **x** as the resolved modes and the two fast modes as the unresolved modes **y**. Therefore, our toy model has the following form:

To ensure energy conservation of the nonlinear operator the coefficients have to satisfy the following relations: *b*_{123} + *b*_{213} + *b*_{312} = 0 and *c*_{134} + *c*_{341} + *c*_{413} = 0. In this particular setup the slow and 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 MTV99 and MTV03 framework [see also Eq. (11)], so we refer to it as a multiplicative triad; whereas the nonlinear interactions involving

*c*produces additive noise in the MTV99 and MTV03 framework [see also Eq. (11)], and so we refer to it as an additive triad.

_{ijk}To systematically vary the strength of various dynamical processes, (15) contains four parameters that control the strength of various interactions: *μ* controls the bare truncation, *θ* the linear coupling, *κ* the additive triad, and *λ* the multiplicative triad interactions.

### b. Model integration details

To explore the parameter space we created a 1000 member ensemble by integrating the model for 10^{4} time units for each of 1000 different initial conditions that were chosen randomly. This has been done for every parameter set we explored. To integrate the model we have used a fourth-order Runge–Kutta scheme for the deterministic part and Euler forward differencing for the stochastic part. We have used a time step of 10^{−3} time units and saved output every ¼ time unit. The following values have been used for the control integration: *b*_{123} = 0.25, *b*_{231} = 0.25, *b*_{312} = −0.5, *c*_{134} = 0.25, *c*_{341} = 0.25, *c*_{413} = −0.5, *L*_{13} = −1.0, *L*_{24} = 1.0, *L*_{12} = *L*_{21} = 1.0, *a*_{1} = 1.0, *a*_{2} = −1.0, *d*_{1} = −2/10, *d*_{2} = −1/10, *F*_{1} = −0.25, *F*_{2} = *F*_{3} = *F*_{4} = 0.0, *γ*_{1} = *γ*_{2} = 1.0, and *σ*_{1} = *σ*_{2} = 1.0.

In all integrations we set the time-scale separation parameter ɛ = 1. This value is consistent with estimates from long integrations of barotropic and quasigeostrophic models with realistic climatology (Franzke et al. 2005; Franzke and Majda 2006). However, all integrations have been carried out also with values of 0.5 and 0.1 for ɛ. These integrations reveal qualitatively very similar results to those for ɛ = 1.

### c. Mean phase space tendencies

Now we investigate the mean phase space tendencies of the low-frequency variables of the stochastic climate model (15). We do this for four different values of the constant forcing parameter *F*_{1}; namely, −0.5, −0.25, −0.18, and −0.1. Those parameter values correspond to different regimes of (15), as will be shown in section 3d. Figure 1 shows that in all four cases the mean phase space tendencies have a very pronounced double-swirl structure. This clearly indicates a strongly nonlinear behavior in the plane spanned by the low-frequency modes. Furthermore, the statistics of the **x** and **y** modes are very close to being Gaussian (not shown). Thus, this model is a perfect tool for numerically testing the predictions of the framework for the Gaussian case (14), according to which only the bare truncation effectively contributes to the nonlinear phase space signature.

To numerically test the predictions in detail results are presented from sensitivity experiments where certain parameters are changed. Figure 2 shows mean phase space tendencies where the strength of the bare truncation is reduced by decreasing the parameter *μ* in (15) in the case where *F*_{1} = −0.25. These results show clearly that by reducing the amplitude of the bare truncation, the double-swirl structures weaken and then disappear. We also study the case where only the amplitude of the nonlinear bare truncation coefficients is reduced. In this case the double-swirl structure also vanishes but the model retains a strong single-swirl structure, stemming from the linear skew-symmetric part (not shown). Because in both cases the PDF remains nearly Gaussian, both are consistent with the theory. We further test the theory by weakening the effect of the triad interactions. Whether we weaken *λ* (Fig. 3a) or *κ* (Fig. 3b), the double swirls remain, further confirming that the unresolved modes are not contributing to the swirls.

Next we investigate how strong the impact is of the linear coupling on the low-frequency tendencies. We have carried out many experiments where we have changed the value of *θ* and also the values of *L*_{13} and *L*_{24}, but we find only a moderate influence that mostly concerns the areal extent of the dynamically active portion of phase space (see, e.g., Figs. 3c,d). However the linear coupling never changes the structure of the double swirls.

To summarize, the numerical experiments show that the behavior of the simple stochastic climate model verifies the predictions of the theory for the Gaussian case. They show that for this simple toy model the bare truncation dynamics have a formative influence on the existence of the double swirls and that the unresolved dynamics do not contribute much to this structure in the mean phase space tendencies of the low-frequency variables.

### d. Fixed points of bare truncation

Our toy model is so simple that we can determine its fixed points and also the bare truncation phase space tendencies. It is dynamically interesting to see how these fixed points are related to the origin of the double-swirl structures.

To do so we set *λ* = *θ* = *κ* = 0 in (15). Thus, we are analyzing a two-dimensional system consisting of only the two **x** modes. A fixed-point analysis has been made in great detail in Majda et al. (2005, their chapter 3) and here we briefly summarize those points of their results that are relevant for this study.

When varying the climatological forcing *F*_{1} and *F*_{2}, as bifurcation parameters, the dynamical system undergoes qualitative changes. The full bifurcation diagram is shown in Majda et al. (2005, their Fig. 3.2). When setting *F*_{2} = 0 and varying *F*_{1}, three bifurcations occur on that line. Inside each regime we choose a representative value for *F*_{1}, namely, −0.5, −0.25, −0.18, and −0.1. For both *F*_{1} = −0.1 and *F*_{1} = −0.5 the system develops unique steady states, and for both *F*_{1} = −0.18 and *F*_{1} = −0.25 the system develops two stable steady states (see Table 1). The bifurcation analysis has been confirmed by Monte Carlo integrations in which, for each of these four values of *F*_{1}, we start from 10^{4} randomly drawn initial conditions and integrate the model forward for 10^{4} time units.

Furthermore, we now present the bare truncation phase space tendencies for all four forcings (Fig. 1, left column). These plots reveal two centers of rotation for the forcing values *F* = −0.18, −0.25, and −0.5. These centers do not necessarily correspond to stable fixed points. One center of rotation (*x ^{L}*) is close to the origin and the other (

*x*) occurs at larger values. For

^{H}*F*= −0.1 there is one center of rotation and a ridge of convergence (Fig. 1a, left column). In the basin of

*x*all trajectories slowly spiral clockwise toward the fixed point for values of

^{L}*F*= −0.1, −0.18, and −0.25, whereas outside this basin, all trajectories rapidly reach

*x*for

^{H}*F*= −0.18, −0.25, and −0.5. Due to the very slow motion toward

*x*and the normalizations of the arrows, the visual impression in Figs. 1a,b,c is misleading. For the forcing values

^{L}*F*= −0.1, −0.18, and −0.25

*x*is indeed a stable fixed point, whereas it is not a stable fixed point for

^{L}*F*= −0.5.

These results give insight into the relationship between the bare truncation, fixed points of the bare truncation, and the mean phase space tendencies. The locations of the stable fixed points are close to some of the centers of the double swirls of the mean phase space. But, not all centers of the double swirls are associated with stable fixed points. This suggests that the whole structure of the bare truncation phase space flow field is responsible for the mean phase space tendencies of (15).

## 4. Mean phase space tendencies of a prototype GCM

After successfully testing the predictions of the framework on a simple toy model, we apply our framework now to a prototype GCM that simulates key aspects of the observed circulation.

### a. Prototype atmospheric GCM

The prototype atmospheric GCM is a spectral three-level quasigeostrophic (QGT21L3) model originally developed by Marshall and Molteni (1993). It is used in this study in the same hemispheric setup as used by Franzke and Majda (2006). The QG equation can be written as

where *q* is potential vorticity (PV), *ψ* streamfunction, *D*(*ψ*) a linear operator that represents dissipative terms, *S* a constant PV source, and *J* the Jacobian operator. The model is triangularly truncated at total wavenumber 21. The constant forcing *S* is derived from the European Centre for Medium-Range Weather Forecasts (ECMWF) data in such a way as to approximate the observed Northern Hemisphere climatological mean state and variance (see Marshall and Molteni 1993 for a more detailed discussion). The model is integrated with a fourth-order Adams–Bashforth numerical scheme with a time step of 10 min. The state vector is saved every 6 h. The model is started from random initial conditions and is integrated for 500 000 days after an initial 1000-day spinup.

### b. Tendency equation

Unlike for the toy model it is impossible to switch on and off various interactions between resolved and unresolved modes in such a model. To determine contributions to mean phase space tendencies we have to calculate the various tendency terms from model output. To do so we write (16), with PV defined in the same way as in Marshall and Molteni (1993) as a streamfunction equation:

where 𝗤 = (∇^{2} + 𝗘)^{−1}, ∇^{2} is the Laplace operator, and 𝗘 is a vertical discretization operator

with the understanding that 𝗘 acts equally on all spectral modes of a vertical level (see Marshall and Molteni 1993 and Franzke and Majda 2006 for more details).

Now we decompose the streamfunction into

where the overbar defines a time mean and the prime the deviation from the time mean. Next we insert this decomposition into the streamfunction Eq. (17),

with

Here *χ _{L}* denotes the linear tendencies, including the planetary and relative vorticity advection, the topographic forcing, and dissipative processes, whereas

*χ*denotes the transient eddy flux. We also have to consider tendencies stemming from the time mean state and the constant forcing

_{NL}*S*. These result in a constant contribution:

### c. Mean phase space tendencies

In this section we describe the mean phase space tendencies in various planes of the QGT21L3 model phase space. These planes are spanned by the leading total energy norm EOFs (Franzke and Majda 2006). But first we focus on the joint PDFs. Figure 4 shows that the joint PDFs are very close to Gaussian. They are weakly skewed and have radial ridges in some planes but not multiple peaks. These PDFs alone do not show signs of strong nonlinear planetary wave dynamics.

In Fig. 5 we display the mean phase space tendencies in various planes. In the plane that is spanned by the first two EOFs there is a clear double-swirl structure visible (Fig. 5a). The physical flow structure corresponding to the two centers of the swirls is similar to opposite phases of the annular mode as in other studies (Selten and Branstator 2004; Branstator and Berner 2005; Berner and Branstator 2007). Also in the plane that is spanned by EOF3 and EOF4 a clear double-swirl structure is visible (Fig. 5d). These structures strongly indicate nonlinear planetary wave behavior. Such double-swirl structures are less dominant but still present in the planes that are spanned by EOF1 and EOF3, and EOF1 and EOF4 (not shown). The planes that are spanned by EOF2 and EOF3, and EOF2 and EOF4 display a single-swirl structure (Figs. 5b,c). As mentioned earlier, single-swirl structures can be produced by linear dynamics, whereas the double-swirl structures are necessarily imprints of nonlinear dynamics. These results also show that even though nonlinear dynamics are evidently present in some planes the corresponding PDF in those planes can still be nearly Gaussian.

These preliminary results indicate that in certain phase planes the planetary wave dynamics are potentially predominantly linear, whereas in other planes there are strong nonlinear signatures. In the following we will mainly focus on the plane that is spanned by EOF1 and EOF2 to see if these structures are created mainly by the bare truncation dynamics, that is, the interaction between EOF1 and EOF2, or if the unresolved EOF modes play a crucial role. Moreover, we will determine through which physical processes they interact with the resolved modes.

For reference we display in Fig. 6a the mean phase space tendencies without the mean forcing (*χ _{F}*), which is not a function of state. As can be seen in Fig. 6a the double-swirl structure is still present even though it is less dominant. Also, the structure has experienced a rotation compared with Fig. 5a. It is this structure in Fig. 6a that we are using below for comparison in the budget analysis in which we ignore the constant mean state terms (

*χ*).

_{F}In Figs. 6b,c we display the bare truncation phase space tendencies; as in (11a) these are the tendencies one would get if only EOF1 and EOF2 interact with each other and not with any of the other EOFs. This gives insight into how important the interaction between EOF1 and EOF2 with all the other EOFs is. Note that in Fig. 6, as in other plots of phase space tendencies, the shading scale varies from panel to panel. As can be seen in Fig. 6b the total bare truncation tendency induces an outwardly directed tendency with a weak rotational part. The bare truncation nonlinear interaction [the nonlinear part of (11a), Fig. 6c] shows a structure that resembles a double swirl, but with a much weaker and different structure from the mean phase space tendencies (Fig. 6a). The comparison of the bare truncation tendencies with the model tendencies (Fig. 6a) clearly reveals that the unresolved modes make a substantial contribution in this model.

### d. Phase space tendencies induced by the unresolved modes

Next we investigate the contributions of the unresolved modes in more detail and reveal which dynamical processes are involved. We are especially interested to reveal the impact of the unresolved modes on the mean phase space tendencies in the EOF1–EOF2 plane.

In Fig. 7 the phase space tendency contributions from the linear coupling [(11b)], the additive triad interactions [(11c)], and the multiplicative triad interactions [(11d)], which are conditioned on the first two EOFs, are displayed. We display the partial sums of the bare truncation [(11a)] and linear coupling [(11b)] in Fig. 8a, the multiplicative triads [(11c)] in Fig. 8b, and additive triads [(11d)] in Fig. 8c. These figures indicate that all three budget terms are needed to recover the full tendency structure; of special importance for recovery of the double-swirl structure are the additive and multiplicative triad interactions, because those tendencies produce nonlinear signatures. These results clearly show the important role the unresolved waves play in the planetary wave dynamics and in imprinting the nonlinear signature. Furthermore, these results indicate that non-Gaussian effects are important for the nonlinear signature in this model.

To estimate the magnitude of deviations from Gaussianity we calculate the skewness and kurtosis of the total energy norm EOFs (Fig. 9). They indicate that several of the first 30 EOFs experience some moderate skewness, with magnitudes of at most 0.25, and also moderate kurtosis. The high values of skewness and kurtosis for the trailing EOFs (EOF 400 and higher) are due to numerical inaccuracies or sampling. Because these trailing EOFs explain very small amounts of variance, their contribution to the mean tendencies of the planetary wave dynamics is likely negligible.

Our main focus was on the EOF1–EOF2 plane. We now want to briefly discuss some other planes. By comparing the other plane that shows a double-swirl structure (EOF3–EOF4, Fig. 5d) with its bare truncation phase space tendency (Fig. 10c), we see the same behavior as for the EOF1–EOF2 plane in that the bare truncation tendency is dominated by a rotational motion in phase space; therefore, the nonlinear signature has to be due to the unresolved degrees of freedom. The phase space tendencies of the planes composed of EOF2 and EOF3 (Fig. 5b), and EOF2 and EOF4 (Fig. 5c) induce a rotation in the opposite direction to their respective bare truncation dynamics (Figs. 10a,b). This indicates that even for linear signatures in the leading phase space planes the unresolved degrees of freedom play a crucial role in the dynamics of the planetary waves in this model.

## 5. Discussion

In this study we investigate the origin of nonlinear planetary wave behavior and put forward a unified view of its dynamics by utilizing theory, a simple stochastic climate model, and a prototype atmospheric GCM. We focus in our study on nonlinearities evident in the mean phase space tendencies of planetary waves because recent studies show that in prototype (Selten and Branstator 2004) and comprehensive GCMs (Branstator and Berner 2005) these tendencies experience a double-swirl structure in some of their dominant planetary waves, indicating the strongly nonlinear behavior of these low-frequency waves. In other planes of phase space the imprint of planetary wave dynamics is rather linear, as shown by Branstator and Berner (2005) and in this study.

Building on the ideas by Charney and DeVore (1979) for a very low-dimensional system, Haines and Hannachi (1995) proposed that nonlinear signatures in phase space tendencies in low-dimensional subspaces of high-dimensional systems could be caused by nonlinearities in the interactions within the subspace. On the other hand, it is possible that interactions involving modes that are not resolved in the low-dimensional subspace may be important. To sort out these possibilities we introduce a framework in section 2 for systematically identifying the relative importance of various interactions between resolved and unresolved modes in producing tendency signatures of interest. One important fact that this framework makes clear is that for systems with Gaussian statistics, the only interactions that can produce nonlinear phase space signatures in low-dimensional subspaces like those seen in earlier studies are interactions between modes of the subspace.

Note that our approach is very different from the approach of Friedrich et al. (2000), Berner (2003, 2005), Sura (2003), Sura et al. (2005), and Crommelin and Vanden-Eijnden (2006a, b) for studying the behavior of a high-dimensional system in a reduced space. These studies approximate the reduced space behavior by directly estimating the drift and diffusion coefficients for a Fokker–Planck equation from data. In their approach the drift contains not only the bare truncation interactions, but also corrections representing the impact of the unresolved degrees of freedom on the resolved modes. Therefore, in these studies an attribution of the causes of observed nonlinear signatures to interactions among resolved and/or to interactions with unresolved modes is not made.

We first test the predictions of our framework on a simple stochastic climate model. This simple model shows a double-swirl motion in its low-frequency subspace, which is qualitatively similar to the prototype GCM. We find that the simple stochastic climate model follows the predictions of the framework closely. The model statistics are very nearly Gaussian and the swirls result from the interactions between modes resolved by the subspace. Furthermore, we find that the centers of rotation of the double swirls of the mean phase space tendencies do not necessarily correspond to stable fixed points of the reduced dynamics of the low-frequency subspace.

We next apply the framework to a prototype GCM whose mean phase space tendencies in low-order subspaces are similar to those by Branstator and Berner (2005) in that certain two-dimensional subspaces are dominated by nonlinear features while others are basically linear. Our results show that the unresolved modes contribute strongly to the mean phase space tendencies of the planetary waves. The unresolved modes interact with the planetary waves by both linear and nonlinear processes. This indicates that previous approaches to understanding planetary wave dynamics in terms of highly truncated planetary wave models are inadequate. And, of course, the presence of nonlinear signatures alone means that the common approach of approximating planetary wave dynamics with linear models driven by additive white noise is also inadequate. The importance of the unresolved modes also implies that, though weak, non-Gaussianities in the distribution of planetary wave states are essential for the production of the GCMs’ nonlinear mean phase space tendencies. Our results imply that in order to understand and model low-frequency planetary waves, procedures are needed that systematically account for the symbiotic relationship between the planetary and unresolved waves. Such systematic procedures have been proposed recently (MTV99, MTV03), and have been successfully applied to a barotropic model (Franzke et al. 2005) and a prototype atmospheric GCM (Franzke and Majda 2006) with realistic climates. The results of our study have important implications for low-order modeling of planetary waves and are consistent with the low-order stochastic modeling study by Franzke and Majda (2006). They highlight the crucial role that linear and nonlinear deterministic correction terms play, as well as both additive and multiplicative noises.

Our results also indicate that of the energetic planetary wave structures most of the non-Gaussianity resides in the first 30 or so EOFs. This suggests that these 30 EOFs are the primary contributors to the planetary wave dynamics of the leading two EOF modes and opens the possibility of systematically deriving a low-order model for the planetary waves involving, at most, 30 modes. This possibility will be explored in the future by the authors.

Our study also shows that one has to be careful in translating the ideas gained from simple models and idealized theoretical reflections to complex models. In contrast to our simple toy model, the prototype atmospheric GCM has a much more complex behavior of the planetary wave dynamics in that the unresolved degrees of freedom have a dominating impact on the mean phase space tendencies of the leading two planetary wave modes. While idealized theoretical considerations (in our case the assumptions of purely Gaussian statistics) and simple toy models are helpful tools that allow insight into otherwise overwhelmingly complex problems (Palmer 1999), such ideas should be tested on a hierarchy of models with increasing complexity as already proposed by Hoskins (1983) and Held (2005).

## Acknowledgments

The authors acknowledge generous support from the National Science Foundation under Grants NSF-CMG DMS-0222133 and DMS-0456713 and the Office of Naval Research under Grants ONR N00014-96-0043 and N00014-05-1-0164.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX

#### Estimating Conditional Mean Phase Space Tendencies

Because the estimation of the integrals over the probability distributions in (10) are computationally challenging, we pursue a different equivalent way to estimate the terms in (11).

First, we define the sum of the unresolved modes as

Now the tendencies of (11b) can be estimated in an appropriate metric as

the tendencies of (11d) can be estimated in an appropriate metric as

and the tendencies of (11c) are estimated as

## Footnotes

* Current affiliation: National Center for Atmospheric Research, Boulder, Colorado

*Corresponding author address:* Dr. Christian Franzke, National Center for Atmospheric Research, P.O. Box 3000, Boulder, CO 80307. Email: franzke@ucar.edu