Michael Ghil
email: ghil@atmos.ucla.edu

Andreas Groth

Dmitri Kondrashov

Andrew W Robertson

Extratropical sub-seasonal-to-seasonal oscillations and multiple regimes: The dynamical systems view

Keywords: bifurcations, cluster analysis, empirical model reduction, hierarchy of models, singular spectrum analysis, predictability

This chapter considers the sub-seasonal-to-seasonal (S2S) prediction problem as intrinsically more difficult than either short-range weather prediction or interannual-to-multidecadal climate prediction. The difficulty arises from the comparable importance of atmospheric initial states and of parameter values in determining the atmospheric evolution on the S2S time scale.

The chapter relies on the theoretical framework of dynamical systems and the practical tools this framework helps provide to low-order modeling and prediction of S2S variability. The emphasis is on mid-latitude variability and the complementarity of the nonlinear-waves vs. multiple-regime points of view in understanding this variability. Empirical model reduction and the forecast skill of the models thus produced in real-time prediction are reviewed.

Introduction and motivation

The John von Neumann point of view [START_REF] Von Neumann | Some remarks on the problem of forecasting climatic fluctuations[END_REF]): short-term numerical weather prediction (NWP) is the easiest -i.e., it is a pure initial-value problem; long-term climate prediction is next easiest -it corresponds to studying the system's asymptotic behavior; intermediateterm prediction is hardest -both initial & boundary values are important. The modeling hierarchy & successive bifurcations as Ariadne's thread through the rungs of the hierarchy.

Intraseasonal time scales, more recently called subseasonal-to-seasonal (S2S), range from the deterministic limit of atmospheric predictability, of about 10 days, up to a season, say 100 days. These time scales occupy a window of overlap between low-frequency variability (LFV) intrinsic to the atmosphere and short climatic time scales that also involve the upper ocean and land-surface features, as well as the stratosphere. These time scales are of particular importance to subseasonal prediction. Theoretical and observational studies of LFV over the past half-century have used two complementary ways of describing atmospheric LFV in the extratropics: (i) episodic, by means of multiple weather [START_REF] Reinhold | Dynamics of weather regimes: Quasi-stationary waves and blocking[END_REF] or flow [START_REF] Charney | Multiple flow equilibria in the atmosphere and blocking[END_REF][START_REF] Legras | Persistent anomalies, blocking and variations in atmospheric predictability[END_REF] regimes; and (ii) oscillatory, by means of broad-peak, slowly modulated oscillations (Ghil and Robertson, 2002, and references therein).

We pursue throughout the chapter the dynamical systems point of view first formulated by John von Neumann [START_REF] Von Neumann | Some remarks on the problem of forecasting climatic fluctuations[END_REF], according to which there are three levels of difficulty in understanding and predicting atmospheric phenomena: (a) short-term numerical weather prediction (NWP) is the easiest, since it represents a pure initial-value problem, as formulated by V. Bjerkness [START_REF] Bjerknes | Das Problem der Wettervorhersage, betrachtet vom Standpunkte der Mechanik und der Physik[END_REF] and L. F. Richardson [START_REF] Richardson | Weather Prediction by Numerical Process[END_REF]; (b) long-term climate prediction is next easiest, since it corresponds to studying the system's asymptotic behavior, i.e., the possible attractors -such as fixed points, limit cycles, strange attractors -and the statistical properties thereof [START_REF] Ghil | Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dynamo Theory, and Climate Dynamics[END_REF][START_REF] Dijkstra | Low-frequency variability of the large-scale ocean circulation: A dynamical systems approach[END_REF][START_REF] Dijkstra | Nonlinear Climate Dynamics[END_REF]; and (c) intermediate-term prediction is hardest, since both the initial and the parameter values are important. We deal with the latter problem by pursuing a full hierarchy of models, from the simplest to the most detailed ones [START_REF] Schneider | Climate modeling[END_REF][START_REF] Held | The gap between simulation and understanding in climate modeling[END_REF], and the tool of successive bifurcations as Ariadne's thread through the rungs of this hierarchy [START_REF] Ghil | Solving problems with GCMs: General circulation models and their role in the climate modeling hierarchy[END_REF][START_REF] Ghil | Hilbert problems for the geosciences in the 21st century[END_REF].

This chapter starts in Sec. 2 with a brief summary of the main characteristics of the oscillatory vs. intermittent viewpoints from a dynamical systems perspective, based on the review article of Ghil and Robertson (2002) and updated using studies published since then. Section 3 then presents a new analysis of the observational characteristics of oscillatory LFV, using a recently developed multivariate spectral analysis tool. With this background and results in hand, Sec. 4 discusses low-order modeling (LOM) approaches and their application to LFV, while outlining how these approaches have developed over the years from so-called linear inverse models (LIMs) or principal oscillation patterns (POPs) to empirical model reduction (EMR) and multi-level stochastic models (MSMs). The chapter concludes in Sec. 5 with the prospects for improving S2S prediction in the extratropics based on the theoretical LFV framework sketched herein and on data-driven LOMs.

Multiple mid-latitude regimes and low-frequency oscillations

Blocking, zonal flow, and teleconnection patterns; phenomenological description of extratropical flows in the S2S band (10-100 days). Mechanical-& thermal-topography effects, i.e., mountains & land-sea contrast. Bifurcations & symmetry breaking as the unifying point of view for the theoretical interpretation of the phenomenology. Markov chains of regimes & predictability.

The case for multiple regimes and their classification

Persistent LFV anomalies, in which the flow patterns differ significantly from the normal climatological circulation and remain stationary for more than a week, have been objectively identified over the North Pacific and North Atlantic in the 1980s [START_REF] Wallace | Teleconnections in the geopotential height field during the Northern Hemisphere winter[END_REF][START_REF] Dole | Persistent anomalies of the extratropical Northern Hemisphere wintertime circulation: Geographical distribution and regional persistence characteristics[END_REF][START_REF] Mo | Cluster analysis of multiple planetary flow regimes[END_REF]. Their onsets and breaks, on the other hand, are rather abrupt. The most familiar of these persistent patterns or "Grosswetterlagen" [START_REF] Namias | Long-range weather forecasting: History, current status and outlook[END_REF], albeit not the only ones, are blocked and zonal flows.

Within the last few decades it has been demonstrated that these patterns can be identified by examining the probability distribution function (pdf) of the corresponding large-scale maps in the atmosphere's phase space [START_REF] Cheng | Cluster analysis of the Northern Hemisphere wintertime 500-hPa height field: Spatial patterns[END_REF]Kimoto and Ghil, 1993a,b;[START_REF] Smyth | Multiple regimes in Northern Hemisphere height fields via mixture model clustering[END_REF]. The resulting patterns resemble those found previously by using correlation analysis.

Many of the methods used to classify weather maps are summarized in Table 1. In such a classification, an individual atmospheric map is thought of as a point in phase space. To achieve a reliable, statistically significant classification, it is necessary to consider a low-dimensional subspace of this phase space that still captures most of the variance in the data set. The usual choice is to compute the analyzed record's empirical orthogonal functions (EOFs), i.e., the eigenvectors of the covariance (or correlation) matrix, and to select a subspace spanned by a few leading eigenvectors [START_REF] Mo | Cluster analysis of multiple planetary flow regimes[END_REF][START_REF] Cheng | Cluster analysis of the Northern Hemisphere wintertime 500-hPa height field: Spatial patterns[END_REF]Kimoto and Ghil, 1993a,b;[START_REF] Smyth | Multiple regimes in Northern Hemisphere height fields via mixture model clustering[END_REF].

Many of the classification methods define the regimes as classes of distinct atmospheric states that have a high probability of occurrence and are separated by regions of lower probability. Some of these methods seek maxima of the pdf by using kernel density estimation (Kimoto and Ghil, 1993b;[START_REF] Corti | Signature of recent climate change in frequencies of natural atmospheric circulation regimes[END_REF] or more ad hoc methods [START_REF] Molteni | Regimes in the wintertime circulation over northern extratropics. I: Observational evidence[END_REF]. Each regime is then formed by the points, or maps, that exceed a given probability threshold in the neighborhood of a pdf maximum. The number of pdf peaks depends on the kernel smoothing parameter used, which can be determined objectively by using a least-squares cross-validation procedure [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF]. [START_REF] Smyth | Multiple regimes in Northern Hemisphere height fields via mixture model clustering[END_REF] have used a mixture model that approximates the pdf by the sum of a small number of multivariate Gaussians. In this case, the regimes are "fuzzy" in the sense that they overlap, and that each particular daily weather map can be assigned a probability of belonging to one or another regime. This was the case also in [START_REF] Mo | Cluster analysis of multiple planetary flow regimes[END_REF], who used a different classification algorithm, and obtained therewith a different number of regimes. This number, as well as the spatial patterns associated with the regimes, are discussed in Sec. 4.

Cluster analysis is a less ambitious approach to classifying atmospheric states: it localizes high concentrations of points, called clusters, but does not pretend to estimate the pdf. There are two main types of clustering algorithm: hierarchical and partitioning. In hierarchical algorithms, one builds a classification tree iteratively, starting from single data points and merging them into clusters according to a similarity criterion. [START_REF] Cheng | Cluster analysis of the Northern Hemisphere wintertime 500-hPa height field: Spatial patterns[END_REF] used Ward's method to do this. In partitioning algorithms, a prescribed number of clusters is chosen, and data points are agglomerated around kernels initially picked as random seeds. The kernels are iteratively modified so as to globally minimize the data scatter about them, as done in the K-means method of [START_REF] Michelangeli | Weather regimes: Recurrence and quasi stationarity[END_REF].

Using low-pass filtered data, based on running means of 5-10 days, introduces some measure of persistence into the above methods, which are based on frequency of occurrence. A second broad class of methods uses quasi-stationarity explicitly. Here, the regimes are defined as comprising states for which large-scale motion is slow in the statistical sense. More precisely, one seeks the large-scale patterns that have, on average, a small time derivative [START_REF] Legras | Persistent anomalies, blocking and variations in atmospheric predictability[END_REF]. This phase-space speed can be computed for maps that do include synoptic-scale motions by a nonlinear equilibration technique [START_REF] Vautard | On the source of midlatitude low-frequency variability. part II: Nonlinear equilibration of weather regimes[END_REF][START_REF] Vautard | Multiple weather regimes over the North Atlantic: Analysis of precursors and successors[END_REF][START_REF] Michelangeli | Weather regimes: Recurrence and quasi stationarity[END_REF]. [START_REF] Huth | Classifications of atmospheric circulation patterns[END_REF] provided a more recent review of classification schemes of atmospheric circulation patterns, and included a table that updates the one of Ghil and Robertson (2002) reproduced herein. [START_REF] Straus | Circulation regimes: Chaotic variability versus SST-forced predictability[END_REF] compared the circulation regimes in the National Center's for Environmental Prediction (NCEP's) reanalyses with the GCM model simulations of the Center for Ocean-Land-Atmosphere Studies (COLA), showing that three out of the four observed clusters have identifiable counterparts in the model simulations. These authors, as well as [START_REF] Straus | Circulation regimes and SST forcing: Results from large GCM ensembles[END_REF], also studied possible changes in the structure of the COLA GCM's regimes in response to changes in tropical sea surface temperature (SST) forcing. [START_REF] Christensen | Simulating weather regimes: impact of stochastic and perturbed parameter schemes in a simple atmospheric model[END_REF], [START_REF] Dawson | Simulating weather regimes: impact of model resolution and stochastic parameterization[END_REF] and [START_REF] Muñoz | A weather-type based cross-timescale diagnostic framework for coupled circulation models[END_REF] used weather regimes to evaluate more broadly the simulation skill of climate models, with [START_REF] Dawson | Simulating weather regimes: impact of model resolution and stochastic parameterization[END_REF], in particular, showing that higher horizontal resolutions in NWP models improves regime simulation.

