## Abstract

A low-order deterministic qualitative model is formulated in order to simulate extratropical low-frequency variability. This deterministic model is based on a filtering of the potential vorticity equation on the 315-K isentrope and a projection onto its leading empirical orthogonal functions. The model has an empirical formulation, and the feedback of unresolved scales is taken into account. The model building procedure is novel, since it is not based on a severe truncation of the physical evolution equations but on an empirical analog averaging of each relevant dynamical process. It can be applied to any geophysical system for which long observational data series are available.

The model is used to diagnose weather regimes with its multiple equilibria and intraseasonal oscillations as periodic orbits. These equilibria result from the balance between large-scale advection, transient feedback, and residual forcing. The authors analyze their forcing budgets and show in particular that transient feedback tends to amplify and advect upstream the regime anomaly patterns. However, the key forcing turns out to be the large-scale advection, since the other forcing terms only reshape the regime anomalies. The maintenance of observed intraseasonal oscillations is also examined by means of forcing budgets. Results show that large-scale advection and transient feedback are also key dynamical factors in the maintenance of their life cycle.

Finally, the low-order model is integrated and qualitatively simulates two of the three identified oscillations, those with periods of 70 and 28 days. The intraseasonal oscillations show up as unstable periodic orbits in the low-order model. This indicates that these oscillations are mostly driven by the internal dynamics of the extratropical atmosphere.

## 1. Introduction

### a. The use of low-order models

Over the last 30 years or so, the theory of nonlinear dynamical systems has facilitated the understanding of major aspects of climate dynamics. The seminal work of Lorenz (1963) showed that a physical system can be at the same time deterministic and unpredictable. Although the three-variable model developed therein was fairly unrealistic, Lorenz’s result was rather general and particularly useful for our terrestrial climate. More than a decade later, Charney and DeVore (1979) also used a simple dynamical system, albeit nonchaotic, to demonstrate the existence of multiple extratropical atmospheric equilibria, given a fixed boundary forcing. Despite the lack of realism of their model and of the equilibria themselves, this work shook the traditional representation of the atmospheric flow as a basic state superimposed with small perturbations. The quest for multiple equilibria, using severely truncated models was also extended to paleoclimate dynamics (LeTreut and Ghil 1983) and to the ocean system (Jiang et al. 1994).

The theory of dynamical systems also allows a neat representation of nearly periodic phenomena, often called “oscillations,” as periodic orbits (Guckenheimer and Holmes 1983). These orbits are usually unstable, which explains that the system’s behavior is not purely periodic in the long run, but exhibit intermittent “spells” of oscillations instead (Plaut and Vautard 1994). Therefore, one is lead to study these orbits with values of the system’s parameters for which weakly nonlinear theory can be applied (Pedlosky 1987), that is, in the vicinity of *Hopf bifurcations.* As an example of application, Hopf bifurcations and their subsequent periodic orbits were proposed as prototypes for oscillatory topographic instabilites in the extratropical atmosphere (Yoden 1985) and intraseasonal oscillations (Jin and Ghil 1990).

Theoretical problems of a different nature were also tackled using low-order models, such as the existence of a “slow manifold” for the extratropical atmosphere (Lorenz 1980; Vautard and Legras 1986; Jacobs 1991) and for testing methodologies, such as in data assimilation (Gauthier 1992; Miller et al. 1994; Pires et al. 1996). The common and essential feature of the models used in the above studies is their simplicity. Above all, these are *qualitative* models: they do not attempt to represent most of the interactions of the underlying complex system, only the key features.

All previous formulations of low-order dynamical systems consisted in simplification and severe truncation, to a reduced set of functions, of the complex flow evolution equations. Recently, there have also been several successful attempts in reducing the dimension of the phase space by projecting model equations onto a reduced set of empirical orthogonal functions (EOFs) (Selten 1993; Mundt and Hart 1994; Kwasniok 1996). These models show increased realism with respect to truncated models based on projection onto harmonic functions. We propose here a different approach to low-order model construction, based on an *empirical* calculation of the flow tendency, combined with projection onto a reduced set of EOFs. The low-order model is therefore constructed from long records of flow maps and their tendencies. The advantage of an empirical formulation is the hope to better fit real data, but its drawback is the difficulty of using it for sensitivity studies: it can only be used for diagnosis of the observed system.

### b. Extratropical low-frequency variability

Yet a general formulation of low-order models is not the only goal of the present study. The particular problem of interest here is the maintenance of the atmospheric extratropical, low-frequency, intraseasonal variability (5–100 days). At these timescales, the extratropical atmosphere is chaotic and unpredictable but can still be described accurately by a limited number of variables. Many investigators showed, using principal component analysis (Preisendorfer 1988), its rotated variant (Richman 1986), or its complex formulation (Horel 1981), that most of the variance is described by a much smaller number of degrees of freedom than any general circulation model needs in order to simulate realistic behavior.

Low-frequency extratropical variability therefore calls for low-order modeling. The zero-order requirement for such a model to be qualitatively realistic is to reproduce climatic means. The first-order requirement is a correct simulation of the so-called extratropical weather regimes (Reinhold and Pierrehumbert 1982), which account for a major part of the variability. Here we impose the condition that the model reproduce other regular phenomena such as intraseasonal oscillations, as identified by many investigators (Branstator 1987; Kushnir 1987; Knutson and Weickmann 1987; Ghil and Mo 1991; Plaut and Vautard 1994).

The major modeling difficulty comes from the interaction between large scales and unresolved synoptic scales. This interaction leads in particular to an equilibration of the weather regimes by the baroclinic transients (Reinhold and Pierrehumbert 1982; Vautard and Legras 1988). Such an equilibration can also exist for the maintenance of low-frequency oscillations but has never been documented previously. Our study investigates this possibility. This equilibration picture would lead one in particular to consider that relevant large-scale phenomena can only be represented by fixed points and periodic orbits of a low-order model if a correct parameterization of transient feedback is included. Other processes are to be taken into account if one aims at reaching some realism. Such is the case of all physical processes as radiative transfer, diabatic processes, and boundary forcing. Our general model formulation incorporates these neglected processes empirically.

### c. A general construction of low-order models

Let us assume that the underlying physical system is defined by the set of variables grouped into a vector *Z* of (possibly infinite) dimension *n.* The evolution equations are given by the dynamical system

The projection of this system onto a low-order phase space separates the orginal state vector into two vectors, *Z* = *X* + *Y,* where *X,* called hereafter the “resolved state vector,” lies in a phase space of dimension *m* much smaller than *n* and *Y* represents the “unresolved variables.” According to the above discussion, the vector *X* can contain the first few leading EOF coefficients of the *Z* vector, called hereafter the principal components (PCs). The following methodology, however, does not depend on the fact that basis functions are specifically EOFs.

