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 midlatitude jet dynamics by embedding atmospheric data and by investigating its properties (bifurcation structure, stability, local dimensions) for different atmospheric flow regimes. The derivation is a three-step process: first, we obtain a 1-D description of the midlatitude jet stream by computing the position of the jet at each longitude using ERA-Interim. 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 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 modified topography. Our model, notwithstanding its simplicity, provides an instructive description of the dynamical properties of the atmospheric jet.

1 Introduction

Jet streams are narrow, fast-flowing westerly air currents near the tropopause. They are a major feature of the large-scale atmospheric circulation and modulate the frequency, severity and persistence of weather events across the extratropics (e.g., Röthlisberger et al., 2016). Their location and intensity also affects commercial aviation and shipping (Reiter and Nania, 1964; Hadlock and Kreitzberg, 1988; Williams and Joshi, 2013). Two types of atmospheric jets can be identified: thermally driven subtropical jets, and eddy-driven jets associated with baroclinic instability at the polar front. In the Northern Hemisphere (NH), the two are not always clearly separated (Lee and Kim, 2003), and when considering monthly or longer time averages, a single, 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 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 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 tested 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., 2016a). Moreover, jet dynamics have been described as a manifestation of multiple equilibria in asymmetrically forced 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 three-dimensional 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 and Ghil, 2005), climate tipping points (Stommel, 1961; Benzi et al., 1982), large-scale atmospheric motions (Lorenz, 1984, 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 behaviors (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 techniques to derive low-dimensional models from experimental data (Letellier et al., 2006). This opposition was motivated by a long sequence 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 data embedding. Recently, Faranda et al. (2017a) 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 (sub-grid) dynamics are represented as stochastic perturbations. Here, we use these results to develop a minimal model of the effective dynamics of the midlatitude 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 et al. (2017a) 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 a given 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 or breaking and the response to topographical features, and we relate the results back to the original reanalysis data.

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

2 Data and methods

2.1 ERA-Interim data and jet position algorithm

The analysis is based on the European Centre for Medium Range Weather Forecasts's ERA-Interim (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 d low-pass Lanczos filter (Duchon, 1979). We then identify the 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 grid points, except at longitude 0. The analysis has been repeated for BRI thresholds between 5 and 15^∘, with no significant qualitative differences.

Figure 1 shows a snapshot of the jet position on 4 February 1979, together with the time series of the daily jet position recorded in 1979 at longitude 120^∘ W. An animation of the jet location for the year 1980 is provided as a supplementary video. Both the time series and the snapshot show large jumps in the jet position. A qualitative analysis of the jet position data suggests that the jet fluctuates around a central latitude (Central Jet, CJ) and seldom shifts to more northerly (NJ) or southerly (SJ) latitudes.

Figure 1 Snapshot of the jet position extracted from the ERA-Interim data set on 4 February 1979 and time series of the jet position for the year 1979, recorded at longitude 120^∘ W.

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 the BRI corresponding to about 10^∘ latitude is $\left|x\right|>\mathrm{1}$.

2.2 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 θ⁻¹. Instantaneity in time corresponds to locality in phase space, such that a value of d and θ⁻¹ can be computed for a given variable (in our case the jet position data) at every time step. 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 how it can evolve from such state. θ⁻¹ describes the persistence of a state in time, thus providing complementary information to d.

2.2.1 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 centered 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:

$$\begin{array}{}\text{(1)}& g\left(x\right(t\left)\right)=-\log \left(\mathrm{dist}\right(x\left(t\right),\mathit{\zeta }\left)\right),\end{array}$$

such that close recurrences of ζ correspond to large values of g(x(t)) (Collet and Eckmann, 2009). Thus, the probability 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 of g(x(t)) in Eq. (1) locks this probability into the exponential member of the generalized Pareto distribution:

$$\begin{array}{}\text{(2)}& \mathrm{Pr}(z>s(q\left)\right)\simeq \exp \left[-\mathit{\vartheta }\left(\mathit{\zeta }\right)\left({\displaystyle \frac{z-\mathit{\mu }\left(\mathit{\zeta }\right)}{\mathit{\beta }\left(\mathit{\zeta }\right)}}\right)\right],\end{array}$$

where z=g(x(t)) and μ and β (obtained via fitting) depend on the point ζ. These are the location and the scale parameters of the distribution. Remarkably, $\mathit{\beta }\left(\mathit{\zeta }\right)=\mathrm{1}/d\left(\mathit{\zeta }\right)$, 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., 2017c, b; 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.

2.2.2 Local persistence

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., 2016b): an adimensional parameter $\mathrm{0}<\mathit{\vartheta }\left(\mathit{\zeta }\right)<\mathrm{1}$ which can be interpreted as the inverse of the mean residence time. We can then compute ${\mathit{\theta }}^{-\mathrm{1}}\left(\mathit{\zeta }\right)=\mathrm{d}t/\mathit{\vartheta }\left(\mathit{\zeta }\right)$, where dt is the time step of our data. Heuristically, if the ith visit to the neighborhood of ζ lasts τ_i consecutive time steps and N such visits are made in total, then ${\mathit{\theta }}^{-\mathrm{1}}\approx (\mathrm{1}/N){\sum }_i{\mathit{\tau }}_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 neighborhood of ζ. θ=1 corresponds to a trajectory residing in the neighborhood 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. (2017c) 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

3.1 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 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, and (iv) smaller-scale phenomena, such as turbulence and baroclinic waves, will 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-dimensionalized 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 return map $x_t^{\left(i\right)}$ vs. $x_{t+\mathrm{1}}^{\left(i\right)}$ (an example is shown in Fig. 2) and searching for a function f⁽ⁱ⁾ such that $x_{t+\mathrm{1}}^{\left(i\right)}=f^{\left(i\right)}\left(x_t^{\left(i\right)}\right)$. The first thing to verify is that the same functional form f⁽ⁱ⁾ may be used at all longitudes i. This is equivalent to asking that there is only one dynamic driving the jet independently of the location. With the choice

$$\begin{array}{}\text{(3)}& f^{\left(i\right)}\left(x\right)=\left\{\begin{array}{ll}-\frac{A(A+x)}{A-c}+r^{\left(i\right)},& x<-c,\\ \sinh \left(\mathit{\beta }x\right)+r^{\left(i\right)},& -c\le x\le c,\\ \frac{A(A-x)}{A-c}+r^{\left(i\right)},& c<x,\end{array}\right.\end{array}$$

where we have dropped the dependencies of x for clarity, the parameters can be fixed at all longitudes as β=0.75, A=3, and $c={\sinh }^{-\mathrm{1}}\left(A\right)/\mathit{\beta }\approx \mathrm{2.4246}$. Even though the functional form of f⁽ⁱ⁾ is independent of longitude, a dependence on i remains in the form of the parameter r⁽ⁱ⁾, which represents the effects of topography in terms of spatial inhomogeneities of the local dynamics. As a first-order approximation, we consider only the difference between land and ocean and assign one of two discrete values to each r⁽ⁱ⁾. The choice of the function f⁽ⁱ⁾ 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):

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{x}_{t+\mathrm{1}}^{\left(i\right)}=& (\mathrm{1}-\mathit{ϵ})f^{\left(i\right)}\left(x_t^{\left(i\right)}\right)+\mathit{ϵ}{f}^{(i-\mathrm{1})}\left(x_t^{(i-\mathrm{1})}\right),\\ & (i=\mathrm{1},\mathrm{2},\mathrm{\dots },N;t=\mathrm{1},\mathrm{2},\mathrm{\dots }).\end{array}\end{array}$$

With this geometry, the dynamics are divided into N=360 cells. Periodic boundary conditions are applied at N=360. The dynamics in each cell i are time-independent but perturbed by the cell i−1 (i.e., its neighbor to the west) with intensity ϵ, which 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 midlatitudes, namely of the order of 100 km.

Figure 2 The average return map extracted from the data at longitude i=1. This is constructed by coarse-graining the state space at i=1 into M partitions $L_k^{\left(\mathrm{1}\right)}$ (k=1, 2, …, M). We then define ${\stackrel{\mathrm{‾}}{x}}^{(\mathrm{1},k)}$ as the midpoint of the partition $L_k^{\left(\mathrm{1}\right)}$, and $Y^{(\mathrm{1},k)}=\mathit{\left\{}{x}_t^{\left(\mathrm{1}\right)}\right|x_{t-\mathrm{1}}^{\left(\mathrm{1}\right)}\in L_k^{\left(\mathrm{1}\right)}\mathit{\}}$ (t=2, 3, …). The black dots represent (${\stackrel{\mathrm{‾}}{x}}^{(\mathrm{1},k)}$, 〈Y^(1,k)〉) for k=1, 2, …, 500, where 〈⋅〉 is the average over the elements of Y^(1,k), computed based on the observed data. The red line represents the approximated average return map $\langle Y^{(\mathrm{1},k)}\rangle =f^{\left(\mathrm{1}\right)}\left({\stackrel{\mathrm{‾}}{x}}^{(\mathrm{1},k)}\right)$ when $\left|x\right|\le c$. In the region $\left|x\right|>c$, we assume linear reflection effects. As a result, we have the return map f⁽¹⁾ in Eq. (3).

3.2 Sub-grid feedbacks to jet dynamics

If we perform a numerical simulation of Eq. (4), the dynamics are fixed to one of the three states (CJ, SJ, NJ), depending on the value of ϵ. This means that the role of small-scale perturbations in triggering the transitions between the states is fundamental. We therefore have to include an additive noise term ${\mathit{\xi }}_t^{\left(i\right)}$ in Eq. (4):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{x}_{t+\mathrm{1}}^{\left(i\right)}=& (\mathrm{1}-\mathit{ϵ})f^{\left(i\right)}\left(x_t^{\left(i\right)}\right)+\mathit{ϵ}{f}^{(i-\mathrm{1})}\left(x_t^{(i-\mathrm{1})}\right)+{\mathit{\xi }}_t^{\left(i\right)},\\ & (i=\mathrm{1},\mathrm{2},\mathrm{\dots },N;t=\mathrm{1},\mathrm{2},\mathrm{\dots }).\end{array}\end{array}$$

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 affect the jet dynamics, such as convection or the interaction between the jet stream and gravity waves (Williams et al., 2003, 2005). Translated to our model with a reference spatial scale of the order of 100 km these phenomena, ranging from a few meters to a few kilometers, imply a perturbation in the range ${\mathrm{10}}^{-\mathrm{4}}<\mathit{\nu }<{\mathrm{10}}^{-\mathrm{3}}$. Several sub-grid 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 describe turbulent eddies (McComb, 1992; Frederiksen and Davies, 1997). Following these ideas, we model ${\mathit{\nu }}_t^{\left(i\right)}\in [-\mathit{\delta }$, δ] as a uniform random variable.

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 are also driven by the effect of an inverse cascade transferring energy to large scales via baroclinic waves (Held and Larichev, 1996).

Baroclinic activity is associated with extratropical cyclones and anticyclones, on scales of the order of 10³ km. These can affect the jet position by several degrees of latitude. Again, there is no 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 of noise η⁽ⁱ⁾, acting at intermediate scales (i.e., between the scale of the jet and the scale of turbulence), is enough to obtain reliable jet breaking dynamics (see Sect. 4.1). The simplest choice for ${\mathit{\eta }}_t^{\left(i\right)}\in [-\mathit{\mu }$, μ] is a block noise taking the same value over bl cells, where bl is the one-dimensionalized size of cyclones or anticyclones (see Appendix B) and obeying the uniform distribution. Another choice for modeling 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 from the embedding procedure used to derive the dynamics of x.

The minimal sub-grid parametrization can thus be written in the following form:

$$\begin{array}{}\text{(6)}& {\mathit{\xi }}_t^{\left(i\right)}={\mathit{\nu }}_t^{\left(i\right)}+{\mathit{\eta }}_t^{\left(i\right)},\end{array}$$

where ${\mathit{\nu }}_t^{\left(i\right)}$ and ${\mathit{\eta }}_t^{\left(i\right)}$ model the effects of small turbulent disturbances and baroclinic eddies, respectively (Fig. 3). The noise terms are discussed further in Appendix B.

Figure 3 Schematic representation of noise contributions to the CML model (Eqs. 3 and 6): ν represents local turbulent disturbances, r topographical features, η baroclinic eddies and i longitudinal positions.

Figure 4 Bifurcation diagrams as a function of κ for (a) land ($r^{\left(i\right)}=-\mathrm{0.02}$) and (b) ocean (r⁽ⁱ⁾=0.0). The gray regions delimit the accessible region of the dynamics with respect to all possible external forcings. A realization of the dynamics is depicted by the black dots. For $r^{\left(i\right)}=-\mathrm{0.02}$ and $\mathrm{0.1574}<\mathit{\kappa }<\mathrm{0.1985}$, there is a small chance of reaching SJ positions and no chance of reaching NJ positions.

3.3 Local dynamics

Owing to the unidirectional coupling in our model and to the large N, the local dynamics can be approximated by a nonautonomous dynamical system:

$$\begin{array}{}\text{(7)}& x_{t+\mathrm{1}}^{\left(i\right)}\simeq f^{\left(i\right)}\left(x_t^{\left(i\right)}\right)+p_t^{\left(i\right)},\end{array}$$

where a nonautonomous external force $p_t^{\left(i\right)}$ is given by

$$\begin{array}{}\text{(8)}& p_t^{\left(i\right)}=\mathit{ϵ}\left[f^{(i-\mathrm{1})}\left(x_t^{(i-\mathrm{1})}\right)-f^{\left(i\right)}\left(x_t^{\left(i\right)}\right)\right]+{\mathit{\nu }}_t^{\left(i\right)}+{\mathit{\eta }}_t^{\left(i\right)}.\end{array}$$

When $\left|p_t^{\left(i\right)}\right|\to \mathrm{0}$, the local dynamics has a stable fixed point at x≈0 and two unstable chaotic sets near $x\approx \pm \mathrm{2}$. When $p_t^{\left(i\right)}>\mathrm{0}$, the resulting perturbed dynamics may exhibit escape behavior 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^{\left(i\right)}=-\mathrm{0.02}$) and ocean (r⁽ⁱ⁾=0; see Appendix B) obtained by approximating the external perturbation $p_t^{\left(i\right)}$ as a random variable obeying the uniform distribution in [−κ, κ]. They both indicate a bifurcation to chaotic and partially chaotic behavior (Sato et al., 2018). The different value of r⁽ⁱ⁾ over land gives rise to an asymmetry in the invariant sets, namely the sets delimiting the accessible region of the dynamics with respect to all possible external perturbations $p_t^{\left(i\right)}$. In Fig. 4, these dynamically reachable regions are depicted in gray, while a realization of the dynamics is depicted by the black dots. With $r^{\left(i\right)}=-\mathrm{0.02}$ over land and $\mathrm{0.1574}<\mathit{\kappa }<\mathrm{0.1985}$, there is a small chance of reaching SJ positions and no chance of reaching NJ positions. This is reflected in the skewed distribution of $x_t^{\left(i\right)}$. For the sake of conciseness, we do not report the detailed bifurcation analysis of the local dynamics here. A brief analysis of the global dynamics is presented in Sect. 5.