Most recently, [START_REF] Hannachi | Low-frequency nonlinearity and regime behavior in the Northern Hemisphere extratropical atmosphere[END_REF] have reviewed the state of knowledge on the low-frequency fluctuations of the extratropical troposphere. These authors showed that retaining periods longer than 10 days yields statistically significant multiple regimes, while [START_REF] Stephenson | On the existence of multiple climate regimes[END_REF], to the contrary, had found [START_REF] Muñoz | A weather-type based cross-timescale diagnostic framework for coupled circulation models[END_REF] Hard (k-means) Model + sectorial [START_REF] Straus | Circulation regimes and SST forcing: Results from large GCM ensembles[END_REF]; [START_REF] Straus | Circulation regimes: Chaotic variability versus SST-forced predictability[END_REF][START_REF] Straus | Atmospheric regimes: The link between weather and the large-scale circulation[END_REF] Hard (k-means)

Hierarchical NH + sectorial [START_REF] Cheng | Cluster analysis of the Northern Hemisphere wintertime 500-hPa height field: Spatial patterns[END_REF] 3 NH clusters PDF estimation Univariate NH [START_REF] Benzi | The statistical properties of general atmospheric circulation: observational evidence and a minimal theory of bimodality[END_REF]; Bimodality [START_REF] Hansen | The probability density distribution of the planetary-scale atmospheric wave amplitude revisited[END_REF] Multivariate NH Kimoto and Ghil (1993a) 3 modes NH + sectorial Kimoto and Ghil (1993b); Multimodal [START_REF] Smyth | Multiple regimes in Northern Hemisphere height fields via mixture model clustering[END_REF] 3 NH clusters Regime classification by persistence Pattern correlations NH Horel (1985) SH [START_REF] Mo | Statistics and dynamics of persistent anomalies[END_REF] 3 regimes Minima of tendencies Models [START_REF] Legras | Persistent anomalies, blocking and variations in atmospheric predictability[END_REF]; [START_REF] Mukougawa | A dynamical model of "quasi-stationary" states in large-scale atmospheric motions[END_REF]; [START_REF] Vautard | On the source of midlatitude low-frequency variability. part II: Nonlinear equilibration of weather regimes[END_REF] Atl.-Eur. sector [START_REF] Vautard | Multiple weather regimes over the North Atlantic: Analysis of precursors and successors[END_REF] 4 regimes

Transition probabilities

Counts

Model + NH [START_REF] Mo | Cluster analysis of multiple planetary flow regimes[END_REF] Elementary

Monte Carlo NH + SH [START_REF] Vautard | Statistical significance test for transition matrices of atmospheric Markov chains[END_REF] Advanced NH + sectorial Kimoto and Ghil (1993b) Advanced NH = Northern Hemisphere; SH = Southern Hemisphere; PDF = probability density function.

no strong evidence in monthly mean reanalysis data for rejecting the single-regime multinormal hypothesis. The obvious reason for the latter result seems to be that the number of maps retained by these authors was too small or that the subspace they retained was too large for a satisfactory estimation of the pdf. Another way of looking at the difference between the results of [START_REF] Hannachi | Low-frequency nonlinearity and regime behavior in the Northern Hemisphere extratropical atmosphere[END_REF] and of [START_REF] Stephenson | On the existence of multiple climate regimes[END_REF] is the fact that the persistent anomalies associated with passage through the distinct regimes in phase space rarely lasts longer than ten days [START_REF] Dole | Persistent anomalies of the extratropical Northern Hemisphere wintertime circulation: Geographical distribution and regional persistence characteristics[END_REF]Kimoto and Ghil, 1993a,b). Thus, [START_REF] Hannachi | Low-frequency nonlinearity and regime behavior in the Northern Hemisphere extratropical atmosphere[END_REF] were separating such passages from the more diffuse flow in phase space outside the regimes, while [START_REF] Stephenson | On the existence of multiple climate regimes[END_REF] were not.

Besides the fairly low degree of controversy that persists on the very existence of multiple regimes, there is still some lack of agreement as to their causes. Thus, [START_REF] Lott | Mountain torques and Northern Hemisphere low-frequency variability. Part I: Hemispheric aspects[END_REF] find in a Hidden Markov Model (HMM) analysis that metastable regime transitions can occur despite a simple model's nearly Gaussian pdf, while [START_REF] Sura | Multiplicative noise and non-gaussianity: A paradigm for atmospheric regimes[END_REF] show that deviations from Gaussianity can also result from linear, stochastically perturbed dynamics with multiplicative noise statistics. [START_REF] Deremble | Multiple equilibria and oscillatory modes in a mid-latitude ocean-forced atmospheric model[END_REF] show that thermal forcing by a time-independent SST front can affect the atmospheric flow patterns and bifurcation sequence obtained in atmospheric models with likewise time-independent, but mechanical, i.e. topographic forcing [START_REF] Legras | Persistent anomalies, blocking and variations in atmospheric predictability[END_REF][START_REF] Jin | Intraseasonal oscillations in the extratropics: Hopf bifurcation and topographic instabilities[END_REF].

One may conclude that the number and variety of methods that have been used to identify and describe LFV regimes are leading up to a tentative consensus on their existence, robustness, and characteristics. We turn therewith to reviewing some of the theoretical foundations for the explanation of multiple regimes.

Theoretical basis of multiple regimes

Rossby wave propagation and interference. The slowly traveling, large-scale wave patterns that were first associated with weather phenomena in the 1930s are solutions of the partial differential equation for the conservation of potential vorticity q along a particle trajectory [START_REF] Gill | Atmosphere-Ocean Dynamics[END_REF][START_REF] Pedlosky | Geophysical fluid dynamics[END_REF]. For the purposes of this expository review, q can be defined as the vorticity ζ of a column of fluid divided by its height h, i.e., q = ζ/h. Conservation of q thus means, for instance, that a column of fluid's anti-clockwise rotation (defined as ζ > 0) will slow down (i.e., ζ decreases to smaller positive values) as the column moves over a mountain range (i.e., h > 0 decreases). This type of vorticity balance leads to slow Rossby waves [START_REF] Rossby | Relation between variations in the intensity of the zonal circulation of the atmosphere and the displacements of the semi-permanent centers of action[END_REF][START_REF] Haurwitz | The motion of atmospheric disturbances on the spherical earth[END_REF] that propagate westward with respect to the mean westerly jet. In this volume, [START_REF] Brunet | Identifying wave processes associated with predictability across time scales: An empirical normal mode approach[END_REF] present a more complete view of the role of such waves in S2S variability.

One view of persistent anomalies in mid-latitude atmospheric flows is that they result simply from the coincidental slowing down or linear interference of such Rossby waves [START_REF] Lindzen | Vacillations due to wave interference: applications to the atmosphere and to annulus experiments[END_REF][START_REF] Lindzen | Stationary planetary waves, blocking, and interannual variability[END_REF]. Another view is that a standing wave induced by topography can lead to a resonant interaction with two separate Rossby waves of distinct wavenumbers and thus produce a long-lived resonant wave triad [START_REF] Egger | Dynamics of blocking highs[END_REF]Ghil and Childress, 1987, Section 6.2). Neither one of these views provides an explanation of the observed clustering of persistent anomalies into distinct flow regimes. But the second one does suggest the more radically nonlinear theory described in the next few paragraphs. [START_REF] Charney | Multiple flow equilibria in the atmosphere and blocking[END_REF] took a major step in formulating a self-consistent atmospheric model for multiple equilibria and connecting it to observations of blocked and zonal flow. They used a highly idealized barotropic model to study the interaction between a zonal flow and simple topography with zonal wavenumber 2. Their model exhibits two stable equilibria for the same strength ψ * A of the prescribed zonal forcing that represents here the strength of the pole-to-equator temperature contrast. [START_REF] Charney | Comparison of a barotropic blocking theory with observation[END_REF] confronted the barotropic theory of [START_REF] Charney | Multiple flow equilibria in the atmosphere and blocking[END_REF] with observations, while [START_REF] Charney | Form-drag instability, multiple equilibria and propagating planetary waves in baroclinic, orographically forced, planetary wave systems[END_REF] extended it to baroclinic flows. In a somewhat complementary vein, [START_REF] Mitchell | Blocking-like solutions of the potential vorticity equation: their stability at equilibrium and growth at resonance[END_REF] followed up on the suggestion (e.g., [START_REF] Mcwilliams | An application of equivalent modons to atmospheric blocking[END_REF], and references therein) that steady-state solutions of the inviscid potential vorticity equation, written as q = G(ψ, p), can simulate blocking patterns, and may be resonantly excited by periodic forcing; here ψ = ψ(x, y, p) is the stream function, p is pressure, and (x, y) are horizontal coordinates pointing east-and northward.