The projection of Eq. (1) onto the resolved variables is

This evolution equation is not closed since *Y* may take different values for a given *X.* Hence, the tendency *G*(*X, Y*) can be considered as a random variable whose probability density distribution depends on *X.* The instantaneous value of *G*(*X, Y*) can be split into a mean 〈*G*〉(*X*), the expectation of this distribution, and a random fluctuation *R*:

The function 〈*G*〉(*X*) therefore represents the conditional average, given *X,* of its instantaneous tendency. It will be called hereafter the *deterministic* part of the tendency. The long-term behavior of the original system (1) leads to a “climatic” distribution of variables *X* and *Y,* described by the density distribution *ρ*(*X, Y*). The deterministic tendency can be expressed in terms of this density:

The fluctuation term *R* is stochastic, and its distribution also depends a priori on *X*, but has a zero expectation. Our approach of closure of (2) consists of omitting the residual term *R*. The low-order deterministic model is therefore

This formulation is different from the classic low-order truncation of complex flow equations, where instead the unresolved scales would be set to zero but parameterized by addition of a new term *D*(*X*), which in general is dissipative, leading to

In Eq. (5), unresolved-scale feedback is included implicitly.

The low-order model (5) can now be integrated starting from any initial condition *X*. By removing the residual, we have presumably removed a large source of variability; hence we do not expect that (5) will be a realistic model. Particularly in the case where the conditional distribution of *G*(*X, Y*) is multimodal, its expectation is a poor representation of the “deterministic tendency.” Instead, we could identify clusters of tendencies and estimate the deterministic tendency as a centroid of one of these clusters. Our purpose here is to examine the simplest approach, that is, the expectation, and to investigate what qualitative features such a simplified model is able to capture.

The calculation of the deterministic part of the tendency [Eq. (4)] is impossible to formulate analytically and is expressed here in a fully empirical fashion: 〈*G*〉(*X*) is calculated using large samples of recorded values of *X* and *Y* and of their tendencies. At a given stage of the integration of the model, *X* only is known, and 〈*G*〉(*X*) is estimated by averaging the observed instantaneous tendencies over a set of selected analogs of *X* found in the data record. A technical description of this calculation has been described by Vautard (1990) and is reviewed in section 3.

### d. Forcing budgets

The system (5) is a low-order dynamical system, for which one can easily determine fixed points (equilibria), in the sense of Charney and DeVore (1979) or Legras and Ghil (1985), or periodic orbits (Guckenheimer and Holmes 1983; Ghil and Childress 1987). Vautard and Legras (1988) used the above formulation in order to define the weather regimes of a quasigeostrophic two-level channel model. The resolved variables were the PCs of the large-scale flow. The above methodology was also used by Michelangeli et al. (1995) in order to identify weather regimes in the Atlantic and Pacific domain, and by Molteni (1996) to assess the role of transient feedback in the maintenance of phase-space probability density peaks and clusters.

When the tendency *G*(*X, Y*) can be expanded into several terms, say *G*_{1}(*X, Y*), *G*_{2}(*X, Y*), . . . , corresponding to different forcing mechanisms, each of these terms can be submitted to the above closure procedure separately. Their resulting deterministic components 〈*G*_{1}〉(*X*), 〈*G*_{2}〉(*X*), . . . , which add up to zero when *X* is an equilibrium of Eq. (5), represent the role of these forcing mechanisms in the maintenance of the equilibrium.

Such a decomposition permits, therefore, an analysis of forcing budgets. In Michelangeli et al. (1995), this decomposition was difficult since the variable used was the geopotential height. Yet, using quasigeostrophic assumptions, Lau and Holopainen (1984), Mullen (1987), and Lau (1988) were able to carry out geopotential tendency budgets of nonequilibrated anomalies. Here, a low-order model is formulated in terms of potential vorticity (PV, hereafter) on an isentropic surface (Hoskins et al. 1985). The main advantage of using PV is the simplicity of the equations governing its evolution. The resulting distinction of transient feedback, advection, and diabatic forcing is particularly easy (Hoskins 1991). We attempt here to assess their roles in the maintenance of weather regimes and intraseasonal oscillations.

### e. Outline

In section 2, we describe the data used. Section 3 is devoted to the technical description of the low-order model’s construction. In section 4, the maintenance of the weather regimes is examined. We examine the role of various forcing mechanisms, such as transient feedback, in their maintenance. In section 5, forcing budgets over the intraseasonal oscillations life cycle are carried out. Section 6 is designed to analyze integrations of the low-order model itself and to determine what part of the variability is captured. Section 7 contains a summary and a short discussion.

## 2. The potential vorticity dataset

The atmospheric data used in this study are twice-daily maps of the Northern Hemisphere streamfunction Ψ and potential vorticity on the 315-K isentropic level. This dataset, described and validated by Brunet et al. (1995), was obtained from analyzed geopotential heights and temperatures from the National Centers for Environmental Prediction (NCEP), formerly the National Meteorological Center, on nine pressure levels and covers the period 1963–88. The procedure followed by Brunet et el. (1995) to obtain the potential vorticity and streamfunction maps consisted of four steps.

Horizontal interpolation from the original geopotential and temperature maps on the NCEP octogonal grid to a T31 Gaussian grid, with the assumption that both temperature and geopotential are symmetrical fields with respect to the equator.

Vertical interpolation of the same quantities to the target isentropic level. This procedure was performed at each grid point of the T31 Gaussian grid.

Calculation of the Montgomery potential on the isentropic level and resolution of the nonlinear gradient-wind balance, in order to obtain the streamfunction.

Calculation of potential vorticity using the hydrostatic approximation.

By comparison with PV maps calculated directly from winds analyzed by the European Centre for Medium-Range Weather Forecasts, Brunet et al. (1995) showed that the nonlinear gradient-wind balance approximation provides more realistic PV maps than that of the simple linear balance. In particular, the latter was shown to overestimate troughs and underestimate anticyclonic structures.

The height of the 315-K level isentropic surface varies significantly in latitude (typically 700 hPa in the Tropics and 200 hPa at high latitudes). It usually crosses the tropopause near 45°, forming a band of large PV gradients due to the abrupt variation of the pseudodensity field. This band demarcates a polar area within which stratospheric polar air is confined. Owing to the large horizontal gradients in the tropopause region, we anticipate that the role of transient motions will be enhanced on this isentropic level at these latitudes.

Throughout this study, we focus on the North Atlantic sector, where low-frequency fluctuations, taking the form of blocking highs, strong zonal flows, or other patterns, are observed (Michelangeli et al. 1995). The North Atlantic sector is defined as the domain lying between 30° and 70°N latitude and 80°W and 40°E longitude.

