Minimal dynamical systems model of the Northern Hemisphere jet stream via embedding of climate data

. We derive a minimal dynamical systems model for the northern hemisphere mid-latitude jet dynamics by embedding atmospheric data, and investigate its properties (bifurcation structure, stability, local dimensions) for different atmospheric ﬂow regimes. The derivation is a three-step process: ﬁrst, we obtain a 1-D description of the mid-latitude jet-stream by computing the position of the jet at each longitude using the ERA-Interim reanalysis. Next, we use the embedding procedure to derive a map of the local jet position dynamics. Finally, we introduce the coupling and stochastic effects deriving from both atmospheric 5 turbulence and topographic disturbances to the jet. We then analyze the dynamical properties of the model in different regimes: one that gives the closest representation of the properties extracted from real data; one featuring a stronger jet (strong coupling); one featuring a weaker jet (weak coupling); and one with modiﬁed topography. Our model, notwithstanding its simplicity, provides an instructive description of the dynamical properties of the atmospheric jet.

spiral-shaped jet structure emerges (e.g. Archer and Caldeira, 2008). In this paper we consider a single NH jet (NHJ), rather than attempting to separate the subtropical and eddy-driven jets (e.g. Belmecheri et al., 2017).
Even though the climatological NHJ is a westerly flow, it can present large meanders on synoptic timescales (e.g. Koch et al., 2006;Röthlisberger et al., 2016). These can cause the local flow to become predominantly meridional, or can even determine a splitting or breaking of the jet (Haines and Malanotte-Rizzoli, 1991). The occurrence of these large meanders in the jet is often 5 associated with events such as temperature and precipitation extremes (e.g. Dole et al., 2011;Screen and Simmonds, 2014).
Although jet dynamics are well understood in a climatological sense, our insights into dynamical features such as jet splitting or meandering are still limited.
The dynamics of meanders and split jets has often been framed in terms of transitions between zonal and blocked flows since the seminal work by Charney and DeVore (1979). Legras and Ghil (1985) and Ghil (1987) used an intermediate complexity barotropic model with dissipation forcing and topography, and observed two distinct equilibria associated with the zonal and blocked flows. Similar mechanisms have been proposed by Mo and Ghil (1988) and then performed in experimental facilities (Weeks et al., 2000). However, there is no consensus about the nature of flow multistability, and a wide range of theoretical explanations and models have been proposed (e.g. Tung and Lindzen, 1979;Simmons et al., 1983;Frederiksen, 1982;Faranda et al., 2016b). Moreover, jet dynamics have been described as a manifestation of multiple equilibria in asymmetrically forced 15 flows (Hansen, 1986) or as a result of soliton-modon structures (McWilliams et al., 1981).
In order to advance our understanding of the jet dynamics, we employ a low-dimensional dynamical systems model derived from reanalysis data. The best-known example of a low dimensional model for atmospheric phenomena is Lorenz' simple threedimensional system representing some features of Rayleigh-Bénard convection (Lorenz, 1963). Thereafter, simple dynamical systems models have been devised to study El Niño (Penland and Matrosova, 1994), ocean-atmosphere interactions (Dijkstra 20 and Ghil, 2005), climate tipping points (Stommel, 1961;Benzi et al., 1982), large scale atmospheric motions (Lorenz, 1984(Lorenz, , 1996 and many other phenomena. The goal of these investigations was not to provide the most realistic representation of the relevant systems, but rather to capture key emerging behaviours (such as chaos, intermittency, multistability). The main drawback of those investigations was the weakness of the connection between models and real data due to the scarcity of observations as well as theoretical limitations. Until very recently, there was a strong case against the use of embedding 25 techniques to derive low dimensional models from experimental data (Letellier et al., 2006). This opposition was motivated by a long sequel of papers that appeared between 1984 and 1991. The initial claim that low dimensional models for complex phenomena could be derived using a very small numbers of variables (see e.g. Nicolis and Nicolis (1984); Fraedrich (1986)) was disproved by rigorous numerical computations by Grassberger (1986) and Lorenz (1991).
Progress in data quality and availability and the advent of stochastic dynamical systems have renewed the attention for 30 data embedding. Recently, Faranda et al. (2017c) have shown that embedding techniques can yield effective low-dimensional dynamics provided that the chosen observables reflect the symmetries of the system and that small-scale (subgrid) dynamics is represented as stochastic perturbations. Here, we use these results to develop a minimal model of the effective dynamics of the mid-latitude jet. This is useful to explore a range of possible behaviors beyond those displayed in the available data, that could have appeared in past climates and could appear again in future climates. In analogy to the model derived by Faranda 35 et al. (2017c) for the von Karman turbulent flow, the jet model is based on a coupled map lattice (CML, see Appendix A).
Each element of the lattice reflects the dynamics of the jet at each longitude. Such a model does not require physical sub-grid terms a priori, but only if they are found to be essential to capture the large-scale phenomenology -which we show is not the case. We then evaluate how this model represents key dynamical features of the jet, namely its stability, the statistics of splitting/breaking and the response to topographical features, and relate the results back to the original ERA-Interim data.