Figure 1 shows the [START_REF] Charney | Multiple flow equilibria in the atmosphere and blocking[END_REF] model's bifurcation diagram, with the strength ψ A of the zonal jet in the model's steady-state solutions plotted against the corresponding strength ψ *

A of the forcing. The two stable equilibria-marked Z and R --are associated with "zonal" and "blocked" flow, respectively, ; see text for the explanation of abscissa and ordinate. The "S"-shaped bifurcation curve is typical of two back-to-back saddle-node bifurcations that give rise to two stable solution branches (solid) separated by an unstable one (dashed). (b) Flow patterns of the zonal (upper) and blocked (lower) branch, corresponding to the two stable equilibria Z and R (after [START_REF] Charney | Multiple flow equilibria in the atmosphere and blocking[END_REF]. (Reprinted with permission from Ghil andChildress (1987): Copyright 1987, Springer Nature.) as illustrated in Fig. 1b. The near-zonal solution is close in amplitude and spatial pattern to the forcing jet and is influenced very little by the topography, whereas the blocked solution is strongly affected by it. In the blocked-flow solution, a high-amplitude ridge is located upstream of the "mountains," similar to the situation during typical observed blocks off the West Coast of North America. This configuration, with a negative zonal pressure gradient on the windward slope of the mountains, corresponds to a negative mountain torque on the atmosphere. [START_REF] Benzi | The statistical properties of general atmospheric circulation: observational evidence and a minimal theory of bimodality[END_REF] and [START_REF] Hansen | The probability density distribution of the planetary-scale atmospheric wave amplitude revisited[END_REF] found evidence of bimodality in a composite index of wave amplitude in the NH mid-latitude flow. Although the statistical significance and robustness of their findings have been subject to criticism [START_REF] Nitsche | Is there evidence of multiple equilibria in planetary wave amplitude statistics[END_REF], direct confrontation of theoretical bimodality with observations has clearly stimulated LFV research during the 1980s. Further comparisons between the various approaches to LFV and to S2S prediction are drawn in Sec. 5.

A reasonable classification of low-pass filtered flow maps into discrete regimes provides only a static view of LFV. The next step is to study the transitions between these regimes over time. A matrix of probabilities for transitions from regime i to regime j is constructed by simply counting the transitions occurring in the data set. This yields an estimated set of conditional probabilities, in line with long-range forecasting experience [START_REF] Namias | Long-range weather forecasting: History, current status and outlook[END_REF][START_REF] Kalnay | Weather predictability beyond a week: an introductory review[END_REF] and the physical intuition that certain pathways of transition are more probable than others.

One kinematic approach to LFV is based on the Markov chain of these transitions. In this approach, knowledge of the system's present state is put to use to make a forecast, rather than using only unconditional probabilities. The Markov-chain view of LFV and, hence, long-range forecasting, is based on the existence of multiple regimes, the expected time of residence in each regime, and the probabilities of transition from one regime to another (Fig. 2).

Finally, several simple and intermediate model studies have addressed the issue of early warning indicators for regime transitions. For instance, [START_REF] Kondrashov | Weather regimes and preferred transition paths in a three-level quasigeostrophic model[END_REF] used the pdf of the exit angle from a given regime in a global baroclinic, quasi-geostrophic, 3-level (QG3) model with topography, originally due to [START_REF] Marshall | Toward a dynamical understanding of atmospheric weather regimes[END_REF], while [START_REF] Deloncle | Weather regime prediction using statistical learning[END_REF] used for the same model a fairly recent statistical learning method called random forests [START_REF] Breiman | Random forests[END_REF]. [START_REF] Kondrashov | Predicting weather regime transitions in Northern Hemisphere datasets[END_REF] then applied the random forests algorithm for studying regime transitions in Northern Hemisphere (NH) data sets, while [START_REF] Tantet | An early warning indicator for atmospheric blocking events using transfer operators[END_REF] used transfer operator theory for the study of transitions from zonal to blocking events in a hemispheric barotropic model. Methods for the study and prediction of regime transitions are discussed in greater detail in Sec. 4.2, after having provided additional insight into LFV dynamics in Sec. 3.2.

Extratropical oscillations in the S2S band

Phenomenological description of the extratropical oscillations in the S2S band. Variations of atmospheric angular momentum (AAM). Topographic instability and Hopf bifurcation; transition to irregular, chaotic behavior. Extratropical oscillations in a GCM with no Madden-Julian Oscillation (MJO) in the tropics.

Phenomenological description

Variations of geopotential height. In this subsection, we analyze the geopotential height field at 500 hPa (Z500) from the ERA-interim reanalysis of [START_REF] Dee | The ERA-interim reanalysis: configuration and performance of the data assimilation system[END_REF] that spans 37 years, January 1979-December 2016). We focus on the region 20 • S -90 • N, so as to include the tropics as well as the NH extratropics. This choice, along with that of the Z500 level, will become clearer as we proceed on to the results.

To identify spectral components of spatio-temporal behavior in the extensive ERA-interim reanalysis datasets, we apply here the multichannel singular spectrum analysis (M-SSA) methodology. For a comprehensive overview of the SSA and M-SSA methodology and of related literature, see Ghil et al. (2002) and (Alessio, 2016, chapter 12). M-SSA essentially diagonalizes the lag-covariance matrix of the multivariate dataset to yield a set of empirical orthogonal functions (EOFs) and the corresponding eigenvalues, which describe the variance captured by each EOF. In contrast to classical principal component analysis (PCA), M-SSA is able to capture oscillatory behavior in pairs of EOFs with approximately equal eigenvalues and dominant frequencies [START_REF] Vautard | Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series[END_REF][START_REF] Plaut | Spells of low-frequency oscillations and weather regimes in the Northern Hemisphere[END_REF]. To improve the separability of distinct frequencies, we rely here on a subsequent varimax rotation of the EOFs, cf. [START_REF] Groth | Multivariate singular spectrum analysis and the road to phase synchronization[END_REF].

Prior to the M-SSA analysis, the composite seasonal cycle is first subtracted from the daily time series. This cycle is computed by averaging the time series at each grid point over all years for each calendar day, and by smoothing this average using a 15-day running mean. Following Feliks et al. (2010, and references therein), we apply next a Chebyshev type I filter to the time series of anomaly maps, i.e., of the differences between the raw maps and the seasonal cycle so obtained, a lowpass filter that removes high-frequency oscillations with a period shorter than 20 days. The lowpass-filtered Z500 anomalies are subsampled with a sampling rate of 10 days, which yields time series of roughly 1300 samples. Finally, the subsampled anomalies are projected onto the 40 leading spatial EOFs of a classical PCA analysis, capturing 90% of the total variance. The corresponding principal components, finally, give the D = 40 input channels for the M-SSA analysis; see [START_REF] Groth | Interannual variability in the North Atlantic ocean's temperature field and its association with the wind stress forcing[END_REF] for further details of M-SSA analyzing extensive reanalysis datasets.

Oscillatory features in time and space. Figure 3(a) shows the eigenvalue spectrum of the subsampled Z500 anomalies. The spectrum exhibits several highly significant spectral peaks in the intraseasonal band; their frequencies agree very well with those identified and described in detail in previous studies [START_REF] Ghil | Intraseasonal oscillations in the global atmosphere. Part I: Northern hemisphere and tropics[END_REF][START_REF] Plaut | Spells of low-frequency oscillations and weather regimes in the Northern Hemisphere[END_REF]. In particular, two eigenvalue pairs stand out clearly above the red-noise spectrum at periods of 50 and 43 days, along with two shorter-period oscillations at 28 and 26 days. While [START_REF] Weickmann | Intraseasonal (30-60 day) fluctuations of outgoing longwave radiation and 250 mb streamfunction during northern winter[END_REF] discussed overall oscillatory behavior in the 30-60-day band, [START_REF] Dickey | Extratropical aspects of the 40-50 day oscillation in length-of-day and atmospheric angular momentum[END_REF] showed that variations in global atmospheric angular momentum (AAM) exhibit spectral peaks with distinct periods near 40 and 50 days. The latter authors demonstrated that the 50-day peak is largely associated with AAM fluctuations in the Tropics, while the 40-day peak is associated primarily with variations in the strength of the mid-latitude westerlies, particularly in the NH. Indeed, the amplitude of the 40-day oscillation in zonal winds is known to be largest during boreal winter, when the winds are strongest in the NH [START_REF] Weickmann | Intraseasonal (30-60 day) fluctuations of outgoing longwave radiation and 250 mb streamfunction during northern winter[END_REF][START_REF] Ghil | Intraseasonal oscillations in the global atmosphere. Part I: Northern hemisphere and tropics[END_REF][START_REF] Strong | Intraseasonal variability in a barotropic model with seasonal forcing[END_REF][START_REF] Strong | Intraseasonal oscillations in a barotropic model with annual cycle, and their predictability[END_REF] Our analysis of the Z500 anomalies goes into greater spatial detail than the zonally averaged one of [START_REF] Dickey | Extratropical aspects of the 40-50 day oscillation in length-of-day and atmospheric angular momentum[END_REF]. Herein, both oscillations, at 50 and at 43 days, make significant contributions to extratropical variability, with a Rossby wavetrain-like pattern across the North Atlantic extending into Eurasia, cf. Figs. 4(a,b). Figure 4(a), however, also shows a larger contribution of the 50-day oscillation to the dynamics in the Tropics, while a pattern that resembles the Pacific-North American (PNA) teleconnection extends northeastward from the Tropical Central Pacific; the 43-day oscillation, cf. Figure 4(b), is mainly limited to the NH, and contributes but little to the tropical LFV, except in the Atlantic sector.