## 3. The low-order potential vorticity model

### a. Potential vorticity evolution equation

Potential vorticity *P* is a material invariant when the (stratified) flow is adiabatic and does not undergo friction or other mechanical forces (Hoskins et al. 1985). The evolution equation of PV on an isentrope therefore reads

where *F _{p}* contains diabatic and mechanical forcing. The advantage of PV for diagnostics is that dynamical and physical forcing [the two terms on the rhs of Eq. (7)] are naturally separated. Admittedly, the physical forcing

*F*

_{p}is dominated by diabatic forcing and, more precisely, by the vertical gradient of diabatic heating (Hoskins et al. 1985; Hoerling 1992).

Since the divergent part of the wind field is not resolved by the nonlinear gradient-wind approach, the model used in the present study is based on a variant of Eq. (7):

where now

### b. Time averaging

The model is designed to simulate only the low-frequency variability. Hence, Eq. (8) is submitted to a 5-day time average (pentad). The timescale chosen here is somewhat arbitrary but is based on the fact that it filters out most of the baroclinic activity. Our results would also hold using 10-day averages. Denoting by an overbar the 5-day average operator, any quantity *Q*(*t*), dependent on the flow, is therefore decomposed over a pentad as

The evolution of PV pentads is therefore given by

The low-frequency dynamical PV tendency is therefore due to nonlinear low-frequency interactions, −*J*(Ψ, *P*), mostly due to the large-scale advection of PV, to the transient feedback, −*J*(Ψ′, *P*′), and to the low-frequency component of the residual forcing *S*_{P}.

A time series of each term of Eq. (10) is calculated from the dataset, for each available pentad, and resampled every 2.5 days. Thus, successive pentads of the resulting series overlap by 2.5 days. The PV pentad tendency at time *t* [lhs of Eq. (11)] is estimated using finite-centered differences between the PV pentad at time *t* + 12 h and that at time *t* − 12 h. Missing dates are not taken into account in the averaging or differentiation process. The two Jacobian terms are estimated using the classic collocation grid method, spatial derivatives being calculated in spherical harmonics (truncation T31). Finally, the PV forcing *S*_{P} is calculated as a residual from the other terms. This term is expected to be noisy, especially because it contains the residual of the approximation of time derivative by finite differences. It is also to be emphasized that truncation error, which can be important for potential vorticity near the tropopause, where steep gradients occur, is contained within this term.

### c. Annual cycle removal

The model is formulated in terms of anomalies with respect to an annual cycle. Let us denote by (∼) the annual cycle averaging operator and by (*) the associated time-dependent anomalies. Any quantity Q can be decomposed as

where *τ* is the calendar time. The annual cycle of each term of Eq. (10) is removed, including that of the time derivative. The resulting anomaly equation is expressed as

In practice, the annual cycle is calculated, at a given calendar time *τ*, by averaging all available data occurring at calendar times within ±15 days. This annual cycle removal procedure has the advantage of filtering out, in principle, all harmonics of the annual cycle also. In the following, we focus on the winter season only; winter being defined as the extended season November–March, except for the study of oscillations (section 5), where data throughout the year are used.

### d. Low-dimensional phase-space projection

Equation (12), which describes the evolution of the PV pentad anomaly, is now projected onto a reduced set of modes as described in the introduction. These modes are the leading *m* EOFs of the North Atlantic PV anomaly field itself. EOFs are calculated using data throughout the year, motivated by the subsequent application of multichannel singular spectrum analysis (see section 5). The leading EOFs are actually insensitive to the presence of summer months. The projection operator is denoted by {·}. This procedure results in *m* equations governing the projected PV anomaly:

Our choice of *m* is discussed in section 3f. The *m*-dimensional vectors *T, LI, TF,* and *RF* represent, respectively, the projected components of the low-frequency PV tendency, the low-frequency interactions, the transient feedback, and the low-frequency residual forcing.

The phase-space projection is associated with spatial filtering since the leading EOFs are primarily of large scale. The first EOF (not shown) corresponds to the North Atlantic oscillation (NAO; Wallace and Gutzler 1981). Table 1 shows the percentage of variance explained by the leading *m* EOFs as a function of *m.* The leading 10 EOFs explain 71% of the variance, which is a smaller fraction than what one would obtain with geopotential heights over the same Atlantic sector [about 95% (Vautard 1990)]. This indicates that PV maps exhibit larger spatial variability than geopotential maps.

Following the notations of the introduction, the *X* vector (resolved variables) is the vector containing the leading *m* PCs, and the *Y* vector is the ensemble of all other variables on which the the PV evolution depends. This large ensemble (probably infinite) does not consist only of PV variables, but also of other variables such as wind, temperature, or even sea and land surface temperatures. Each tendency term in Eq. (14) depends on *X* implicitly, but this dependence is not exclusive and cannot be formulated analytically. For instance, transient feedback (*TF*) depends on *X* by the simple fact that, given a large-scale PV environment, transients generally follow the largest PV gradient areas, the storm tracks (Dole 1986; Lau 1988; Nakamura and Wallace 1990).

### e. Model closure

Equation (14) has the same form as Eq. (2) in the introduction, with *G*(*X, Y*) being replaced by the sum of low-frequency interactions (*LI*), *TF,* and residual forcing (*RF*). The deterministic part of each of these three forcing terms is estimated as a function of *X*, as described by Eq. (4). We approximate the conditional average 〈*Q*〉(*X*) [see Eq. (4)] of any quantity *Q* as

where the weighting function is given by

For the calculation of deterministic tendency terms, *Q* in Eq. (15) is successively replaced by *LI, TF,* and *RF.* In Eq. (14), the *X*_{i}s are the available values of the state vector, that is, the projected pentad PV anomalies, and the *Q*_{i}s are the simultaneous values of *Q.* The pool of possible analogs *X*_{i} is also restricted to the winter season. By definition of the weighting function, only analogs of *X* are taken into account in the averaging process. Note also that 〈*Q*〉(*X*) is a differentiable function of *X.*

In Eq. (16), *D* is a reference radius estimated, as in Michelangeli et al. (1995), as a given percentile of the distribution of distances between observed state vectors. Michelangeli et al. (1995) took the value below which 10% of this distribution falls. Here we use the value of *D* obtained for 15% in order to keep approximately the same number of analogs in the neighborhood of *X,* since our PV record is about 1.5 times shorter than the geopotential dataset used therein. In any event, the results presented below depend weakly on *D*: taking 10%, 20%, or even 30% yields the same qualitative results.