5
First, we provide the details of the ERA-Interim data and of the jet detection algorithm (Section 2). We then present the stochastic coupled lattice map model and compute its bifurcation structure (Section 3). Next, we introduce some instantaneous dynamical indicators (Section 4) and use them to relate the conceptual model to more complex climate models and reanalysis data (Section 5). Finally, we highlight the open questions our results can answer and the new questions they pose (Section 6).

ERA-Interim data and jet position algorithm
The analysis is based on the European Centre for Medium Range Weather Forecasts' ERA-Interim reanalysis (Dee et al., 2011).
We consider daily data with a 1 • horizontal resolution over the period 1979-2016.
The jet position is diagnosed through a modified version of the approach by Woollings et al. (2010). We take daily mean wind-speed averaged over 200-400 hPa and apply a 10-day low-pass Lanczos filter (Duchon, 1979). We then identify the 15 latitudinal position of the jet at every longitude as the location of the strongest wind, over the band 15 • -75 • N. This approach is intended to provide a "raw" measure of the jet variability. We then consider the longitude and time dependence of the latitude of the jet to monitor its waviness.
We define an index of large jet meanders, or breaks (Breaking Index, BRI), as the daily number of meridional variations in jet position of more than 10 • of latitude across adjacent longitude gridpoints, except at longitude 0. The analysis has been 20 repeated for BRI thresholds between 5 • and 15 • , with no significant qualitative differences. In order to embed the data and derive the effective maps of the dynamics, we remove the seasonal cycle from the data by subtracting, longitude by longitude, the average meridional position for each calendar day and dividing by the standard deviation. For the deseasonalized data, the dimensionless threshold for the computation of BRI corresponding to about 10 • latitude is |x| > 1.

Local dynamical systems metrics
Our analysis leverages two recently developed dynamical systems metrics, namely: the local dimension of the attractor d and the stability or persistence of phase-space trajectories θ −1 . Instantaneity in time corresponds to locality in phase-space, such that a value of d and θ −1 can be computed for a given variable (in our case the jet position data) at every timestep. d is a proxy for the system's active number of degrees of freedom. It provides information on how the system can reach a given state, and 5 how it can evolve from such state. θ −1 describes the persistence of a state in time, thus providing complementary information to d.

Local Dimension
The local dimension is estimated by making use of extreme value statistics applied to Poincaré recurrences. The Freitas et al. (2010) theorem, modified by Lucarini et al. (2012), states that the probability of entering a hyperball with a small radius 10 centred on a state ζ on a chaotic attractor obeys a generalized Pareto distribution (Pickands III, 1975). In order to compute this probability empirically, we first calculate the series of distances dist(x(t), ζ) between the point on the attractor ζ and all other points x(t) on the trajectory. This series is transformed via the distance function: such that close recurrences of ζ correspond to large values of g(x(t)) (Collet and Eckmann, 2009). Thus, the probability 15 of entering a small hyperball around ζ is transformed into the probability of exceeding a high threshold s(q), where q is a percentile of the series g(x(t)) itself. In the limit of an infinitely long trajectory, it can be shown that the choice g(x(t)) in Eq.
(1) locks this probability into the exponential member of the generalized Pareto distribution: where where z = g(x(t)), and µ and β (obtained via fitting) depend on the point ζ. These are the location and the scale 20 parameters of the distribution. Remarkably, β(ζ) = 1/d(ζ), where d(ζ) is the local dimension around the point ζ. This result has been recently applied to a range of atmospheric and oceanic fields (e.g. Faranda et al., 2017b, a;Messori et al., 2017;Faranda et al., 2019a, b). In this paper, we use the quantile 0.975 of the series g(x(t)) to determine q. Our results are robust with respect to reasonable changes in this quantile.