Z = Zonal B = Blocked C = Wave Train τ Z τ Z τ B τ C
Furthermore, the 28-day and 26-day oscillations that show up as significant in Fig. 3(a) have dynamical patterns (not shown) that resemble a westward-traveling wave of this periodicity [START_REF] Branstator | A striking example of the atmosphere's leading traveling pattern[END_REF][START_REF] Kushnir | Retrograding wintertime low-frequency disturbances over the North Pacific ocean[END_REF][START_REF] Dickey | Extratropical aspects of the 40-50 day oscillation in length-of-day and atmospheric angular momentum[END_REF][START_REF] Ghil | Intraseasonal oscillations in the global atmosphere. Part I: Northern hemisphere and tropics[END_REF]. Figure 3(a) also points to a 71-day oscillation, whose period associates it with the 70-day oscillation emphasized by [START_REF] Plaut | Spells of low-frequency oscillations and weather regimes in the Northern Hemisphere[END_REF]. The 71-day oscillation found in Z500 here, however, falls below the stringent 95% confidence level of the other four peaks.

To visualize the dynamics underlying the variability pattern in Figs. 4(a,b), we have also calculated phase composites of the 43-day and 50-day mode. The phase composites of the 43-day mode are shown in Fig. 5. This mode's spatio-temporal structure is characterized by a retroregressing wave train, and it bears a remarkable resemblence to the 30-35-day oscillation of [START_REF] Plaut | Spells of low-frequency oscillations and weather regimes in the Northern Hemisphere[END_REF]. In phase category 1, we observe an amplification of the North Atlantic dipole structure, while in category 3, we observe a blocking pattern over Scandinavia. The overall pattern also shows greater variability in the Atlantic sector, with very little participation of the Pacific one.

The phase composites of the 50-day oscillation in Fig. 6 share certain features with the 43-day oscillations in the Atlantic sector, but exhibit much stronger variations in the Pacific sector. The 50-day tropical amplitude, although quite weak, suggests that the Pacific sector is much more involved in the tropicalextratropical interaction, whereas the 43-day mode in Fig. 5 is more focused on the Atlantic extratropical sector, in good agreement with the results of [START_REF] Plaut | Spells of low-frequency oscillations and weather regimes in the Northern Hemisphere[END_REF]. The present results thus concur, overall, with the conclusion of [START_REF] Plaut | Spells of low-frequency oscillations and weather regimes in the Northern Hemisphere[END_REF] that tropicalextratropical interactions in the intraseasonal band are much more active in the Pacific sector and with a 50-day periodicity, whereas the 43-day mode is more active in the Atlantic sector's extratropics. A similar distinction between the MJO and a separate 50-day mode in the Indian monsoon was made by [START_REF] Krishnamurthy | Intraseasonal oscillations of the monsoon circulation over South Asia[END_REF] and found to help S2S predictability there by [START_REF] Krishnamurthy | Predictability at intraseasonal time scale[END_REF]. Even so, [START_REF] Cassou | Intraseasonal interaction between the Madden-Julian Oscillation and the North Atlantic Oscillation[END_REF] showed that the MJO does affect two of the four weather regimes defined by [START_REF] Vautard | Multiple weather regimes over the North Atlantic: Analysis of precursors and successors[END_REF] in the North Atlantic sector and Western Europe. Cassou's result is consistent, still, with Fig. 4, in which we observe a significant contribution of the 50-day oscillation in the NH.

In Fig. 3(b), we have furthermore calculated the local variance [START_REF] Plaut | Spells of low-frequency oscillations and weather regimes in the Northern Hemisphere[END_REF][START_REF] Groth | Impacts of natural disasters on a dynamic economy[END_REF] for some of the intraseasonal oscillations, which shows a strong episodic modulation. This modulation illustrates one of the interesting relationships between the episodic, multiple-regime description of atmospheric LFV and the oscillatory one, based on intraseasonal oscillations (Kimoto and Ghil, 1993b;[START_REF] Plaut | Spells of low-frequency oscillations and weather regimes in the Northern Hemisphere[END_REF][START_REF] Koo | Multiple regimes and low-frequency oscillations in the Southern Hemisphere's zonal-mean flow[END_REF][START_REF] Crommelin | Regime transitions and heteroclinic connections in a barotropic atmosphere[END_REF][START_REF] Crommelin | Observed nondiffusive dynamics in large-scale atmospheric flow[END_REF][START_REF] Kondrashov | Weather regimes and preferred transition paths in a three-level quasigeostrophic model[END_REF]. The complementarity of these two points of view, episodic and oscillatory (Ghil and Robertson, 2002), will be discussed in Sec. 5.

Topographic instability and Hopf bifurcation

We outline here how a hierarchy of models [START_REF] Ghil | Hilbert problems for the geosciences in the 21st century[END_REF]Ghil andRobertson, 2000, 2002) can be used to formulate and test the hypothesis that the 40-day oscillation is an intrinsic mode of the NH extratropics, associated with the interaction between the jet stream and mid-latitude mountain ranges. The rudiments of this hypothesis originate in the highly idealized barotropic model of [START_REF] Charney | Multiple flow equilibria in the atmosphere and blocking[END_REF], which we discussed in Sec. 2.2.

More complex models-both barotropic (i.e., single-layer) and baroclinic (i.e., multilayer), with more spatial degrees of freedom than the [START_REF] Charney | Multiple flow equilibria in the atmosphere and blocking[END_REF] model-exhibit multiple flow patterns, some of which are similar to the blocked and zonal ones found in Fig. 1 above, even for fairly realistic values of the forcing. The crucial difference is that the equilibria found in the more complex models are no longer stable, and the system oscillates around the blocked solution or fluctuates between the zonal and blocked solutions in an irregular, chaotic way [START_REF] Legras | Persistent anomalies, blocking and variations in atmospheric predictability[END_REF]Ghil and Childress, 1987, ch. 6).

The [START_REF] Legras | Persistent anomalies, blocking and variations in atmospheric predictability[END_REF] model, for instance, has 25 spherical harmonics and its bifurcation diagram also appears as Fig. 8 in Ghil and Robertson (2002). In this diagram, the branch of blocked equilibria is destabilized by an oscillatory instability as the intensity of the forcing jet increases. The transition from a stable stationary solution -called a fixed point in dynamical systems theory -to a stable periodic solution, called a limit cycle in this theory, is termed a Hopf bifurcation. The limit cycle that arises from this bifurcation increases in amplitude as the forcing jet becomes stronger, and it has a period of roughly 40 days. This limit cycle loses its stability in turn as forcing is increased further and the flow becomes chaotic, but the 40-day periodicity survives in the power spectrum of the chaotic flow. [START_REF] Jin | Intraseasonal oscillations in the extratropics: Hopf bifurcation and topographic instabilities[END_REF] showed that, when a sufficiently realistic meridional structure of the solutions' zonal jet is allowed, the back-to-back saddle-node bifurcations of Fig. 1a are indeed replaced by a Hopf bifurcation and thus transition to finite-amplitude periodic solutions can occur. Eigenanalyses of the unstable equilibria in a barotropic model with higher horizontal resolution, as well as its time-dependent solutions, indicate oscillatory instabilities with intraseasonal (35-50 days) and biweekly (10-15 days) time scales [START_REF] Strong | Intraseasonal variability in a barotropic model with seasonal forcing[END_REF]. Floquet analysis of this model's limit cycles [START_REF] Strong | Intraseasonal oscillations in a barotropic model with annual cycle, and their predictability[END_REF] confirms that the 40-day oscillations that arise in it by oscillatory topographic instability are stronger in winter than in summer, like the observed NH oscillations [START_REF] Ghil | Intraseasonal oscillations in the global atmosphere. Part I: Northern hemisphere and tropics[END_REF][START_REF] Knutson | 30-60 day atmospheric oscillations: Composite life cycles of convection and circulation anomalies[END_REF].

To test the theory of NH extratropical oscillations developed in simpler models, [START_REF] Marcus | The extratropical 40-day oscillation in the UCLA general circulation model. Part I: Atmospheric angular momentum[END_REF][START_REF] Marcus | The extratropical 40-day oscillation in the UCLA general circulation model. Part II: Spatial structure[END_REF] studied a 3-year perpetual-January simulation that had been performed with a version of the UCLA GCM in which no self-sustained MJO was apparent in the tropics. A robust 40-day oscillation in AAM is found to arise in the model's NH extratropics when standard topography is present. Three shorter runs with no topography produced no intraseasonal oscillation; this result is consistent with a topographic origin for the NH extratropical oscillation in the standard model. The spatial structure of the circulation anomalies associated with the model's extratropical oscillation is shown in an 8-minute video, https://AAM_the_Video, as displayed in several model-simulated fields, namely 500-hPa geopotential heights, 250-hPa streamfunction, and surface pressure torques.

The oscillation is dominated by a standing wavenumber-2 pattern, which undergoes a predominantly barotropic, tilted-trough vacillation. High values of AAM are associated with low 500-hPa heights over the northeast Pacific and the North Atlantic oceans, and vice versa. These flow patterns resemble the configurations seen in the [START_REF] Charney | Multiple flow equilibria in the atmosphere and blocking[END_REF] simple model (see Fig. 1b here). The NE-SWtilting and NW-SE-tilting phases in the video's 500-hPa fields are strongly reminiscent of the extremes and intermediate phases of the 40-day oscillation that arises by Hopf bifurcation from the blocked equilibrium in the [START_REF] Legras | Persistent anomalies, blocking and variations in atmospheric predictability[END_REF] model; see discussion of Ghil et al. (2003, Fig. 3).

Low-order data-driven modeling, dynamical analysis, and prediction

Linear and nonlinear LOMs and the role of memory effects. Empirical model reduction (EMR) & EMR-based prediction. Partial observations and stochastic closure in LOMs. Slow changes in the forcing and random effects: non-autonomous and random dynamical systems.

Background and methodological LOM developments

State-of-the-art, highly resolved GCMs, while able to simulate detailed interactions within the climate system over a wide range of scales, generate detailed 4-D climate variability that is visually as complex as currently available observational datasets, and is hence no less challenging to interpret. Dynamical analysis of climatic phenomena typically involves a set of multiple GCM simulations that are designed to isolate physical processes governing the simulated, and by inference, observed climate variability. These simulations are, however, computationally expensive, and their interpretation is hindered by the presence of model biases due to incomplete or imperfect parameterizations of the unresolved physical processes.