4 Validation of the model against ERA-Interim data

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 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 μ=0.6, bl = 15, ϵ=0.33, κ=0.5, and $\mathit{\delta }=\mathrm{0.5}\times {\mathrm{10}}^{-\mathrm{4}}$. We further compare the results of model runs containing all noise terms with runs where individual terms are suppressed: coupling (ϵ=0), topography (r=0) and baroclinic waves (bl = 1).

4.1 Spatiotemporal dynamics

We first consider the latitudinal distribution of the yearly median jet positions at each longitude (dots in Fig. 5) and their interannual mean (solid lines in Fig. 5). The ERA-Interim data set (Fig. 5a) presents a negative interannual mean jet position at almost all longitudes, with noticeable zonal asymmetries and a marked interannual variability. The best model run (Fig. 5b) captures both the interannual variability and, thanks to the term r, the longitudinal variations in average location. A run without coupling (ϵ=0) is shown in Fig. 5c). In this simulation the dynamics are local, except for the presence of block noise, resulting in a discontinuous jet profile. Unlike the ERA-Interim data, the run with no topography (Fig. 5d) has median values which are roughly symmetric around zero. Finally, the run with suppressed baroclinic activity (Fig. 5e) has a smaller interannual variability than the ERA-Interim data and sharp changes in the median values of x following the geographic constraints.