25
Extreme value statistics also allow estimating the persistence of a given state ζ, by inspecting the temporal evolution of the dynamics around it. A measure of persistence around ζ can be obtained from the mean residence time of the trajectory within the neighborhood of ζ. To quantify this, we employ the so-called extremal index ϑ (Freitas et al., 2012;Faranda et al., 2016a): an adimensional parameter 0 < ϑ(ζ) < 1 which can be interpreted as the inverse of the mean residence time. We can then compute θ −1 (ζ) = dt/ϑ(ζ), where dt is the timestep of our data. Heuristically, if the ith visit to the neighbourhood of ζ lasts τ i consecutive time steps, and N such visits are made in total, then θ −1 ≈ (1/N ) ∑ i τ i . In practice, instead of this naive estimator, we compute the extremal index using the likelihood estimator of Süveges (2007). θ = 0 corresponds to a stable fixed point of the dynamics, so that the trajectory resides an infinite amount of time in the neighbourhood of ζ. θ = 1 corresponds 5 to a trajectory residing in the neighbourhood of ζ for only one time step per visit. The estimate of θ is thus sensitive to the dt used. If dt is too large, the time dependence structure is unresolved and θ will be close to 1. Conversely, if dt is too small, θ will be close to 0. Faranda et al. (2017b) observed that θ varies between 0.3 and 0.5, when dt = 1 day, for sea-level pressure fields over the North Atlantic. In this work, we use the same dt.
3 Derivation of the lattice jet model

Model Framework
While not directly issued from the Navier-Stokes equations, our framework builds on concrete physical hypotheses, namely that: (i) the physics of the jet is the same at every longitude and it is only slightly modified by the presence of topographical constraints, (ii) the jet can experience sudden breaks and shifts from its central position (CJ) to northern (NJ) or southern latitudes (SJ), (iii) the jet must propagate to the west, (iv) smaller scale phenomena, such as turbulence and baroclinic waves, will 15 be introduced in the model only if necessary to reproduce the effective dynamics in the data. This latter point is fundamentally different from the philosophy of direct numerical simulations. We construct our model starting from the local time series of the non-dimensionalised jet position x measured at each longitude i and time t. We use the simplest possible embedding procedure (see Appendix B), which consists of plotting the (an example is shown in Figure 2) and searching for a function f (i) such that x The 20 first thing to verify is that the same functional form f (i) may be used at all longitudes i. This is equivalent to asking that there is only one dynamics driving the jet independently of the location. With the choice: where we have dropped the dependecies of x for clarity, the parameters can be fixed at all longitudes as: β = 0.75, A = 3, and c = sinh −1 (A)/β ≈ 2.4246. Even though the functional form of f (i) is independent of longitude, a dependence on i remains 25 in the form of the parameter r (i) , which represents the effects of topography in terms of spatial inhomogeneities of the local dynamics. As first order approximation, we consider only the difference between land and ocean, and assign one of two discrete values to each r (i) . The choice of the function f (i) is not unique; however the one we propose here is a suitable option that satisfies hypotheses (i) and (ii) above. In order to reproduce the eastward propagation of the jet (hypothesis iii), we introduce the coupled map lattice (CML, see Kaneko (1983) and Appendix A): With this geometry, the dynamics is divided into N = 360 cells. Periodic boundary conditions are applied at N = 360. The dynamics in each cell i is time-independent, but perturbed by the cell i−1 (i.e. its neighbour to the west) with intensity ϵ, which 5 we estimate and scale based on the observed data. This further implies that our reference length-scale in the model corresponds to that of 1 • longitude in the mid latitudes, namely order 100km.