Figure 1 shows a scatter diagram of the pattern correlation, as a function of the distance, between pairs of PV state vectors projected onto the leading six EOFs. This cloud of points is obtained from about 250 data vectors resampled every 15 days in order to avoid time correlation. The distance *D* achieving 15% of the distance distribution is represented by the vertical bar in the figure. Most of the pairs less than a distance *D* apart have relatively large pattern correlations. The typical pattern correlation associated with the worst analogs taken into account in the weighted average (15) is about 0.5. The negatively correlated pairs are mostly low-amplitude pairs, for which pattern correlation is not a suitable measure of distance.

The removal of the nondeterministic part of the tendency leads to the final formulation of the low-order model, given schematically by

### f. Choice of the phase-space dimension m

Throughout this article (unless specified otherwise), we use *m* = 6, but all our conclusions remain valid for *m* values varying between 5 and 10. In order to justify more clearly the choice of this range, we compute the average pattern correlations between the terms of Eq. (14), *T, LI, TF,* and *RF* and their deterministic component 〈*T*〉(*X*_{k}); 〈*LI*〉(*X*_{k}), 〈*TF*〉(*X*_{k}), and 〈*RF*〉(*X*_{k}) are calculated for the data state vectors *X*_{k}. When estimating deterministic components using Eq. (15), we remove from the pool of possible analogs *X*_{i} those occurring within ±15 days about the date of *X*_{k}, in order not to artificially overestimate the average pattern correlation.

These pattern correlations are shown as a function of *m* in Fig. 2. All curves increase with *m,* obtain a maximum, and then decrease. The increase for low values of *m* is due to the high information content of the very first PCs, while the final decrease is due to the poorer quality of analogs for large values of *m.* The optimal value of *m* for the parameterization of the transient feedback is *m* = 6; it is even lower for the residual forcing (*m* = 4), but higher for the low-frequency interactions (*m* = 9). The most deterministic term is *LI*, with pattern correlations of about 0.45. Next, the residual forcing has a maximal pattern correlation of about 0.30, and the transient feedback is poorly parameterized with pattern correlations of 0.2 at best. The PV tendency itself has a maximal pattern correlation of about 0.3, occurring strikingly at the larger value *m* = 14, but the corresponding curve undergoes fluctuations beyond *m* = 8 that are probably statistically nonsignificant. The model should perform best for values of *m* between 5 and 10.

## 4. Quasi-stationary weather regimes

### a. Weather regime patterns

As presented in the introduction, weather regimes *X** are the equilibria—also called the “fixed points”—of the dynamical system (17), satisfying thus

Note that this definition differs from the classic definition of multiple equilibria as proposed by Charney and Devore (1979). In this work, the low-order model describes the evolution of the largest scales but does not include transient feedback or residual forcing parameterization. Roughly speaking, Charney and Devore–type multiple equilibria would be found by seeking vectors *X** for which only *LI* vanishes in Eq. (18).

Solutions of Eq. (18) are obtained using a quasi-Newton descent algorithm (Nocédal 1980). Since Eq. (18) results from a statistical calculation, we do not expect to find exactly the same weather regimes as if another independent sample of “observed data” were used. Vautard (1990) devised a parametric index of statistical significance of the solutions, based on the distribution of the norm of the lhs of Eq. (18) that would obtain when changing the sample. Our purpose here is not to document all solutions obtained. Instead, we focus on the two most significant solutions, that is, the ones having the largest significance index. These solutions are also qualitatively recovered when using only the first or second half of the dataset. For each value of the dimension *m* between 5 and 10, two major groups of solutions are found, which are similar to those obtained using *m* = 6.

The anomaly pattern of these two solutions is represented in Fig. 3. Both exhibit a dipole anomaly structure. The first one is clearly related to the Euro–Atlantic blocking phenomenon, with a typical reversed potential vorticity gradient over western Europe. It will be called hereafter the “blocking regime” (BL). The other solution, having almost exactly the opposite anomaly, will be called the “zonal regime” (ZO). These two solutions are also among the significant solutions found using geopotential heights as the state vector (Vautard 1990; Michelangeli et al. 1995). Other solutions are found, similar also to the patterns found in these previous studies, but bear a much smaller statistical significance.

### b. Maintenance of the weather regimes

The question of the maintenance of the weather regimes, mainly the blocking regime, has been addressed thoroughly in previous studies. Charney and Devore (1979), and later Legras and Ghil (1985), explained the maintenance of high-amplitude eddies by the wave-mean flow resonant interactions in the presence of topography, a description that is somewhat inconsistent with the surface boundary conditions in the Euro–Atlantic sector. Transient feedback, which was absent in their models, was suggested to be a key factor by Shutts (1983), Hoskins et al. (1983), and Pierrehumbert (1985) in the maintenance of modon-type stationary structures. Many observational studies (Illari and Marshall 1983; Mullen 1987; among others) and modeling studies (Reinhold and Pierrehumbert 1982; Vautard and Legras 1988; Molteni 1995) supported this standpoint. The still debatable question is whether transients are really necessary to maintain multiple weather regimes such as blocking, or whether they only reshape or reinforce large-scale equilibria of a model that would not incorporate unresolved-scale parameterization. Arguments in favor of the second hypothesis are that Euro–Atlantic blocking strongly resembles the classic modon structure (McWilliams 1980; Haines and Marshall 1987), which is stationary when embedded in a zonal westerly flow compensating for the westward self-advection of the dipole. We present below some observational evidence for this hypothesis.

Figure 4 shows the deterministic forcing anomalies associated with the three terms leading to the balance of Eq. (18), for BL and ZO. These terms sum up to zero by definition. For both BL and ZO, *TF* is constructive, but mostly propagative in the upstream direction, while *LI*s tend to advect the anomalies downstream. The residual forcing is of the same order of magnitude as the other terms and mainly dissipates the anomaly.

These results are fairly consistent with those of Mullen (1986, 1987) who studied the composite transient forcing for high-amplitude anomalies, but without assuming quasi-stationarity. The net dissipative effect of residual forcing can be understood by physical considerations on diabatic heating: anticyclonic anomalies are usually associated with dry weather during which the earth’s radiative cooling leads to a positive vertical gradient of diabatic heating, resulting in a net positive PV forcing. On the other hand, cyclonic vorticity anomalies are associated with middle-atmosphere excess of diabatic heating, and therefore a negative gradient of diabatic heating near the tropopause, resulting in a net positive PV forcing. These arguments are, however, somewhat speculative since the residual forcing also contains an important contribution from the divergent forcing.