Figure 5 Single-year median location (dotted points) and multiyear average of medians (solid lines) of the meridional jet position for ERA-Interim (a) and model (b–e) data. (b) Best-fit model, obtained with μ=0.6, bl = 15, ϵ=0.33 and $\mathit{\delta }=\mathrm{0.5}\times {\mathrm{10}}^{-\mathrm{4}}$; (c) as in (b) but with ϵ=0; (d) as in (b) but with r=0; (e) as in (b) but with bl = 1. The simulations consist of 37 years of daily jet positions.

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 $\left|x\right|>\mathrm{1}$. Note that this corresponds to breaks in the jet position with the same threshold defined by the BRI, although there is not exact correspondence. We then assign “0” to all the observations with $\left|x\right|<\mathrm{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 dynamical phases in complex systems (Kaneko, 1990). The so-obtained binary spatiotemporal dynamics are shown in Fig. 6a–e for all the previously described runs. In the ERA-Interim data (Fig. 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 a westerly propagation of the clusters. The best model fit captures the qualitative aspects of this behavior, although the longitudinal coherence is weaker (see Sect. 4.2 below). In the remaining model simulations, the suppression of different noise terms alters either the cluster size or the westerly propagation of the clusters (Fig. 6c–e). A quantitative analysis of the cluster size spectra is presented in Fig. 6f for space clusters and Fig. 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 are turbulent, with energy at all scales. Despite its simplicity, our model 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 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 consists of discrete peaks with the energy concentrated at precise 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.

Figure 6 (a–e) Space–time daily representation of the binarized jet dynamics: 1 (black) corresponds to a NJ or SJ shift ($\left|x\right(t\left)\right|>\mathrm{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 Fig. 5. Space (f) and time (g) cluster spectra for the binarized ERA-interim data (black) and the different model runs (colors).

4.2 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 (Fig. 7a), θ (Fig. 7b) and the BRI (Fig. 7c) for every day in the data set. 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., 2017b). 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 midlatitude atmospheric circulation are presented by Madonna et al. (2017). The model runs can produce average dimensions comparable to those observed in the ERA-Interim data, except for the bl = 1 case. There, the fragmented dynamics lead to a much higher dimension. This is consistent with the spatiotemporal diagrams shown in Fig. 6. The models' inverse persistence θ is slightly larger than that observed in reanalysis but still of the order of 2 d. 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. (2017c), who found that high d matches blocking-like atmospheric configurations in the North Atlantic region. For the limiting bl = 1 case, the BRI is also correlated with θ: the more breaks, the lower the persistence of the flow.

Figure 7 Box plots of the local dimension d (a), inverse persistence θ (b) and breaking index BRI (c) for the ERA-Interim data and four numerical simulations as in Fig. 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.

5 Bifurcation diagram and jet regimes

The bifurcation diagram in Fig. 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 look symmetric with respect to x=0 if r=0 everywhere, but the addition of topography via r⁽ⁱ⁾ alters the relative proportion 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 behavior in the stochastic bifurcation obtained from the approximated local dynamics (Fig. 4). By analyzing the bifurcation structure of the conceptual model as a function of the coupling coefficient – which mimics the coherence of the jet – we identify three behaviors: (i) a strong and uniform jet 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 ($\mathrm{0.6}<\mathit{ϵ}<\mathrm{0.9}$); and (iii) an intermediate state of transition between the two. These jet regimes are broadly consistent with those obtained in idealized atmospheric simulations (Lachmy and Harnik, 2016; Son and Lee, 2005), although here we do not delve into the physical mechanisms underlying the different behaviors. 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.

6 Conclusions

We have derived a minimal model of the jet stream position dynamics, based on a stochastic coupled map lattice, by embedding data extracted from ERA-Interim. This procedure innovates over earlier studies (e.g., Faranda et al., 2017a) 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 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 jet's 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 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 BRI defined here.

Figure 8 Bifurcation diagram of the global dynamics obtained for μ=0.6, bl = 15, $\mathrm{0}<\mathit{ϵ}<\mathrm{1}$ and $\mathit{\delta }=\mathrm{0.5}\times {\mathrm{10}}^{-\mathrm{4}}$. The diagram represents the density of states ρ(x) obtained by varying ϵ. The vertical gray line indicates the value used as best fit to the ERA-Interim data.

The analysis we have conducted can, however, already answer some of the questions left open in Faranda et al. (2017b) and Madonna et al. (2017) concerning the possibility of reducing the complex midlatitude circulation dynamics to low-dimensional 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 investigations observed relatively low dimensions when considering the full sea-level pressure fields (Faranda et al., 2017c, b). It also suggests that breaks in the jet are responsible for higher dimensions.

Appendix A: Coupled map lattice

A coupled map lattice (CML; Kaneko, 1983) is given by

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{x}_{t+\mathrm{1}}^{\left(i\right)}=& (\mathrm{1}-\mathit{ϵ})f\left(x_t^{\left(i\right)}\right)+{\displaystyle \frac{\mathit{ϵ}}{\mathrm{2}}}\left[f\left(x_t^{(i-\mathrm{1})}\right)+f\left(x_t^{(i+\mathrm{1})}\right)\right],\\ & (i=\mathrm{1},\mathrm{2},\mathrm{\dots },N;t=\mathrm{1},\mathrm{2},\mathrm{\dots }),\end{array}\end{array}$$

where ϵ∈[0, 1], $x_t^{\left(i\right)}\in \mathbb{R}$ and $f\left(x\right):\mathbb{R}⟶\mathbb{R}$. For our jet dynamics, we adopt the open flow model, which is a class of CML with unidirectional coupling (Kaneko, 1985):

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{x}_{t+\mathrm{1}}^{\left(i\right)}=& (\mathrm{1}-\mathit{ϵ})f\left(x_t^{\left(i\right)}\right)+\mathit{ϵ}f\left(x_t^{(i-\mathrm{1})}\right),\\ & (i=\mathrm{1},\mathrm{2},\mathrm{\dots },N;t=\mathrm{1},\mathrm{2},\mathrm{\dots }).\end{array}\end{array}$$

The CML is a phenomenological model to study complex spatiotemporal 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.

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_k^{\left(i\right)}$ (k=1, 2, …, M) and let ${\stackrel{\mathrm{‾}}{x}}^{(i,k)}$ be the midpoint of $L_k^{\left(i\right)}$. Then, we construct a set $Y^{(i,k)}=\mathit{\left\{}{x}_t^{\left(i\right)}\right|x_{t-\mathrm{1}}^{\left(i\right)}\in L_k^{\left(i\right)}\mathit{\}}$ (t=2, 3, …) and a return map f⁽ⁱ⁾ via the return plot of (${\stackrel{\mathrm{‾}}{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:

$$\begin{array}{}\text{(B1)}& \langle Y^{(i,k)}\rangle =f^{\left(i\right)}\left({\stackrel{\mathrm{‾}}{x}}^{(i,k)}\right),i=\mathrm{1},\mathrm{2},\mathrm{\dots },N,k=\mathrm{1},\mathrm{2},\mathrm{\dots },M,\end{array}$$

where $\left|x\right|\le c$. In the region $\left|x\right|>c$, we assume linear reflection effects. As a result, we have the return map f⁽ⁱ⁾ in Eq. (3). Figure 2 illustrates the construction for i=1 and M=500.

An important ingredient of the jet dynamics is the presence of topographic obstacles to the midlatitude zonal flow. Mountain ranges and land–sea boundaries cause meridional deviations in the mean jet location (Tibaldi et al., 1980). This inhomogeneity can be modeled via a parameter r⁽ⁱ⁾ that mimics this “spatial noise”. Since the topography is at most a few kilometers in height, this translates to a perturbation of the order of 10⁻³ in the model. Reasonable geographical constraints are therefore $r^{\left(i\right)}=-\mathrm{0.02}$ (i ∈ land) and r⁽ⁱ⁾=0.0 (i ∈ ocean), where “land” spans the ranges $\mathrm{0}\le i<\mathrm{161}$ and $\mathrm{239}\le i<\mathrm{301}$ and “ocean” spans the ranges $\mathrm{161}\le i<\mathrm{239}$ and $\mathrm{301}\le i<\mathrm{360}$. The negative sign for the jet shifts over land is justified by the negative median values of the ERA-Interim jet position anomalies over land (compare Fig. 5a and b with Fig. 5d where no topography is present).

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 of a few meters to a few kilometers, implying a perturbation in the range ${\mathrm{10}}^{-\mathrm{4}}<\mathit{\nu }<{\mathrm{10}}^{-\mathrm{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 (where bl is the one-dimensionalized size of cyclones or anticyclones in our model) with an amplitude of the order of 1. The latter value is determined empirically as it is an indicative magnitude of the large shifts midlatitude 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 100 km, and assuming a typical scale for extratropical 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 deviated northwards. We therefore take the block perturbation to be of the size bl = 15 blocks.

Owing to the unidirectional coupling in our lattice jet model and to the large N, the local dynamics can be approximated by a nonautonomous dynamical system $x_{t+\mathrm{1}}^{\left(i\right)}\simeq f^{\left(i\right)}\left(x_t^{\left(i\right)}\right)+p_t^{\left(i\right)}$, where the external force $p_t^{\left(i\right)}=\mathit{ϵ}\left[f^{(i-\mathrm{1})}\right(x_t^{(i-\mathrm{1})})-f^{\left(i\right)}(x_t^{\left(i\right)}\left)\right]+{\mathit{\nu }}_t^{\left(i\right)}+{\mathit{\eta }}_t^{\left(i\right)}$. Assuming that the time averages $\langle f^{(i-\mathrm{1})}\left(x_t^{(i-\mathrm{1})}\right)-f^{\left(i\right)}\left(x_t^{\left(i\right)}\right)\rangle$, $\langle {\mathit{\nu }}_t^{\left(i\right)}\rangle$ and $\langle {\mathit{\eta }}_t^{\left(i\right)}\rangle$ are all 0 by symmetry, we have $\langle p_t^{\left(i\right)}\rangle \approx \mathrm{0}$. Thus, we recover the average return map given in Eq. (B1).