Moreover, GCMs represent a broad range of time and space scales and use a state vector that has many millions of scalar variables. While detailed weather prediction out to a few days does require such high numerical resolution, considerable work, both theoretical (Ghil and Childress, 1987, Sec. 6.5) and databased [START_REF] Toth | Degrees of freedom in Northern Hemisphere circulation data[END_REF], has shown that important aspects of observed atmospheric LFV can be represented and predicted with a substantially smaller number of degrees of freedom, as encapsulated by LOMs.

This statement extends beyond the issues of simulating and understanding LFV: it clearly includes prediction as well, and it is also consistent with the fact that climate concerns statistical time and space aggregates of weather. In fact, while there has been a recent surge of interest in S2S prediction through coordinated international modeling activity, most of this activity involves GCM-based ensemble prediction systems. Empirical, data-driven prediction, though, has a long history in forecasting the El Niño-Southern Oscillation (ENSO) and the best LOMs (e.g., [START_REF] Kondrashov | A hierarchy of data-based ENSO models[END_REF] were still outperforming -no longer than six years ago -most GCM-based ENSO forecasts, even for lead times beyond the S2S range [START_REF] Barnston | Skill of real-time seasonal ENSO model predictions during 2002-2011 -is our capability improving[END_REF]. Unfortunately, since the [START_REF] Barnston | Skill of real-time seasonal ENSO model predictions during 2002-2011 -is our capability improving[END_REF] review paper, the plume of real-time Niño-3.4 predictions -which includes those of stochastic-dynamic models, as well as of GCMs -that is continuously monitored at the International Research Institute for climate and society (IRI) no longer allows one to isolate the skill of individual LOMs or GCMs. Still, inspection of the existing documentation at https://iri.columbia.edu/our-expertise/climate/forecasts/enso/ 2017-October-quick-look/?enso_tab=enso-sst_table does indicate that -over the last seven boreal winters, 2012-13 to 2017-18 -ENSO prediction has remained a tough problem for both GCMs and simpler models. Hence, the contribution of the latter to real-time S2S forecasting can still be quite useful, in and of itself, as well as in stimulating improvements in GCM predictions.

A theoretical approach to LOM development relies on time scale separation in the full governing equations to derive a reduced set of differential equations for the LFV and to parameterize the fast, unresolved variability -such as atmospheric convection on synoptic and mesoscales -by stochastic forcing (e.g., [START_REF] Epstein | Stochastic dynamic prediction[END_REF][START_REF] Fleming | On stochastic dynamic prediction: I. The energetics of uncertainty and the question of closure[END_REF][START_REF] Majda | Models for stochastic climate prediction[END_REF]). An alternative approach to climate diagnosis and prediction is based on empirical modeling of the LFV, as observed in selected fields, under the assumption that it can be described as a spatio-temporal, nonlinear and stochastic process. This statistical data-driven approach does not require knowledge of the governing equations: even though it lacks the immediate dynamical interpretability of climate models that are derived from first principles, it does allow one to reproduce detailed aspects of the observed statistics by inverse modeling techniques, albeit in the subset of fields being simulated by the LOM. The available datasets -whether given by direct observations, reanalyses or GCM simulations -are used to estimate both the model's low-order, deterministic part and its driving noise.

Development of data-driven LOMs can thus be cast as a closure problem: These models can simulate and predict LFV provided they are able to properly account for (a) linear and nonlinear interactions between the resolved components, i.e., the high-variance LFV modes; and (b) the interactions between the resolved components and the huge number of unresolved ones that represent the unobserved small-scale processes and are not explicitly included in the LOM. The key steps are thus estimating interactions between the macroscopic variables, identifying the hidden variables, and modeling the cross-interactions between the macroscopic and hidden variables.

Linear inverse models (LIMs) and Principal Oscillation Pattern (POP) analysis [START_REF] Penland | Random forcing and forecasting using principal oscillation pattern analysis[END_REF][START_REF] Penland | Forecasting Northern Hemisphere 700-mb geopotential height anomalies using empirical normal modes[END_REF][START_REF] Penland | A stochastic model of IndoPacific sea surface temperature anomalies[END_REF][START_REF] Penland | The optimal growth of tropical sea surface temperature anomalies[END_REF] assume a linear dynamical model for the macroscopic variables, while their interactions with the small-scale, hidden processes are approximated by spatially correlated white noise. The Empirical Model Reduction (EMR) approach [START_REF] Kravtsov | Multi-level regression modeling of nonlinear processes: Derivation and applications to climatic variability[END_REF][START_REF] Kravtsov | Empirical model reduction and the modeling hierarchy in climate dynamics and the geosciences[END_REF][START_REF] Kondrashov | Low-order stochastic model and "past-noise forecasting" of the Madden-Julian oscillation[END_REF] generalizes LIMs by allowing (i) quadratic interactions among macroscopic variables; (ii) time-delayed dynamical effects via memory terms; and (iii) a richer temporal structure of the noise. The memory effects that may appear in both the stochastic forcing and the dynamical operator are conveyed by hidden variables arranged into a stacked system of levels. Each additional level includes a new hidden variable, which is less auto-correlated than the one included on the previous level, until the last-level variable can be approximated by spatially correlated white noise.

Chen et al. ( 2016) conducted a comprehensive suite of EMR experiments using a 4000-yr-long preindustrial control simulation dataset from the CM2.1 coupled model of NOAA's Geophysical Fluid Dynamics Laboratory (GFDL), with the goal of better understanding ENSO diversity, nonlinearity, seasonality and the effects of memory in the simulation and prediction of the Tropical Pacific's SST anomalies. The results show that multilevel nonlinear EMR models that account for SST history improve ENSO prediction skill dramatically.

The EMR model coefficients introduced on each level are estimated by multilevel regression techniques; see [START_REF] Kravtsov | Empirical model reduction and the modeling hierarchy in climate dynamics and the geosciences[END_REF] for a comprehensive review of the EMR methodology. [START_REF] Strounine | Reduced models of atmospheric low-frequency variability: Parameter estimation and comparative performance[END_REF] have systematically compared various model reduction methods, as applied to wintertime geopotential height anomalies, and have demonstrated that, in the absence of clear scale separation, EMR methodology's success is rooted in its multilayer structure that accounts for the memory effects needed to achieve optimal closure [START_REF] Kravtsov | Multi-level regression modeling of nonlinear processes: Derivation and applications to climatic variability[END_REF]. [START_REF] Newman | A study of subseasonal predictability[END_REF] showed that, at the beginning of the past decade, the forecast skill of a LIM model of NH extratropical, weekly averaged circulation is comparable to that of NCEP's medium-range forecast (MRF) model at week 2 (days 8-14) and week 3 (days 15-21). [START_REF] Zhang | Cracking the MJO nut[END_REF] and [START_REF] Vitart | Madden-Julian oscillation prediction and teleconnections in the S2S database[END_REF] have documented the progress in MJO simulation and prediction by GCMs since then. Thus, the best GCM predictions have achieved a bivariate correlation skill near 0.6 at four weeks, as measured between the two observed and forecast Real-Time Multivariate MJO Indices known as RMM1 and RMM2.

At the same time, LOMs have also progressed beyond LIMs, as outlined above. [START_REF] Kondrashov | Low-order stochastic model and "past-noise forecasting" of the Madden-Julian oscillation[END_REF] have improved upon the EMR-based predictions used in real-time ENSO forecasting [START_REF] Kondrashov | A hierarchy of data-based ENSO models[END_REF] by including in the forecast ensemble information on the "weather noise" that prevailed at particular phases of the MJO cycle, via the past-noise forecasting (PNF) methodology of Chekroun et al. (2011b). These authors have demonstrated that, in retrospective forecasting, PNF-enhanced EMR methodology yields bivariate-RMM correlation results that are quite comparable to those of the best GCMs (Kondrashov et al., 2013, Fig. 2).

More broadly speaking, EMR-based climate predictions have proven to be highly useful and competitive as statistical benchmarks against physics-based models for prediction on the S2S time scales associated with the MJO [START_REF] Kondrashov | Low-order stochastic model and "past-noise forecasting" of the Madden-Julian oscillation[END_REF], as well as with ENSO [START_REF] Kondrashov | A hierarchy of data-based ENSO models[END_REF].

The Multilayer Stochastic Model (MSM) framework introduced by Kondrashov et al. (2015) establishes, in a data-driven context, a clear connection between the EMR methodology and the Mori-Zwanzig (MZ) formalism of statistical physics. It is this sound physical basis that helps derive a closed system of stochasticdynamic equations for a subset of the original system's variables (Kondrashov et al., 2015, Secs. 4 & 5). This closed system of equations includes Markovian and non-Markovian deterministic terms that are not necessarily quadratic, as well as stochastic noise terms [START_REF] Chorin | Optimal prediction with memory[END_REF][START_REF] Chorin | Stochastic Tools in Mathematics and Science[END_REF].

The major methodological problem is that, in the highly nonlinear climate system, there is typically a continuum of scales, from the fastest and shortest to the slowest and longest. In this case, the MZ formalism allows one to treat the interactions between LFV and smaller scales as memory effects. These memory effects represent a significant deviation from the markovianity of the dynamical models that are most often used in modeling atmospheric, oceanic and climate phenomena and processes. Markovianity means dependence on the initial state alone, as opposed to states that preceded it. Thus, ordinary and partial differential equations, as well as the Markov chains mentioned in Sec. 2 herein, are markovian: delay-differential and other functional differential equation [START_REF] Hale | Introduction to Functional Differential Equations[END_REF] are not. [START_REF] Bhattacharya | Internal variability of an energy-balance model with delayed albedo effects[END_REF] introduced the latter into the climate sciences and they have been used since extensively in ENSO modeling, with memory terms that account for travel times of Kelvin and Rossby waves across the Tropical Pacific; see also (Ghil et al., 2015, and references therein).

The non-Markovian terms become particularly important in the absence of scale separation, whereas [START_REF] Stinis | A comparative study of two stochastic mode reduction methods[END_REF] showed that, for relatively large scale separation, the results of the MZ methodology are similar to those of the model reduction formalism of [START_REF] Majda | Models for stochastic climate prediction[END_REF]. The MSM framework thus improves upon EMR by its greater generality of the deterministic terms, as well as by the proper understanding of the memory effects. The improvements in simulation and prediction achieved by nonlinear (EMR and MSM) vs. linear (LIM and POP) data-driven models (e.g., [START_REF] Barnston | Skill of real-time seasonal ENSO model predictions during 2002-2011 -is our capability improving[END_REF] seem to further buttress the credibility of the highly nonlinear Devil's staircase scenario for ENSO [START_REF] Jin | El Niño on the devil's staircase: Annual subharmonic steps to chaos[END_REF][START_REF] Tziperman | El Niño chaos: Overlapping of resonances between the seasonal cycle and the Pacific ocean-atmosphere oscillator[END_REF] vs. stable noise-forced dynamics [START_REF] Neelin | ENSO theory[END_REF]. [START_REF] Kondrashov | Stochastic modeling of decadal variability in ocean gyres[END_REF] have shown that a change of basis by M-SSA can help in applying the MSM methodology to capture the LFV in multiscale oceanic turbulence, because M-SSA utilizes time-lagged information and therefore implicitly conveys the memory effects. Data-adaptive harmonic decomposition (DAHD: [START_REF] Chekroun | Data-adaptive harmonic spectra and multilayer Stuart-Landau models[END_REF] further advances the dynamic interpretation of nonlinear datadriven models. DAHD provides, for a broad class of time-evolving datasets, reduction coordinates that can be efficiently modeled by systems of paired stochastic differential equations (SDEs), while using a fixed set of predictor functions. The methodology has proven itself already in the modeling, simulation and prediction of the Arctic's sea ice extent and concentration field (Kondrashov et al., 2017a,b), but remains to be tested for other climatic time series.

In summary, the field of low-order, data-driven stochastic-dynamic models has undergone remarkable development in recent years. These models offer, therewith, new opportunities to further our understanding and prediction on S2S time scales. The following two subsections present examples of their successful application to mid-latitude LFV variability and to ENSO phenomena, respectively.

Dynamical diagnostics and empirical prediction on S2S scales

LOMs can be very useful for the dynamical interpretation of atmospheric oscillations and multiple regimes. Such models not only help compact the dataset's information content but can also provide insights into the dynamics and predictability of climatic LFV, via the analysis of the reduced model's mathematical structure. Thus, [START_REF] Kondrashov | Empirical mode reduction in a model of extratropical low-frequency variability[END_REF][START_REF] Kondrashov | Signatures of nonlinear dynamics in an idealized atmospheric model[END_REF] constructed and analyzed a LOM of extratropical atmospheric LFV by applying EMR to the output of a long simulation of the QG3 model mentioned in Sec. 2.2. This model has a fairly realistic climatology and variability, and it has been used quite extensively to study NH mid-latitude flows (e.g., D'Andrea and [START_REF] D'andrea | Extratropical low-frequency variability as a low-dimensional problem I: A simplified model[END_REF][START_REF] D'andrea | Extratropical low-frequency variability as a low-dimensional problem. II: Stationarity and stability of largescale equilibria[END_REF]Deloncle et al., 2007, and references therein). The full QG3 model's phase space has a dimension of O(10 4 ), while the derived EMR model has only 45 variables. The regimes of both the EMR model and the QG3 model were computed using a Gaussian mixture model by following Smyth et al. (1999, and references therein). The centroids of the four regimes for the full QG3 model are plotted in Fig. 7, and they correspond to opposite phases of the Arctic Oscillation (AO + and AO -) and of the North Atlantic Oscillation (NAO + and NAO -), in good agreement with the earlier results of [START_REF] Molteni | Long-term fluctuations in the statistical properties of low-frequency variability: Dynamical origin and predictability[END_REF]. The EMR model captures very well the non-Gaussian features of the full QG3 model's probability density function (PDF) and shares with the latter its four anomalously persistent flow patterns, as well as the Markov chain of transitions between these regimes (not shown).

180 ° W 150 ° W 120 ° W 90 ° W 60 ° W 30 ° W 0 ° 30 ° E 60 ° E 90 ° E 120 ° E 150 ° E 180 ° E a) NAO + 180 ° W 150 ° W 120 ° W 90 ° W 60 ° W 30 ° W 0 ° 30 ° E 60 ° E 90 ° E 120 ° E 150 ° E 180 ° E b) NAO - 180 ° W 150 ° W 120 ° W 90 ° W 60 ° W 30 ° W 0 ° 30 ° E 60 ° E 90 ° E 120 ° E 150 ° E 180 ° E c) AO + 180 ° W 150 ° W 120 ° W 90 ° W 60 ° W 30 ° W 0 ° 30 ° E 60 ° E 90 ° E 120 ° E 150 ° E 180 ° E d) AO -
The spatial correlations between the anomaly patterns in Fig. 7 and those of the EMR model (not shown) all exceed 0.9. This good match between the regime centroids in the two models, full and reduced, as well as between the corresponding Markov chains, emphasizes the intrinsic dynamical nature of both the NAO and the AO, and strongly suggests that the NAO is not merely a regionalization of the AO, at least not in the fairly realistic QG3 model.