In order to quantify relationships between anomaly and deterministic forcing patterns, we show in Table 2 the pattern correlations between regime anomalies and the forcing terms, and we compare these correlations with those obtained on average over the whole set of observed state vectors. As for the construction of Fig. 2, nearly simultaneous analogs are removed from the pool of possible analogs in the calculation of the deterministic part of forcing terms taken at data points. On average, low-frequency interactions are slightly constructive but almost orthogonal to the weather regimes’ anomalies. On the other hand, transients are dissipative on average but constructive for weather regimes. Therefore, any parameterization of the Reynolds stresses as a damping would be in error for the particular case of the weather regimes. In all cases, residual forcing is dissipative.

### c. The large-scale balance of weather regimes

In order to investigate further the contribution of each term in the observed regime patterns, we calculate the stationary solutions that one would obtain when one or more of the three forcing terms are omitted in Eq. (18). First of all, when both transient forcing and residual forcing are omitted, two statistically significant solutions are obtained (Figs. 5a and 5b), which are distorted versions of the ZO and BL anomalies (cf. with Fig. 3). The existence of such solutions indicates that there are phase-space areas where the low-frequency interactions are vanishing or weak and that transient and residual forcings actually displace quasi-stationary solutions from these equilibria, in agreement with the standpoint of Reinhold and Pierrehumbert (1982), and the modeling analysis of Molteni (1996). This result also supports the expansion about the stationary solutions of Pierrehumbert and Malguzzi (1984), who analyzed the departures of stationary solutions in the presence of weak transient or boundary forcing.

If one omits only the transient feedback or the residual forcing, zonal and blocking solutions are still found (Figs. 5c–f), also with some distortions in the patterns themselves, confirming the previous hypothesis. By contrast, when the low-frequency interaction term is removed, the balance is achieved for solutions having very different patterns than the above blocking or zonal solutions (not shown), indicating once again that low-frequency interactions are a key factor for the balance to be achieved.

Another way to confirm the large-scale balance is to compare the amplitude (the rms of the anomaly) of each term of the balance for the quasi-stationary flows with the amplitudes obtained for the data points. More precisely, at each observed PV state vector *X*_{j}, as well as for BL and ZO, we compute the norm of the deterministic part of *LI*, *TF*, and *RF*, and report it as a function of their number of analogs in Fig. 6. In fact, the horizontal axis of this graph is not exactly the number of analogs of *X*_{j}, but the *equivalent number of degrees of freedom N*_{e}(*X*_{j}) (Vautard 1990) used to compute the conditional average [Eq. (14)]. For any state vector *X, N*_{e}(*X*) is given by

The distribution of the amplitude decreases as the number of degrees of freedom increases, not surprisingly, since large statistical fluctuations dominate the amplitude when a small number of data points lie in the neighborhood of the data point *X*_{j}. Another reason could be that points with few neighbors are actually transient states where phase-space velocity is high, hence the large amplitudes. The values of *N*_{e}(*X*) for the weather regimes are high, a necessary condition for them to be statistically significant. However, Fig. 6a shows clearly that the amplitude of the anomalous *LI* forcing is much weaker than “normal” for a state vector having that many degrees of freedom. Such is not the case for the other forcing terms (Figs. 6b and 6c). Therefore, the equilibration of the two weather regimes is achieved, to first order, by the balance of the *LI* term.

The above equilibration experiments and budgets are only *suggestive* of causal relationships; should one be able to remove physically one of the above three forcings, the atmospheric system would find, in the long term, a different climatic balance. In particular, removal of a physical process could modify qualitatively other physical processes, since forcing terms are themselves significantly correlated (see Table 3). In particular, the primary large-scale balance found here could be deeply modified or destroyed by removal of transients.

## 5. Intraseasonal oscillations

Fixed points are not the only generic features of dynamical systems. Periodic orbits are associated with oscillatory behavior. Depending on their stability, they can lead to intermittent spells of oscillations (Plaut and Vautard 1994). The objectives of the remainder of this article are

to identify intraseasonal oscillations in the PV anomaly fields using multichannel singular spectrum analysis (this section),

to calculate forcing budgets as in the previous section (this section), and

to show that periodic orbits are relevant prototypes for these oscillations (section 6).

### a. Multichannel singular spectrum analysis of the PV anomaly fields

Multichannel singular spectrum analysis (MSSA) has been used in several recent studies (Kimoto et al. 1991; Plaut and Vautard 1994; Lau et al. 1994). It is a data analysis technique based on a lag-space principal component analysis. As such, it does not differ from the classic extended EOF analysis of Weare and Nasstrom (1982), but the use of long windows, which was not the original design of extended EOF analysis, allows the identification of oscillations (with periods shorter than the lag window) and the derivation of important spectral properties (Plaut and Vautard 1994). The major property of MSSA is that it is able to identify highly intermittent oscillations that would not necessarily correspond to a peak in the power spectrum of the multichannel time series. These oscillations show up as pairs of eigencomponents in phase quadrature with each other. For technical details about the MSSA technique, the reader is referred to Plaut and Vautard (1994).

We apply MSSA to the 5-day average PV anomaly field over the North Atlantic sector. We follow the procedure of Plaut and Vautard (1994) by taking as “channels” the six leading PCs of the PV pentads calculated in section 3. We use a lag window of *W* = 200 days, which, in principle, allows one to study oscillations with periods in the range of 20–100 days, roughly speaking. Data are sampled every 2.5 days, and therefore the size of the matrix to diagonalize is 6 × 80 = 480. Data throughout the calendar year are used because the length of the window does not allow one to use winter data only. The eigenfunctions of the lag-covariance matrix are called the space–time EOFs (ST-EOFs), and the corresponding time-dependent projection coefficients are called the space–time principal components (ST-PCs).

We select oscillatory pairs of ST-PCs, using the criterion developed by Plaut and Vautard (1994). This criterion requires that the correlation function between ST-PC *k* and ST-PC *k* + 1 has at least one maximum above 0.5 and one minimum below −0.5 on each side of the lag axis. Under the above conditions, three pairs are selected from our analysis, pairs 7–8, 20–21, and 23–24. The associated periods are estimated from a maximum entropy (ME) spectral analysis of the ST-PCs. ME spectra peak at 70 days for pair 7–8, 33 days for pair 20–21, and 28 days for pair 23–24. We did not analyze ST-PCs of a rank larger than 30. Note the good correspondence between the periods of the oscillations found in Plaut and Vautard (1994), who did not analyze periods smaller than 30 days. The robustness of these oscillations was confirmed by performing many experiments with different values of the method’s parameters. All results presented below remain valid as long as the number of channels is larger than 3 and when the window length is extended from 150 days to 1 yr.

The space–time behavior of these oscillations is studied by a phase-composite analysis of the *reconstructed oscillations,* similar to that of Plaut and Vautard (1994). Phase and amplitude indices *θ*(*t*) and *A*(*t*), defining the “instantaneous state” of the oscillation, are calculated in three steps.