Subgrid feedbacks to jet dynamics
If we perform a numerical simulation of Eq. (4), the dynamics is fixed to one of the three states (CJ, SJ, NJ), depending on the 10 value of ϵ. This means that the role of small scales perturbations in triggering the transitions between the states is fundamental.
We therefore have to include an additive noise term ξ The noise is a fundamental ingredient for the breaking of the jet and the transition between zonal and blocked states, as shown in tank experiments and numerical simulations (Jacoby et al., 2011). Physically, noise arises from key sub-grid processes that 15 affect the jet dynamics, such as convection or the interaction between the jet stream and gravity waves (Williams et al., 2003(Williams et al., , 2005. Translated to our model with reference spatial scale of order 100km these phenomena, ranging from a few meters to a few kilometers, imply a perturbation in the range 10 −4 < ν < 10 −3 . Several subgrid parametrization of turbulence exists: the seminal works of Kraichnan (1961) and Thomson (1987) showed that if large scales are represented by a deterministic term, a single random variable can drive the turbulence term. This means that Langevin model representations are appropriate to 20 describe turbulent eddies (McComb, 1990;Frederiksen and Davies, 1997). Following these ideas, we model ν However, considering small scale turbulent disturbances to the jet dynamics is not sufficient to reproduce the blocking and breaking of the jet. Even if the introduction of ν as a stochastic term can account for the direct Kolmogorov turbulent cascade (Kolmogorov, 1941), the jet dynamics is also driven by the effect of an inverse cascade transferring energy to large scales via 25 baroclinic waves (Held and Larichev, 1996).
Baroclinic activity is associated with extra-tropical cyclones and anticyclones, on scales of order 10 3 km. These can affect the jet position by several degrees of latitude. Again, there is not a unique way to model baroclinic waves in our framework.
We follow the rationale of multiscale parametrizations as they can be theoretically justified (e.g. Wouters and Lucarini (2013); Kitsios and Frederiksen (2019)) and are numerically efficient (Faranda et al., 2014). The simple introduction of another source 30 of noise η (i) , acting at intermidiate scales (i.e. between the scale of the jet and the scale of turbulence), is enough to obtain a reliable jet breaking dynamics (see Sect. 4.1). The simplest choice for η (i) t ∈ [−µ, µ] is a block noise taking the same value over bl cells (the one-dimensionalized size of cyclones/anticyclones, see Appendix B), and obeying the uniform distribution.
Another choice for modelling baroclinic disturbances to the jet could be to introduce a second deterministic equation, weakly coupled with the jet position. However, this choice requires additional hypotheses and parameters and does not emerge naturally 5 from the embedding procedure used to derive the dynamics of x.
The minimal sub-grid parametrization can thus be written in the form: where ν

Local dynamics
Owing to the uni-directional coupling in our model and to the large N , the local dynamics can be approximated by a nonautonomous dynamical system: where a non-autonomous external force p (i) t is given by: When |p (i) t | → 0, the local dynamics has a stable fixed point at x ≈ 0, and two unstable chaotic sets near x ≈ ±2. When p (i) t > 0, the resulting perturbed dynamics may exhibit escape behaviour from the fixed point to the chaotic regions with positive Lyapunov exponents. Figure 4 shows the bifurcation diagrams as a function of κ over land (r (i) = −0.02), and ocean (r (i) = 0, see Appendix B) obtained by approximating the external perturbation p t . For the sake of conciseness, we do not report the detailed bifurcation analysis of the local dynamics here. A brief analysis for the global dynamics is presented in Sect. 5.
In this section we compare the ERA-Interim deseasonalized jet position data with numerical simulations of our model. In order to have the same statistical sample as as for the reanalysis, we simulate 37 years of daily snapshots of the jet position. The best fit of our model to the data is obtained, by a trial and error procedure, for the parameters: η = 1.2, bl = 15, ϵ = 0.33, δ = 10 −4 .
We further compare the results of model runs containing all noise terms with runs where individual terms are suppressed in 5 turn: the coupling (ϵ = 0), topography (r = 0), and baroclinic waves (bl = 1).

Spatiotemporal Dynamics
We first consider the latitudinal distribution of the yearly median jet positions at each longitude (dots in Figure 5) and their interannual mean (solid lines in Figure 5). The ERA-Interim dataset (Figure 5 We next consider the NHJ's shifts from CJ towards NJ or SJ positions. We binarize the dynamics by the detecting all the events such that |x| > 1. Note that this corresponds to breaks in the jet position with the same threshold defined by the breaking index BRI, although there is not exact correspondence. We then assign '0' to all the observations with |x| < 1 (CJ) and '1' to all the others (NJ or SJ). This procedure, known as coding, is widely used in dynamical systems analysis to identify different 20 dynamical phases in complex systems (Kaneko, 1990). The so-obtained binary spatio-temporal dynamics are shown in Figure   6a-e for all the previously described runs. In the ERA-Interim data (Figure 6a), the switch from CJ to NJ and SJ phases occurs in clusters displaying a characteristic longitudinal extent and temporal persistence. There is also some indication of westerly propagation of the clusters. The best model fit captures the qualitative aspects of this behavior, although the longitudinal coherence is weaker (see section 4.2 below). In the remaining model simulations, the suppression of different noise terms 25 alters either the cluster size or the westerly propagation of the clusters (6c-e). A quantitative analysis of the cluster size spectra is presented in Figure-6f for space clusters and Figure-6g for time clusters. There is a clear power law behavior, reminiscent of a multiscale structure (Schertzer et al., 1997). This is coherent with the claim that the underlying jet dynamics is turbulent, with energy at all scales. Despite its simplicity, our model is reproduces this power law behavior. The theoretical reasons are nontrivial and can be related to the possibility of building turbulent cascades starting from simple Langevin equations (Wouters 30 and Lucarini, 2013;Faranda et al., 2014). We underline here the necessity of adding ϵ and having bl > 1. Indeed when ϵ = 0 the spatial cluster spectrum consist of discrete peaks with the energy concentrated at precises scales. These are a resonance of the block noise size. When instead we impose bl = 1, we still recover a power-law behavior, but the slope for the temporal clustering strongly deviates from that observed in the ERA-Interim data.

Dynamical indicators
We further assess our model by means of the d and θ metrics described in Sect. 2.2, computed on both ERA-Interim data and the coupled map lattice. We also compare here the statistics of the spatial breaks of the jet, detected via the indicator BRI. Figure 7 show the box-plots of d (Figure 7-a), and θ (Figure 7-b) for each day in the data set and the breaking index BRI (Figure 7-c). The ranges of values of d and θ for the ERA-Interim data resemble closely those found for sea-level pressure fields over the Northern Hemisphere (Faranda et al., 2017a). This supports the claim that the position of the jet is indicative of large-scale features of the NH atmospheric circulation. Similar claims about the relevance of low dimensional projections in describing the mid-latitude atmospheric circulation are presented by Madonna et al. (2017). The model runs can produce 10 average dimensions comparable to those observed in the ERA-Interim data, except for the bl = 1 case. There, the fragmented dynamics leads to a much higher dimension. This is consistent with the spatio-temporal diagrams shown in Figure 6. The models' inverse persistence θ are slightly larger than those observed in reanalysis, but still of order 2 days. Here we can clearly see the effect of the noise suppression (ϵ = 0 and bl = 1) in modifying the dynamical properties by leading to lower persistence.
Finally, we remark that the number of breaks is correlated to the local dimension. This result is consistent with Faranda et al. 15 (2017b), who found that high d match blocking-like atmospheric configurations in the North Atlantic region. For the limiting bl = 1 case, BRI is also correlated with θ: the more breaks, the lower the persistence of the flow.

Bifurcation diagram and jet regimes
The bifurcation diagram in Figure 8 is constructed by plotting the empirical density ρ(x) of the jet position at all longitudes as a function of ϵ. The vertical gray line corresponds to the value of ϵ that best fits the ERA-Interim data. The diagram would 20 look symmetric with respect to x = 0 if r = 0 everywhere, but the addition of geography via r (i) alters the relative proportions of time spent in SJ versus NJ. Specifically, our asymmetric land/ocean distribution implies a southward shift of the average CJ position with increasing coupling. This is reminiscent of the behaviour in the stochastic bifurcation obtained from the approximated local dynamics (Figure 4). By analysing the bifurcation structure of the conceptual model as a function of the coupling coefficient -which mimics the coherence of the jet -we identify three behaviours: (i) a strong and uniform jet 25 where large meridional excursions in the jet location are relatively rare events (ϵ < 0.35), which is close to the jet dynamics as inferred from the ERA-Interim data; (ii) a state with sharp meridional excursions in which the jet is very unstable and on average shifted far to the South (0.6 < ϵ < 0.9); and (iii) an intermediate state of transition between the two. These jet regimes are broadly consistent with those obtained in idealised atmospheric simulations (Lachmy and Harnik, 2016;Son and Lee, 2005), although here we do not delve into the physical mechanisms underlying the different behaviours. It is also noteworthy that our model qualitatively reproduces a southern jet configuration, even though we provide it with a single NHJ and do not distinguish between eddy-driven and subtropical jets.
We have derived a minimal model of the jet stream position dynamics, based on a stochastic coupled map lattice, by embedding data extracted from the ERA-Interim reanalysis. This procedure innovates over earlier studies e.g. (Faranda et al. (2017c)) by making use of a coupled map lattice derived from a local embedding of the data, and could be adapted to systems with several degrees of freedom. Instead of embedding the data of a global observable in a high-dimensional space, we have constructed 5 the return map for the local position of the jet and then added, via coupling and noise, the physical ingredients identified in previous studies as drivers of the jet dynamics. The conceptual model is then validated and tuned using dynamical indicators of the jets dimension and persistence in the reanalysis data.
Future analyses could apply this approach to the southern hemisphere, where the role of topography is less important than in its northern counterpart. This would allow us to better constrain the influence of topography on the dimension-persistence 10 diagrams. Another possibility would be to use the low-dimensional model to build a surrogate data set of the jet positions and then apply this to atmospheric analogues, so as to construct realistic atmospheric dynamics. Finally, it would be interesting to study whether further projections of the atmospheric dynamics to a lower dimensional space are possible, beyond the model developed here, and to test possible relations between different atmospheric blocking indices and the breaking index BRI defined here.

15
The analysis we have conducted can however already answer some of the questions left open in Faranda et al. (2017a) and Madonna et al. (2017) concerning the possibility of reducing the complex mid-latitude circulation dynamics to lowdimensional representations given by blocking indices or conceptual models. The fact that the dimension-persistence diagram of our minimal model qualitatively matches many features obtained for the ERA-Interim jet position and sea-level pressure fields shows that a substantial part of the dynamics projects along a single line (the jet position). This may explain why previous 20 investigations observed relatively low dimensions when considering the full sea-level pressure fields (Faranda et al., 2017b, a).
It also suggests that breaks in the jet are responsible for higher dimensions.
The CML is a phenomenological model to study complex spatio-temporal dynamics in systems with large numbers of degrees of freedom. The idea is to discretize the dynamics in space and time, while capturing the global phenomenology of the system. CMLs have been successfully applied to processes such as turbulence in thermal convection (Yanagita and Kaneko, 1993) and turbulent puff in pipe flow (Avila and Hof, 2013). It is convenient for us to model the jet dynamics leveraging the CML approach because we can extract a local one-dimensional map from the observed time series.

5
Appendix B: Average return map and noise To extract the local jet dynamics, we construct an average return map. We first coarse-grain the state space into M partitions L (i) k (k = 1, 2, . . . , M ) and letx (i,k) be the midpoint of L k } (t = 2, 3, . . .) and a return map f (i) via the return plot of (x (i,k) , ⟨Y (i,k) ⟩), where ⟨·⟩ is the average over the elements of Y (i,k) at each longitude i, and at each partition k: 2, . . . , , N, k = 1, 2, . . . , M, (11) where |x| ≤ c. In the region |x| > c, we assume linear reflection effects. As a result, we have the return map f (i) in Eq. (3). As discussed in Sect. 3.2, noise is a fundamental ingredient in the jet dynamics. The "turbulent noise" term ν relates to physical phenomena in the range from a few meters to a few kilometers, implying a perturbation in the range 10 −4 < ν < 10 −3 , where ν is a random variable obeying the uniform distribution. The second noise term, η, relates to baroclinic activity and we model it as a block noise taking the same value over bl blocks (the one-dimensionalized size of cyclones/anticyclones in our 25 model) with an amplitude of order 1. The latter value is determined empirically as it is an indicative magnitude of the large shifts mid-latitude baroclinic systems can induce in the jet. To determine a realistic length for bl, we reason as follows: given that our model has a reference scale of about 100km, and assuming a typical scale for extra-tropical cyclones of about 3000 km, we then have that bl ≈ 30 blocks. However, the perturbations are associated with the cyclone radius rather than diameter: upstream of the cyclone, the jet will mostly be deviated southwards, while downstream of the cyclone, the jet will mostly be Owing to the uni-directional coupling in our lattice jet model and to the large N , the local dynamics can be approximated by a non-autonomous dynamical system x where the external force p t ⟩ are all 0 by symmetry, we have ⟨p (i) t ⟩ ≈ 0. Thus we recover the average return map given in Eq. (11).       . a-e) Space-time daily representation of the binarized jet dynamics: 1 (black)corresponds to a NJ or SJ shift (|x(t)| > 1) and 0 (white)corresponds to a CJ position. The results are for the ERA-Interim data (a) and model runs (b-e). The latter are the same as in Figure   5. Space (f) and time (g) cluster spectra for the binarized ERA-interim data (black) and the different model runs (colors).  Figure 5. In each box, the central mark is the median, the edges of the box are the 25th and 75th percentiles, the whiskers extend to the most extreme data points not considered outliers, and outliers are plotted individually.