The four regimes that are best supported by synoptic experience, as well as by the statistical analysis of the upper-air data for the past half-century, are the zonal and blocked phases of westerly flow in the Atlantic-Eurasian and Pacific-North-American sector, respectively [START_REF] Cheng | Cluster analysis of the Northern Hemisphere wintertime 500-hPa height field: Spatial patterns[END_REF][START_REF] Smyth | Multiple regimes in Northern Hemisphere height fields via mixture model clustering[END_REF]Ghil and Robertson, 2002). The AO, also called the NH annular mode [START_REF] Wallace | North Atlantic Oscillation/annular mode: Two paradigms-one phenomenon[END_REF], seems to be only statistically annular: it representa a redistribution of mass between the poles and subtropics, but actual "sloshing" events are sectorial, by-and-large. Hence the correlation between the subtropical highs over the North Atlantic and the North Pacific in observational and reanalysis data is quite low at sub-monthly time scales (e.g., Kimoto and Ghil, 1993b).

Both the AO and NAO indices are, however, dominated by Arctic sea-level pressure variations, and hence the NAO and AO are highly correlated in observations. Moreover, both the QG3 model [START_REF] Kondrashov | Empirical mode reduction in a model of extratropical low-frequency variability[END_REF] and the no-MJO version of the UCLA GCM [START_REF] Marcus | The extratropical 40-day oscillation in the UCLA general circulation model. Part I: Atmospheric angular momentum[END_REF][START_REF] Marcus | The extratropical 40-day oscillation in the UCLA general circulation model. Part II: Spatial structure[END_REF] seem to possess some form of a 40-day sloshing mode that originates in mid-latitudes. Seeing the hemispheric AO + and AO -patterns of Figs. 7(c,d) replace the PNA-sector patterns of Ghil and Robertson (2002, and references therein) is thus not too surprising: the hemispheric AO ± centroids overlap largely in the Atlantic-Eurasian sector with NAO ± , respectively, from which they differ mainly by much stronger centers of action in the PNA sector. In addition, M-SSA identifies intraseasonal oscillations with a period of 35-37 days and of 20 days in the data generated by both the full QG3 model and its EMR version. These oscillatory modes are similar to those seen in the observations, as plotted in Fig. 3 here, and in previous studies (e.g., Ghil andRobertson, 2000, 2002, andreferences therein). The former one is clearly these models' version of the extratropical, 40-day mode found by [START_REF] Legras | Persistent anomalies, blocking and variations in atmospheric predictability[END_REF], while the QG3 model is expected to lack the MJO-related 50-day mode seen in Fig. 3. [START_REF] Lott | Mountain torques and atmospheric oscillations[END_REF]Lott et al. ( , 2004a,b) ,b) have shown that observational NH data confirm the importance of mountain torques in an extratropical oscillation that has a shorter periodicity than the MJO. Mountain torque anomalies were also found to anticipate the phases of the AAM oscillation, as well as the onsets and breaks of certain flow regimes.

EOF-1 EOF-2 a) -2 0 2 -2 -1 0 1 2 EOF-1 EOF-2 b) -2 0 2 -2 -1 0 1 2 EOF-1 EOF-2 c) -2 0 2 -2 -1 0 1 2 EOF-1 EOF-2 d) -2 0 2 -2
The QG3 model's LFV features can be interpreted via dynamical analysis of the reduced model. In particular, the AO -regime arises from the unique steady state of the EMR model's three-level deterministic operator, at the larger positive values of PC-1, while the PDF ridge that is present in the QG3 model, cf. Kondrashov et al. (2006, Fig. 6), coincides with the location of a plateau of quasi-stationary states of the reduced model, cf. here. The increase of the extent of the PDF plateau in both the full and the reduced model -first along the EOF-1 direction and then, at lower negative values of PC-1, along the EOF-2 direction -appears to be due to preferential instabilities in the subspace of the QG3 model's leading modes. These instabilities are to be equilibrated by the full model's transient feedbacks that are captured in the EMR by the noise and memory effects.

Changes in the number and amplitude of the PDF's modes were also observed by [START_REF] Molteni | Long-term fluctuations in the statistical properties of low-frequency variability: Dynamical origin and predictability[END_REF] when varying the potential vorticity forcing in the QG3 model to represent a subtropical Rossby wave source appropriate for cold ENSO events. The changes plotted in Fig. 8 here can be interpreted as changes that might arise in a GCM subject either to stochastic parametrizations with varying noise intensities (e.g., Palmer and Williams, 2009, and references therein) or to fluctuating changes in external forcing. The rather surprising succes of EMRs, in their various guises, at simulating and predicting S2S variability may be largely attributed to the presence and proper intensity of the noise processes they include.