- We calculate the reconstructed oscillation,
*R*(*t*), from ST-PCs (*k,**k*+ 1). The reconstruction process consists essentially of finding the linear combination of the ST-PC pair that best fits the original signal. The oscillation*R*(*t*) is a vector of the same dimension as the “state vector,” that is, six in our case. It is noteworthy that the sequence of operations to obtain the*m*-dimensional time series*R*from the original series of PV anomalies*P** is a sequence of linear filters, which is denoted in the following by Γ:(20)*R*= Γ(*P**). The time derivative

*d*[*R*_{n}(*t*)]/*dt*of the reconstructed oscillation in the channel*n,*which has the largest reconstructed variance, is calculated by finite-centered difference. Both the reconstructed channel*R*_{n}(*t*) and its time derivative are normalized by their standard deviations.The modulus

*A*(*t*) and argument*θ*(*t*) of the complex number*Z*(*t*) =*d*[*R*_{n}(*t*)]/*dt*+*iR*_{n}(*t*) define, respectively, the instantaneous amplitude and phase of the oscillation.

The mean state of the oscillation at phase *θ* can be calculated by averaging the reconstructed oscillation *R*(*t*) over all occurences *t* when |*θ*(*t*) − *θ*| ≤ *ε*, where *ε* is a suitably chosen constant angular parameter. We use *ε* = 22.5°. The resulting composite values *R*(*θ*) are called the *phase composites.* It is now a periodic six-dimensional time series that represents the “average life cycle” of the oscillation. By multiplying the six coefficients of *R*(*θ*) by the corresponding leading EOFs, one obtains the average cycle in physical space.

### b. Space–time patterns of the intraseasonal oscillations

In Figure 7, we display half the life cycle of the 70-day oscillation by plotting the contours of the reconstructed PV anomalies over phases differing by an angle of 45°, that is, four maps for one half-cycle (phases 1–4). As found in Plaut and Vautard (1994), it consists essentially of the poleward propagation of a dipole pattern over the western part of the Atlantic domain. In phase 1, the oscillation projects mainly onto the first EOF, the NAO. During the other half of the cycle (phases 5–8, not shown), the anomalies of the reconstructed PV field are almost perfectly symmetrical. Owing to filtering effects, the amplitude of the reconstructed oscillation is weak compared with standard anomalies of the PV field.

Figure 8 shows the half-cycle of the 33-day oscillation, also identified by Plaut and Vautard (1994). The oscillation has a westward propagating space–time pattern, with a phase resembling the Euro–Atlantic pattern (phase 4). Figure 9 shows the last of the three oscillations, with a 28-day period. The oscillation is a wave train propagating eastward across the North Atlantic Ocean.

### c. Maintenance of the intraseasonal oscillations

In order to calculate forcing budgets over the life cycle of the intraseasonal oscillations, the reconstruction filter Γ is applied to Eq. (14). Therefore, the time derivative of the reconstructed oscillation writes

Each term is then averaged over constant phases *θ* of the oscillation, yielding the average forcing budget of the oscillation,

where, for the sake of simplicity, we denote by *LI*(*θ*) [respectively, *TF*(*θ*) and *RF*(*θ*)] the reconstructed phase composite of the forcing terms. Figure 10 shows the average cycle of the reconstructed PV anomalies *P*(*θ*) projected onto the plane spanned by the first two EOFs (closed curve), as well as the associated reconstructed phase composites of the three forcing terms, for each oscillation (the arrows). The sum of the three arrows is the reconstructed tendency *T*(*θ*) (not shown) and is, by definition, tangent to the average cycle.

Low-frequency interactions clearly play a dominant role in the dynamics of the 70-day oscillation. Both *TF* and *RF*, which are small, tend to dissipate the oscillation anomaly. Transients therefore have, in the maintenance of the dominant intraseasonal oscillation, a role that is opposite to that found for the maintenance of the weather regimes. The low-frequency interactions tend, on the contrary, to amplify, but mainly propagate the oscillation.

For the 33-day oscillation, all terms are of equal importance. Low-frequency interactions mainly amplify the anomaly, whereas both *TF* and *RF* have a dominant propagative component. Transient feedback is almost conservative, that is, tangent to the average cycle, and accounts for about 1/3 of the total tangent velocity. Since the oscillation is westward propagating, TF accelerates the westward propagation, as for the blocking structure. One easily understands that this enhanced acceleration effect may be very important for nearly resonant waves, since it inhibits in situ growth of the associated anomalies. A general circulation model having a bad representation of transient feedback, for example, underestimation, would produce large errors in the dynamics of nearly resonant phenomena. For this oscillation, the residual forcing has a major contribution, mostly propagative and somewhat dissipative.

For the eastward propagating 28-day oscillation, all terms are nearly tangent to the phase-space average cycle. The oscillation is mostly driven by the *LI,* whereas both *TF* and *RF* decelerate the eastward propagation but are one order of magnitude smaller than *LI.* Once again, transients act as an upstream propagator, which seems to be a quite general behavior.

### d. Phase locking

MSSA usually does not gather into a single pair an oscillation and its harmonics, but creates an oscillatory pair for each significant harmonic. Plaut and Vautard (1994), and Lau et al. (1994), showed that the 70- and 33-day oscillations are frequently phase locked. By calculating joint histograms of phase occurrences of the two oscillations, as in Plaut and Vautard (1994), we also find the same result here (not shown). Phase locking can be sketched as in Lau et al. (1994), by plotting, for each of the complete 25 winters, ST-PC 20 against ST-PC 7 (Fig. 11). For many winters, the trajectory displays a shape reminiscent of the Lissajoux curve. Note, however, that phase locking is not systematic and can only be demonstrated from statistics.

## 6. Time-dependent solutions of the low-order model

### a. Conservation principles and the corrector

The aim of this section is to show that the low-order model formulated by Eq. (17) provides time-dependent solutions having qualitatively realistic dynamics. By definition, the low-order model has the same regimes as those identified in section 4. Such is not necessarily the case for intraseasonal oscillations, since nondeterministic components of the tendencies may have a major and systematic effect on their maintenance. Therefore, we focus here on the recovery of intraseasonal oscillations by model simulations. The major problem that occurs when integrating the deterministic system (17) is the violation of energy and enstrophy conservation principles. This problem is inherent to our formulation since the geographical domain covered by the model is only the North Atlantic domain. An undesirable consequence is that time-dependent variables may blow up within finite time.

This problem is enhanced by the fact that time-dependent solutions may occasionaly spend some time in phase-space areas where observed data are missing or highly scattered. In such a case, the small number of analogs taken into account in the conditional averaging process [Eq. (15)] provides a largely erroneous estimate of the average tendency, which may point toward an area of even lower density of data state vectors.

To ensure that the model’s solutions avoid low-density areas, we add to the deterministic equation a correction term, which is active only in those regions, and points toward higher-density areas. The density of neighboring points can be estimated as the equivalent number of degrees of freedom *N*_{e}(*X*) defined by Eq. (19), which is a differentiable function of the current state vector *X.* The corrector term is thus calculated as

where

In Eq. (23a), *λ* is a suitably chosen constant. In Eqs. (23b,c) the maximal number of degrees of freedom for the corrector to be active, *N*_{o}, is chosen as small as possible, at least smaller than the values of the number of degrees of freedom of the weather regimes identified in section 4. The corrector points in the direction of the gradient of *N*_{e}(*X*), that is, toward higher-density areas. However, if *N*_{o} is too small, the estimate of this gradient can be totally wrong, due to undersampling, and the corrector points toward a random direction. By trial and error, we found that values of *N*_{o} between 25 and 100 provided equivalent and satisfactory results. Unless specified, the following results are obtained with *N*_{o} = 50. The value of *λ*, to which the model is fairly insensitive, was also chosen by trial and error, *λ* = 0.5 PVU^{2} day^{−1}.

Note that fluctuations of *N _{e}*(

*X*) can, in principle, arise from two sources. The first one is the distribution of analogs within the ball; if many poor or high quality analogs have a fairly uniform distance to the center, then

*N*

_{e}(

*X*) is equal to the number of analogs

*N,*while it is much smaller than

*N*when distance distribution is very inhomogeneous. The second source is the fluctuation of the number of analogs itself. In practice,

*N*

_{e}(

*X*) is almost uniformly found to be about half of the number of analogs

*N*; hence, its small values occur exclusively when the phase-space trajectory visits “desertic” areas.

### b. Time-dependent solutions

The model (17) is integrated with a predictor–corrector scheme, using a time step of 1/2 day, starting from the vicinity of the zonal solution. In practice, possible observed analogs are cataloged in a fixed array, as well as the associated instantaneous observed tendencies. Then at each time step the analogs falling within a ball of radius *D* are selected, and the associated tendencies are averaged to estimate the deterministic tendency and the corrector. The model is first run for 2500 days, in order to remove transient behavior. Starting from the final condition of this preliminary run, the model is integrated for 25 000 days, and we display in Fig. 12a the time behavior of the first principal component (the first variable of the state vector) during the first 250 days. In order to show the sensitivity to the corrector’s parameters values, we also display the time evolution of the solution obtained with (*N*_{o} = 25, *λ* = 0.5) and (*N*_{o} = 50, *λ* = 1). In Fig. 12b, the time behavior of the ratio *K*(*t*) of the norm of the corrector by the norm of the uncorrected tendency is displayed, for the above three sets of parameter values. Clearly, the time evolution is fairly insensitive to corrector’s parameters values, since the three trajectories in Fig. 12a remain close to each other up to a lead time of about 100 days, although in the beginning of this time sequence there is a long period of active corrector. The relative amplitude of the corrector rarely exceeds 50% of that of the uncorrected tendency.

Figure 12c shows a longer time sequence of the evolution of PC 1, and Fig. 12d shows the corresponding corrector relative amplitude ratio. The solution, which is chaotic, has an intermittent oscillatory behavior, with fluctuations of short and longer periods. These oscillations cannot be produced by the artificial corrector, since they also show up during long periods when it is not active.

### c. Intraseasonal oscillations simulated by the low-order model

In order to check whether the low-order chaotic model reproduces the intraseasonal oscillations of section 5, we apply MSSA to the long 25 000-day simulation, using the same MSSA parameters as in section 5. Two oscillations, with respective periods of about 200 days and 50 days are identified in the low-order model run. Figure 13 shows the phase-composite analysis of the two oscillations, as in Fig. 10. The structure of the simulated 200-day oscillation, as well as the structure of its associated forcing, is qualitatively similar to those of the observed 70-day oscillation. Note, however, that the amplitude of the simulated oscillation is much larger than that of the observed one. Despite the erroneous period and amplitude, the low-order model reproduces qualitatively the 70-day oscillation, since anomaly patterns and forcing budgets are consistent with the observed ones. The same conclusions hold for the similarity between the simulated 50-day oscillation and the observed 28-day oscillation. However, the observed 33-day oscillation is not reproduced by the low-order model.

The fact that forcing budgets of the two simulated oscillations are qualitatively correct is suggestive that these phenomena are of internal origin. Indeed, the low-order model tendencies are expressed only as functions of the extratropical Atlantic flow. Should these oscillations be, for instance, a direct response of some tropical oscillations, such as the Madden and Julian (1971) oscillation, the deterministic part of the tendencies would not contain the oscillatory signal.

The differences between simulated and observed oscillations may arise from two main factors.

The deterministic-composite operator picks out many poor analogs, leading to a smoothing of the actual tendency vector field and removal of important dynamical features.

The tendencies actually depend on other important independent factors, such as surface fluxes, tropical flow, or unresolved scales. That is, our “deterministic” formulation misses important dynamical effects that cannot really be parameterized only as a function of the North Atlantic PV anomalies only and are therefore rejected in the nondeterministic forcing.

In order to address these possibilities, we plot in Fig. 14 the forcing budgets of the *deterministic* part of tendency terms along the average life cycle of the three *observed* average oscillations, again in the same format as in Fig. 10. For the 70- and 28-day oscillations, the deterministic forcing composites are fairly similar to the actual forcing composites of Fig. 10, with slight direction displacements but large-amplitude underestimation of the *LI* and *TF* forcing terms. For these oscillations, deterministic tendency composites are systematically underestimated by a factor of 2–3, the same factor as the simulated–observed period ratio. By contrast, deterministic forcing budgets for the 33-day oscillation are totally erroneous, especially for the low-frequency interaction term.

Figure 14 does not allow one to draw definite conclusions about the error sources. Nevertheless, the *LI* term of the 33-day oscillation is clearly badly simulated, with systematic direction bias. Such a bias cannot result from averaging effects (1) exclusively. Indeed, averaging over poor analogs probably leads to an artificial reduction of the amplitudes of the composites, but certainly not to a systematical bias of their phase. We are left with the conclusion that nondeterministic forcing is a key component of the dynamics of the 33-day oscillation. This nondeterministic forcing could arise, for instance, from the interactions with the tropical atmosphere and the so-called Madden–Julian oscillation.

### d. Periodic orbits and intermittency

When the dimension of the model is decreased to 3, the model becomes nonchaotic, and time-dependent trajectories converge to either a stable periodic orbit with a period of 75 days or another periodic orbit with a period of 225 days. These periodic orbits are displayed in Fig. 15. Along the 75-day orbit, the corrector is not active, whereas along the 225-day orbit it is active in some places, but variations in its parameters do not change crucially the shape of this orbit, as shown by the orbit obtained with *N _{o}* = 25 and

*λ*= 0.5. Budgets of forcings along these orbits (not shown) display the same character as the observed ones, meaning that, again, despite the amplitude and period differences, they correspond qualitatively to the 70- and 28-day observed oscillations.

When the model is integrated in dimension 4, the dynamical system becomes chaotic. Figure 16 shows a time sequence of the evolution of PC 1. This sequence does not correspond to a spinup sequence, but has been obtained after a long time of integration. The intermittency phenomenon is quite clear from this sequence. The trajectory spirals along a weakly unstable 60-day periodic orbit (sequence 1), before being ejected off its neighborhood and undergoing longer-period fluctuations around a second periodic orbit (sequence 2). A careful investigation of these two sequences show that they indeed correspond to time periods when the phase-space trajectory lies, respectively, close to two unstable orbits qualitatively similar to that displayed in Fig. 15. As the dimension of the low-order model is increased, clean intermittent sequences, as in Fig. 16, become hard to identify, due to the larger instability of underlying periodic orbits.

The above results suggest that the 70- and 28-day observed intraseasonal oscillations are due to the existence of periodic orbits in the atmospheric dynamical system, resulting from the equilibration between large-scale interactions, transient feedback, and physical processes. However, these orbits must be unstable in the real atmosphere, which makes their identification impossible in practice.

## 7. Summary and conclusions

This article proposes a general methodology to build low-order deterministic, albeit empirical, models given long observational series of the system under study. These models are based on the projection of the evolution equations onto a small number of leading empirical orthogonal functions. The model is made deterministic by an empirical closure of the equations including the feedback effects of unresolved processes. The model formulation is applied to the diagnostics of the extratropical variability, with focus on the Euro–Atlantic domain. In this domain, the generation of low-frequency variability by alternation between several weather regimes cannot be due solely to topographic instability. Reinhold and Pierrehumbert (1982) and Vautard and Legras (1988) proposed instead an internal mechanism of equilibration between large-scale interaction and transient feedback (and other physical processes), but these were verified only on model simulations. Here, a long series of potential vorticity on an isentropic surface (315 K) is analyzed. We also attempt to extend the equilibration mechanism to the maintenance of intermittent intraseasonal oscillations. In the parameterized large-scale model formulation, equilibrated weather regimes are found as unstable stationary solutions while oscillations are found as periodic orbits.

The main diagnostic results of the weather regime study are as follows.

The two most statistically significant weather regimes correspond to symmetrical blocked and zonal dipolar anomalies located over western Europe. These regimes were also identified by Vautard (1990) and Michelangeli et al. (1995) using the same methodology but geopotential height data, making diagnostics rather difficult.

The budget of forces maintaining the weather regimes shows that all three terms have comparable magnitudes. Large-scale interactions mostly advect the two regime anomalies downstream, whereas transient feedback advects anomalies upstream, a result consistent with the findings of Illari and Marshall (1983) and Mullen (1987). Transients also have an amplification effect. This is in contrast with the classic eddy diffusivity principles, where actually unresolved turbulence has dissipative effects. On average, however, transient feedback indeed acts as an enstrophy dissipator, as expected from the two-dimensional turbulence theory. Finally, residual forcing is mostly dissipative, which is consistent with our knowledge of the vertical structure of diabatic heating.

Weather regimes lie in phase-space areas where the large-scale interactions are already weak, achieving a “zero-order balance.” The principal mechanisms of formation of multiple weather regimes seem therefore already present in the large-scale balance, but regime patterns are reshaped by the transients and residual processes.

Intraseasonal oscillations of periods of 70, 33, and 28 days are identified using multichannel singular spectrum analysis (MSSA). The space–time behavior of the first two Euro–Atlantic oscillations is similar to that obtained by Plaut and Vautard (1994). New results about these phenomena are the following.

All oscillations result from a combination of all three forcing terms. However, while low-frequency interactions play a dominant role in the maintenance of the 70- and 28-day oscillations, the three forcing terms are equally important for the 33-day oscillation.

Transient feedback has a dissipative effect in the maintenance of the 70-day oscillation, an upstream propagative effect for the 33-day oscillation, and an upstream counterpropagative, albeit very small, effect for the 28-day oscillation. Transients, as for the blocking regime, have a general tendency to advect the anomalies upstream.

Finally, time-dependent solutions of the low-order model are calculated, and two of the three oscillations are qualitatively reproduced. The periods of the simulated oscillations, however, are too long and their amplitudes too large. In spite of that, their space–time behavior and forcing budgets are in agreement with those of observed oscillations.

The discrepancy between observed and simulated periods has two likely sources. One is a technical problem: the available amount of data enforces the empirical parameterization of tendency terms, based on an analog procedure, to be too smooth and to take large neighborhoods, reducing the amplitude of parameterized tendencies. The other source is the lack of information given in the low-dimensional state vector for the parameterization. We restricted this state vector to PV anomalies over the North Atlantic, but additional dynamical and physical forcing can also originate from outside the sector and/or from the vertical structure of the flow, which is not taken into account here.

When the dimension of the model’s state vector is decreased to three, the two simulated oscillations appear as two stable periodic orbits. In dimension four, these two orbits become unstable, and time-dependent trajectories, which are chaotic, undergo transitions from one orbit to the other, leading to a well-marked intermittent behavior. In higher dimensions, it is likely that these two periodic orbits still exist, but the model shows a very chaotic behavior and the oscillations can only be detected from statistical methods.

Above all, the methodology used here is mainly empirical. All terms of the model are estimated from the analog procedure. The degree of realism of the present modeling approach would certainly gain from the inclusion of analytical expressions for the various forcing sources, in particular for the advection. By extending the domain both horizontally and vertically, one would probably increase the deterministic part of forcing terms, that is, include more information in the flow structure and consequently in the flow-dependent part of the forcings.

## Acknowledgments

We are particularly thankful to Christopher Strong for a careful review of the manuscript. We also thank Guy Plaut for his participation in the development of the MSSA code. Two anonymous reviewers helped a lot in improving the initial manuscript. This work was supported by the French power company, Electricité de France, and the European contract “Short Term Climate Variability.” The stay of Eduardo Da Costa at LMD was granted by the Junta Nacional Investigacao Cientifica e Tecnologica.

## REFERENCES

## Footnotes

*Corresponding author address:* Dr. Robert Vautard, Laboratoire de Météorologie Dynamique, Ecole Normale Superieure, 24, rue Lhomond, 75231 Paris Cedex 05, France.

Email: vautard@lmd.ens.fr