The dominant intraseasonal oscillation in both the QG3 and the optimal EMR model has a period of about 35-37 days and it is associated with the least-damped eigenmode of the latter, when linearized about its climatological state. While there is no clear scale separation in the QG3 simulations, it is the interactions between the model's largest-scale modes -whose variance is concentrated in the subspace of its four leading EOFs -and the intermediate scales, captured by EOFs 5-15, that appear to be responsible for the slowdown of the intraseasonal oscillation's trajectory in the full QG3 model, as well as in the optimal EMR model. This slow-down is associated with the emergence of the AO + and NAO -regimes.

The EMR model's stochastic forcing represents the unresolved smaller-scale processes in the QG3 model. The effect of this forcing on the full dynamics is revealed by a bifurcation analysis with respect to the noise's amplitude . In particular, Fig. 8 shows that, as increases from zero to its optimal value, the model trajectory is initially confined to a small region near the EMR model's unique steady state (panel (a)) and that it gradually fills up the quasi-stationary ridge along the EOF-1 axis (panels (b)-(d)). When becomes large enough, the intraseasonal oscillatory mode is excited, resulting in the model's PDF expanding in the EOF-2 and EOF-3 directions, while additional quasi-stationary states appear. The sequence of panels in Fig. 8 can be taken to imply that the stochastic forcing needs to be sufficiently strong in order to obtain realistic LFV regimes. The role of stochastic processes in nonlinear dynamics will be discussed more fully at the end of this chapter, in Sec. 5. [START_REF] Kondrashov | Signatures of nonlinear dynamics in an idealized atmospheric model[END_REF] have also analyzed this EMR model's nonlinear dynamical operator and shown that it gives rise to "swirls" in the time-mean tendencies of the model's state-averaged trajectories, when projected onto a low-dimensional subspace. These swirls are dominated by the deterministic interactions between the large-scale modes. More recently, [START_REF] Tantet | An early warning indicator for atmospheric blocking events using transfer operators[END_REF] have associated similar swirl patterns with persistent weather regimes in the output of a barotropic atmospheric model with an intermediate resolution of 231 spherical modes. These authors applied transfer operator methods (Chekroun et al., 2014, and references therein) to develop early warning indicators for regime transitions, in particular from zonal to blocked flows, as mentioned in Sec. 2.2.

The idea of considering preferred directions of instability of a regime's centroid goes back to the "ghost equilibria" of [START_REF] Legras | Persistent anomalies, blocking and variations in atmospheric predictability[END_REF]. In particular, [START_REF] Deloncle | Weather regime prediction using statistical learning[END_REF] and [START_REF] Kondrashov | Predicting weather regime transitions in Northern Hemisphere datasets[END_REF] have applied the random forests algorithm [START_REF] Breiman | Random forests[END_REF] to learn the predictors connected to these preferred directions in a low-dimensional space. This statistical approach worked rather well in medium-to-long-range prediction of regime breaks, in both the QG3 model [START_REF] Deloncle | Weather regime prediction using statistical learning[END_REF] and NH atmospheric observational data [START_REF] Kondrashov | Predicting weather regime transitions in Northern Hemisphere datasets[END_REF]. The results of [START_REF] Tantet | An early warning indicator for atmospheric blocking events using transfer operators[END_REF] for why that appears to be the case. In the context of S2S prediction, the dynamically motivated predictors might turn out to correspond to preferred instability directions of either fixed points that "anchor" multiple regimes, such as El Niño and La Niña (Dijkstra, 2005, Ch. 7, and references therein), or to similar instabilities of limit cycles that are associated with oscillatory modes, such as the MJO or the extratropical, 40-day mode discussed in Sec. 3.1 and the present one. Instability of limit cycles in atmospheric LFV was studied, for instance, by [START_REF] Strong | Intraseasonal oscillations in a barotropic model with annual cycle, and their predictability[END_REF], and is mentioned here in Sec. 3.2.

LFV and multilayer stochastic closure: A simple illustration

In this subsection, we provide an illustrative example that will help understand the reasons for the effectiveness of MSM modeling, and hence of EMR models, in simulating and potentially predicting extratropical S2S variability. This example illustrates intraseasonal ENSO prediction.

Consider the following nonlinear, periodically and stochastically forced two-variable (x 1 , x 2 ) model from Chekroun et al. (2011b, Supplementary Information):

   dx 1 = {(r + σdW t )x 1 (α + x 1 )(1 -x 1 ) -cx 1 x 2 + a sin(2πf 0 t)}dt, dx 2 = {-mαx 2 + (c -m)x 1 x 2 }dt. (1)
Here the variable x 1 is subject to deterministic, additive forcing a sin(2πf 0 t), while the stochastic term σdW t represents a random perturbation of the parameter r by white noise, and we assume that only the variable x 2 is observed. The default values of the model parameters are σ = 0.3, m = r = 1, c = 1.5 and α = 0.3. The model is integrated from time t = 0 to t = T f = 2000 (in dimensionless units), with the initial state (x 1 (0), x 2 (0)) = (0.5, 0.5), by using a classical stochastic Euler-Maruyama scheme with step size ∆t = 0.1. When both periodic forcing and noise are absent, the model has only one globally stable equilibrium.

When turning on the periodic forcing with amplitude a = 0.05 and frequency f 0 = 0.25, the system exhibits only one periodic orbit of period 4, which is globally stable. In the presence of noise, with σ = 0.3, an LFV mode with a period equal to approximately 25 units -i.e., a frequency of f = 0.04 = f 0 -becomes dominant, as seen in Fig. 9. This mode (blue in the figure) is captured by SSA reconstruction, and it can be attributed to a damped non-normal mode that is sustained by the noise (Chekroun et al., 2011b). In the present ENSO-related illustration, the periodic forcing T 0 = 1/f 0 = 4 corresponds to the seasonal cycle, while the intrinsic periodicity of T 25 is associated with ENSO's low-frequency mode of roughly 6 years.

By applying EMR solely to a time series of the observed variable x 2 , cf. [START_REF] Kravtsov | Multi-level regression modeling of nonlinear processes: Derivation and applications to climatic variability[END_REF], a univariate two-level EMR model is obtained for x 2 alone. This scalar model includes a cubic polynomial, as well as multiplicative interaction with the periodic forcing, on its main level (not shown). Note that this parametric form is entirely different from the correct form of the right-hand side of the x 2 -equation in Eq. ( 1), as it does not explicitly include interactions with the periodically and stochastically forced unobserved variable x 1 . These interactions are parameterized by the memory effects conveyed in the two-level EMR model.

Prediction skill is plotted in Fig. 10 in terms of both root-mean-square errors (RMSE) and anomaly correlations (Corr). The figure shows that the out-of-sample prediction skill of the ensemble mean for the estimated EMR model is significantly better than that of a linear single-level LIM model, with the EMR's maintaining useful prediction skill, i.e. Corr≥ 0.5, up to a lead time of roughly 6 time units, which represents roughly a quarter of a period of the LFV mode.

This example illustrates the benefits for prediction of including nonlinear and memory terms in empirical stochastic models, as demonstrated by comparing the performance of UCLA's EMR model in real-time ENSO forecasting with other statistical forecasts; see [START_REF] Barnston | Skill of real-time seasonal ENSO model predictions during 2002-2011 -is our capability improving[END_REF].

Concluding remarks

Considerable progress has been made in the 15 years since the review article of Ghil and Robertson (2002) in low-order predictive modeling of what was called at that time extended or long-range forecasting and is now more specifically referred to as S2S prediction. Four of the approaches that were considered in this modeling are summarized in Fig. 11.

One approach to persistent anomalies in mid-latitude atmospheric flows on 10-100-day time scales is to consider them simply as due to slowing down of Rossby waves or to their linear interference [START_REF] Lindzen | Vacillations due to wave interference: applications to the atmosphere and to annulus experiments[END_REF][START_REF] Lindzen | Stationary planetary waves, blocking, and interannual variability[END_REF]. This approach is illustrated in the sketch labeled (c) within the figure: zonal flow Z and blocked flow B are simply slow phases of an harmonic oscillation, like the neighborhood of t = π/2 or t = 3π/2 for a sine wave sin(t); or else they are due to an interference like that occurring for a sum A sin(t) + B sin(3t) near t = (2k + 1)π/2. A more thorough, quasi-linear version of this approach is to study long-lived resonant wave triads between a topographic Rossby wave and two free Rossby waves [START_REF] Egger | Dynamics of blocking highs[END_REF][START_REF] Trevisan | Stationary response of barotropic weakly non-linear Rossby waves to quasi-resonant orographic forcing[END_REF]Ghil and Childress, 1987, Section 6.2). Neither version of this approach, though, explains the anomalies' organizing into distinct flow regimes. Rossby et al. ( 1939) initiated a different, genuinely nonlinear approach by raising the possibility of multiple equilibria as an explanation of preferred atmospheric flow patterns. These authors drew an analogy between such equilibria and hydraulic jumps, and formulated simple models in which similar transitions between faster and slower atmospheric flows could occur. This multiple-equilibria approach was then pursued quite agressively in the 1980s [START_REF] Charney | Multiple flow equilibria in the atmosphere and blocking[END_REF][START_REF] Charney | Comparison of a barotropic blocking theory with observation[END_REF][START_REF] Legras | Persistent anomalies, blocking and variations in atmospheric predictability[END_REF]Ghil and Childress, 1987, Sections 6.3-6.6) and it is illustrated in Fig. 11 by the sketch labeled (a): one version of the sketch illustrates models that concentrated on the B-Z dichotomy [START_REF] Charney | Multiple flow equilibria in the atmosphere and blocking[END_REF][START_REF] Charney | Comparison of a barotropic blocking theory with observation[END_REF][START_REF] Benzi | The statistical properties of general atmospheric circulation: observational evidence and a minimal theory of bimodality[END_REF], the other on models (e.g., [START_REF] Legras | Persistent anomalies, blocking and variations in atmospheric predictability[END_REF] that allowed for the presence of additional clusters, like those found by Kimoto and Ghil (1993a) or [START_REF] Smyth | Multiple regimes in Northern Hemisphere height fields via mixture model clustering[END_REF], viz. opposite phases of the NAO and PNA anomalies (P N A, RN A and BN AO N AO -in sketch (a) of Fig. 11). The LFV dynamics in this approach is given by the preferred transition paths between the two or more regimes; see again Table 1 in Sec. 2.1 and references therein.

A third approach is associated with the idea of oscillatory instabilities of one or more of the multiple fixed points that can play the role of regime centroids. Thus, [START_REF] Legras | Persistent anomalies, blocking and variations in atmospheric predictability[END_REF] found a 40-day oscillation arising by Hopf bifurcation off their blocked regime B, as illustrated in sketch (b) of the figure. An ambiguity arises, though, between this point of view and a complementary possibility, namely that the regimes are just slow phases of such an oscillation, caused itself by the interaction of the mid-latitude jet with topography, cf. Sec. 3.2. Thus, Kimoto and Ghil (1993b) found, in their observational data, closed paths within a Markov chain whose states resemble well-known phases of an intraseasonal oscillation. Such a possibility was confirmed in the QG3 model by [START_REF] Kondrashov | Weather regimes and preferred transition paths in a three-level quasigeostrophic model[END_REF]. Furthermore, multiple regimes and intraseasonal oscillations can coexist in a two-layer model on the sphere within the scenario of "chaotic itinerancy" (Itoh andKimoto, 1996, 1997).

Finally, sketch (d) in the figure refers to the role of stochastic processes in S2S variability and prediction, whether it be noise that is white in time -as in [START_REF] Hasselmann | Stochastic climate models. I: Theory[END_REF] or in LIMs [START_REF] Penland | Random forcing and forecasting using principal oscillation pattern analysis[END_REF][START_REF] Penland | A stochastic model of IndoPacific sea surface temperature anomalies[END_REF][START_REF] Penland | Forecasting Northern Hemisphere 700-mb geopotential height anomalies using empirical normal modes[END_REF][START_REF] Penland | The optimal growth of tropical sea surface temperature anomalies[END_REF] -or red in time, as in EMRs and MSMs [START_REF] Kravtsov | Multi-level regression modeling of nonlinear processes: Derivation and applications to climatic variability[END_REF][START_REF] Kravtsov | Empirical model reduction and the modeling hierarchy in climate dynamics and the geosciences[END_REF][START_REF] Kondrashov | Empirical mode reduction in a model of extratropical low-frequency variability[END_REF][START_REF] Kondrashov | Low-order stochastic model and "past-noise forecasting" of the Madden-Julian oscillation[END_REF]Kondrashov et al., , 2015) ) or even non-Gaussian [START_REF] Sardeshmukh | Understanding the distinctively skewed and heavy tailed character of atmospheric and oceanic probability distributions[END_REF]. Stochastic processes may enter into models situated on various rungs of the modeling hierarchy, from the simplest conceptual models to high-resolution GCMs. In the latter, they may enter via stochastic parametrizations of subgrid-scale processes (e.g., Palmer and Williams, 2009, and references therein), while in the former they may enter via stochastic forcing, whether additive or multiplicative, Gaussian or not (e.g., Kondrashov et al., 2015, and references therein).

Figure 11 simply summarizes some of the key dynamical mechanisms of mid-latitude S2S variability, as discussed in this chapter, without providing a definitive answer to which approach to low-order modeling and prediction is or will be most productive in the near future. Such an answer is likely to be long in coming, if ever. On the more practical aspects of S2S prediction, though, it is clear that the ideas presented herein -on the role of LOMs vs. more detailed and highly resolved models, and on extratropical sources of variability and hence of predictability -are certain to be of use. The most promising LOMs appear to be those that do include nonlinear dynamics, memory effects, and colored noise. EMRs and MSMs that do so have proven quite competitive in ENSO forecasting, according to [START_REF] Barnston | Skill of real-time seasonal ENSO model predictions during 2002-2011 -is our capability improving[END_REF], and are likely to provide important prediction benchmarks for S2S forecasting, as well as important tools for a more complete understanding of the dynamics of extratropical S2S variability and teleconnections.

We conclude with a few comments about the role of climate change in the S2S variability and predictability problem. As mentioned in Sec. 1, the difficulty of the S2S problem is largely due to its being mid-way between the short-term weather problem and the long-term climate problem [START_REF] Von Neumann | Some remarks on the problem of forecasting climatic fluctuations[END_REF]. A mathematically self-consistent way out of this difficulty is provided by the theory of non-autonomous and random dynamical systems, as introduced recently into the climate sciences (e.g., [START_REF] Ghil | Climate dynamics and fluid mechanics: Natural variability and related uncertainties[END_REF]Chekroun et al., 2011a;[START_REF] Drótos | Probabilistic concepts in a changing climate: A snapshot attractor picture[END_REF]Ghil, 2017, and references therein). This theory -as opposed to the classical theory of autonomous dynamical systems, in which coefficients and forcing do not depend on time -allows one to include very fast forcing, by stochastic processes, as well as much slower forcing. Examples of the latter include not only changes in radiative forcing due to anthropogenic changes in aerosol and greenhouse gas concentrations, but also forcing by ocean-atmosphere (e,g., Ghil, 2001, 2017, andreferences therein) or stratosphere-troposphere interactions (e,g., Holton et al., 1995, and references therein). There is no room here to either provide a more detailed account of the appropriate concepts and methods or to actually apply it to NH LFV. Fairly readable accounts for the climate scientist do exist, though [START_REF] Ghil | Climate variability: Nonlinear and random aspects[END_REF][START_REF] Moron | Weather types across the Maritime Continent: From the diurnal cycle to interannual variations[END_REF]Chang et al., 2015, Ch. 2), and might encourage the interested reader to use them in the S2S context.

Figure 1 :

 1 Figure 1: Multiple equilibria of a three-mode quasi-geostrophic model with simplified forcing and topography. (a) Bifurcation diagram showing model response to changes in forcing; see text for the explanation of abscissa and ordinate. The "S"-shaped bifurcation curve is typical of two back-to-back saddle-node bifurcations that give rise to two stable solution branches (solid) separated by an unstable one (dashed). (b) Flow patterns of the zonal (upper) and blocked (lower) branch, corresponding to the two stable equilibria Z and R (after Charney and DeVore, 1979). (Reprinted with permission from Ghil and Childress (1987): Copyright 1987, Springer Nature.)

Figure 2 :Figure 3 :

 23 Figure 2: Schematic Markov chain with three regimes, B, C, and Z, for 'Blocked', 'Wave Train' and 'Zonal'. (a) Some preferential paths between pairs of regimes are shown, along with the corresponding transition probabilities, for instance pZB . (b and c) The distribution of residence times in log-linear coordinates differs from the straight line associated with a red-noise process; the mean residence time for each regime is denoted by τ . (Reprinted with permission from[START_REF] Ghil | Dynamics, statistics and predictability of planetary flow regimes[END_REF]): Copyright 1987, Springer Nature.)

Figure 4 :

 4 Figure 4: Spatial structure of intraseasonal Z500 variability with (a) periods of 50 days, and (b) of 43 days; based on the M-SSA analysis in Fig. 3. Shown is the relative standard deviation (in %) at each grid point, i.e. relative to the total standard deviation of overall Z500 variability.

Figure 5 :

 5 Figure 5: Phase composites of the 43-day oscillation from the M-SSA analysis of Z500 in Fig. 3. Four phases that correspond to one-half of the cycle are shown in color (in m). The phase categories are labeled in the panel legends by their midpoint (in days); also shown are the extreme values (min and max, in m).

Figure 6 :

 6 Figure 6: Same as Fig. 5, but showing the phase composites of the 50-day oscillation.

Figure 7 :

 7 Figure 7: Gaussian mixture-model centroids, showing streamfunction anomaly maps at 500 hPa, for the QG3 model: (a) NAO + ; (b) NAO -; (c) AO + ; and (d) AO -. Positive contours are solid and land masses are shaded; twenty contour levels between maximum and minimum values are used, with the following intervals (in 10 6 m 2 s -1 ): (a) 1.1; (b) 0.8; (c) 0.8; and (d) 1.1. After[START_REF] Kondrashov | Empirical mode reduction in a model of extratropical low-frequency variability[END_REF]. ( c American Meteorological Society. Reproduced with permission.)

Figure 8 :

 8 Figure 8: Mixture-model clusters for EMR model simulations of the full QG3 atmospheric model. The four panels differ by the amplitude of the EMR's stochastic forcing, which is multiplied by a factor : a) = 0.2; b) = 0.4; c) = 0.6; d) = 1(optimal case). Note that the number of the clusters and their separation in phase space increase with . The four regimes are only labeled in panel (d). After[START_REF] Kondrashov | Empirical mode reduction in a model of extratropical low-frequency variability[END_REF] ( c American Meteorological Society. Reproduced with permission.)

Figure 10 :

 10 Figure 10: Prediction skill for the optimal EMR model obtained solely from the x 2 (t) time series. The best cubic two-level EMR model clearly outperforms LIM both in terms of anomaly correlations (red curve vs. magenta) and root-mean-square (RMS) errors (blue vs. cyan).

Figure 11 :

 11 Figure 11: Schematic overview of atmospheric LFV mechanisms.

Table 1 :

 1 Classification methods for weather maps. Based onGhil and Robertson (2002) and updated.

	Approach	Method	Data sets	References	Comments
			Regime classification by position
	Cluster analysis	Categorical	NH	Mo and Ghil (1988)	Fuzzy
			NH + sectorial	Michelangeli et al. (1995)	Hard (k-means)
			Model	Dawson and Palmer (2014)	Hard (k-means)
			Model+NH+sectorial

  provide further dynamical understanding Time series of the component x 2 (red curve) and its LFV mode obtained by SSA reconstruction (blue curve). The latter has a dominant period of 25 nondimensional units, and it captures 36% of the total variance.

			LFV of X 2
	0.32
						data
						SSA-RC
	0.3
	0.28
	0.26
	0.24
	0.22
	0.2
	0.18
	0.16
	1400	1500	1600	1700	1800	1900
				time
	Figure 9:

Acknowledgments. It is a pleasure to acknowledge detailed and constructive comments by Fréderic Vitart and an anonymous reviewer. This work was partially supported by U.S. Office of Naval Research under MURI grants N00014-12-1-0911 and N00014-16-1-2073.