A Dynamical Systems Characterisation of Atmospheric Jet Regimes

. Atmospheric jet streams are typically separated into primarily "eddy-driven", or "polar-front", jets and primarily "thermally-driven", or "subtropical" jets. Some regions also display "merged" jets, resulting from the (quasi) co-location of the regions of eddy generation with the subtropical jet. The different locations and driving mechanisms of these jets issue from very different underlying mechanisms, and result in very different jet characteristics. Here, we link the current understanding of the dynamical jet maintenance mechanisms, mostly issuing from conceptual or idealised models, to the phenomena observed 5 in reanalysis data. We speciﬁcally focus on developing a unitary analysis framework, grounded in dynamical systems theory, which may be applied to both idealised models and reanalysis, and allow for direct intercomparison. Our results illustrate the use of dynamical systems indicators to diagnose jet regimes


Introduction
An idealised model fit for the task at hand should be able to maintain three distinct jet structures: a thermally-driven subtropical jet, an eddy-driven polar-front jet and a merged jet. Just as in the real atmosphere, these regimes should differ in the location and variability of the jet stream, as well as in the structure, zonal wavenumber and phase speed of the dominant modes. 55 Here, we use the two-layer modified quasi-geostrophic (QG) spherical model of Lachmy and Harnik (2014). Unlike other QG models, our setup includes advection of the zonal-mean momentum by the ageostrophic mean meridional circulation. This enables the model to resolve the momentum balance of the subtropical jet. The model can therefore reproduce the three different jet regimes and their distinct wave-mean flow feedback mechanisms (Lachmy and Harnik, 2016).
To elucidate the intrinsic dynamical characteristics of the different regimes, we interpret atmospheric flows as representative 60 of the evolution of chaotic atmospheric attractors. Recent advances in dynamical systems theory have demonstrated that any instantaneous state of a chaotic system may be described by two metrics: the local dimension -related to the system's active degrees of freedom around that particular state -and a local measure of persistence (Lucarini et al., 2016;Faranda et al., 2017b). This approach can easily be applied to a variety of datasets, including suitably processed reanalysis data, and may thus be used to provide a direct analogy between modelled and observed flows. Here, we apply it for the first time to the study of atmospheric jets in both idealised model and reanalysis data. In doing so, we underscore the links that can be made between the dynamical systems metrics and the physical characteristics of the atmospheric flow.
The connection between the dynamical systems characteristics of the flow and the jet regimes issues from the wave spectrum and its interaction with the zonal-mean flow. According to Lachmy and Harnik (2016), the flow in the idealized two-layer modified QG model transitions from a subtropical jet regime to a merged jet regime and then to an eddy-driven jet regime as 70 the eddy energy is increased. At the transition between the merged and eddy-driven jet regimes, the wave energy spectrum becomes turbulent and an inverse energy cascade takes place. At this transition, an increase in the jet's latitudinal variability and a decrease in its characteristic variability time scale are also seen (Lachmy and Harnik, 2020). We expect these changes to be reflected in the dynamical systems metrics, which capture the active degrees of freedom and the persistence of the flow.
We first provide a brief overview of the data, model and analysis approaches in Section 2. Sections 3 and 4 outline the 75 dynamical characteristics of the different flow regimes in the model and reanalysis, respectively. Finally, we discuss these results in the context of both idealised models and studies of the observed atmospheric jet, and draw our conclusions in Section 5. This study uses the numerical model of Lachmy and Harnik (2014), designed as a minimal complexity representation of jet dynamics. The model is based on the QG approximation for a sphere and includes the interactions between the zonal-mean zonal wind, the mean meridional circulation and the eddies, quantified as deviations from zonal-mean values. To balance the heat and momentum budgets, radiative damping to an equilibrium profile and surface friction are also included. The model has two vertical layers, which are meant to represent the lower and upper troposphere. Here, we refer to these with the subscripts l and u, respectively. Two key aspects of the model are: a numerical hyperdiffusion scheme which dissipates energy away from the smallest scales and a representation of the advection of zonal-mean momentum by the mean meridional circulation through an ageostrophic term. The waves are treated separately from the zonal-mean flow to allow for a clear conceptual separation between them. This modified QG framework allows the study of the complex interactions between the mean flow, the waves, and the mean meridional circulation, while retaining the simplicity of an idealised model. Since the QG assumptions are not 90 valid close to the equator, the model dynamics are most relevant to the extratropics. The crude vertical resolution leads to sharp regime transitions, which would likely be smoother in a system with a higher resolution. Despite its limitations, the model nonetheless captures the qualitative characteristics of the observed jet regimes (Lachmy and Harnik, 2016). For a detailed description of the model equations, we refer the reader to Appendix A.
The model setup, radiative equilibrium profile, and fixed parameter values used here are the same as in Lachmy and Harnik 95 (2016) and are provided in Appendix A. These are meant to mimic wintertime conditions. We focus our analysis, both in the model and in the reanalysis data (Sect. 2.2), on the SH because it is closer to zonal symmetry than its northern counterpart.
In the model, the different flow regimes are obtained through different combinations of the layer thickness (H) and the wavedamping (r) parameters. r specifically represents the ratio between the damping parameters for the eddies and for the zonalmean flow. As H and r are increased, the flow becomes more stable and the equilibrated eddy amplitude decreases. We consider 100 a parameter sweep of 27 different combinations of 7 km ≤ H ≤ 10 km and 0.5 ≤ r ≤ 2. The increments of H and r are 0.5 km and 0.5, respectively. We removed the simulation with H = 8 km and r = 2 from the analysis, because it was dominated by unrealistically regular oscillations of the eddy amplitude. The numbering of these runs is shown in Table A1. The parameter values were chosen so that all the observed jet regimes are captured, while the eddies are not completely stabilized.
In the model, we diagnose jet characteristics using: barotropic zonal-mean zonal wind (0.5 × (U l + U u ), where overbars 105 denote zonal means) and barotropic wave vorticity (0.5 × (q l + q u )). The two give complementary information on the flow, by considering the zonal-mean and the wave fields, respectively. We choose to analyze the barotropic (i.e., vertical mean) variables, since they capture most of the variability of both the upper and lower layer components. The analysis was repeated using the 3-D (upper and lower levels) full potential vorticity and wave potential vorticity fields. The results support the conclusions drawn from barotropic zonal-mean zonal wind and barotropic wave vorticity concerning the relative differences in d and θ 110 between the different jet regimes (not shown).

Reanalysis Data
Part of the analysis is conducted on data from the European Centre for Medium-Range Weather Forecasts' ERA-Interim reanalysis (Dee et al., 2011). We use daily average winds interpolated at a horizontal resolution of 0.5 • over the austral winters (June, July and August, JJA) of 1979-2017. We focus on a South Pacific domain spanning 120 • W -120 • E, 15 • S -75 • S 115 (see Fig. 1). This is chosen based on Bals-Elsholz et al. (2001); Nakamura and Shimpo (2004) and Koch et al. (2006) to focus on a longitudinal region with coherent jet characteristics and displaying a range of jet regimes (see also Sect. 4). The jet is diagnosed using the 300 hPa and 850 hPa zonal wind. These are widely used isobaric surfaces to study upper and lowertropospheric flows (e.g. Meleshko et al., 2016), and we select them here in analogy to the upper and lower-level flows in the QG model. The dynamical systems metrics (Sect. 2.3 below) are then computed on barotropic zonal wind, defined as the average of the zonal wind at these two levels. We further calculate eddy kinetic energy (EKE) from daily data. To separate the eddy flow field from the mean component, we apply a six-day high-pass Lanczos filter with 61 weights. The use of EKE allows an effective comparison with model results, as this quantity provides both a clear separation between model jet regimes (Fig. 2a) and can be used to characterise SH jet variability (e.g. Inatsu and Hoskins (2004); Shiogama et al. (2004)).

125
The above data is analysed by applying a recently developed approach grounded in dynamical systems theory. This evolved from the findings of Lucarini et al. (2016) and Faranda et al. (2017b), and allows computing the instantaneous (in time) or local (in phase-space) properties of a dynamical system by combining extreme value theory with Poincaré recurrences. A given succession of latitude-longitude maps of an atmospheric variable of interest is interpreted as a long trajectory in a reduced phase-space of the atmospheric flow. Each map corresponds to both a specific point in this phase-space and a specific time.

130
Instantaneity in time is therefore equivalent to locality in phase space. Local (instantaneous) properties are then computed for all points (timesteps) in our dataset. We specifically compute two metrics, namely the local dimension d and the persistence The local dimension is the local counterpart to the attractor of a dynamical system, namely a geometrical object defined in the phase space hosting all the possible states of the system. It describes the geometry of the system's trajectory in a small region 135 of the phase space around a state of interest ζ, which in our case could be a latitude-longitude map of barotropic zonal wind from reanalysis data. d may be taken as a proxy for the number of active degrees of freedom of the system about ζ. In other words, it measures how many directions all possible trajectories originating from a certain state in phase-space can take locally on the attractor. The persistence θ −1 of a state ζ is a measure of the system's typical residence time in the neighbourhood of ζ. That is, a measure of how long the system persists in states that closely resemble ζ. Both d and θ may be related to a state's 140 intrinsic predictability (Messori et al., 2017;Hochman et al., 2019a). Indeed, the attractor dimension is related to error growth and thus predictability. If trajectories are constrained on a low-dimensional manifold, they will be more easily predictable than if they evolve on a high-dimensional phase space, with many different directions available (Buschow and Friederichs, 2018). This argument may be extended to d when considering a specific point on the attractor. Thus, a state with a low d will have a higher intrinsic predictability than one with a high d. Similarly, one may interpret highly persistent states (high θ −1 ) as having 145 a high intrinsic predictability. Indeed, a high θ −1 implies that a persistence forecast -the simplest possible forecast one can make -will roughly capture the short-term evolution of the state. Forecasting the evolution of highly transient states (low θ −1 ) will typically require an understanding of the dynamics of the system, and thus correspond to a lower intrinsic predictability.
However, since the degree of predictability is jointly determined by d and θ, specific dynamical regimes such as travelling linear waves may be comparatively predictable and display a low d, yet also a low persistence. The fact that both d and θ −1 are local 150 metrics, implies that the intrinsic predictability is conceptually different from the predictability inferred from the performance of a numerical weather prediction model. Although the information provided by the two partly overlaps (Scher and Messori, 2018), the former depends on the local geometry of the attractor and on the characteristic time scales of the dynamics in the neighborhood of the state of interest. On the other hand, the numerical forecasts depend on the specific numerical model used and, for forecasts with a long lead-time, on non-local properties of the trajectory.

155
To estimate the local dimension we leverage the Freitas-Freitas-Todd theorem (Freitas et al., 2010), modified by Lucarini et al. (2012), which characterises the system's recurrences around the state of interest ζ. The theorem specifically indicates that the cumulative distribution function of suitably defined recurrences of the system about ζ converges to the exponential member of the Generalised Pareto Distribution (GPD). As detailed in Appendix B, the local dimension may then be estimated directly from the parameters of the GPD. Here, we follow the estimation procedure of (Faranda et al., 2019b).

160
The persistence is instead obtained from the extremal index, which we can interpret as the inverse of the mean cluster size of recurrences about ζ (Moloney et al., 2019). Further details on the estimation of the extremal index and the derivation of θ are provided in (Süveges, 2007) and Appendix B here. The persistence is bounded in [1, ∞] and is in units of the time step of the data. It is therefore essential to use a dataset whose time step is smaller than the typical timescale of the physical processes of interest. An overly long time step would indeed result in all instantaneous states tending to θ −1 = 1. In line with previous 165 studies, we deem daily data sufficient to capture the salient features of the large-scale jet variability (e.g. Woollings et al., 2010;Madonna et al., 2017).
As final product of our analysis, one obtains a value of d and θ −1 for every time step and variable in the dataset. Estimating d and θ −1 for real-world data is subject to a number of caveats, as discussed in Appendix B. There are nonetheless both formal and empirical arguments supporting the application of our framework to data that deviate from the theoretical case, and indeed 170 the two dynamical systems metrics have been successfully applied to a variety of climate datasets (e.g., Rodrigues et al., 2018;Scher and Messori, 2019;Hochman et al., 2019b;Faranda et al., 2019a, b;Brunetti et al., 2019;De Luca et al., 2020b, a), and more generally to a range of chaotic dynamical systems Pons et al., 2020). The simple model we adopt can reproduce the three jet regimes: subtropical, merged and eddy-driven. We classify our simulations according to the structure and driving mechanism of the zonal-mean zonal wind and according to the variability and spectral properties of the flow, following Lachmy and Harnik (2016). The subtropical jet is located at the edge of the Hadley cell, displays weak eddy kinetic energy and is maintained by zonal-mean advection of planetary momentum. The merged jet is located inside the Ferrel cell and is maintained mainly by eddy momentum flux convergence. It has a narrow latitudinal 180 structure and a very low latitudinal variability. The eddy-driven jet is also maintained by eddy momentum flux convergence, but is much wider than the merged jet and displays large fluctuations between a single and double jet structure. The characterization of the different flow regimes is further detailed in Appendix A. As the wave energy is increased by decreasing H and/or r (Sect. 2.1), the flow transitions from a subtropical jet regime to a merged jet regime and finally to an eddy-driven jet regime. In Fig. 2a, the different simulations are marked according to the flow regime and sorted by their time mean EKE. The 185 different spectral properties of the flow in the three regimes are apparent in Fig 2b. In the merged jet regime the wave spec-trum is dominated by a wavenumber 5 mode, whereas in the eddy-driven jet regime the spectrum is much wider and extends to lower wavenumbers. The model further reproduces flows which are a mixture of the eddy-driven and merged jet regimes, meaning that the simulations either vacillate between the two regimes or that the flow displays characteristics of both regimes, depending on the chosen diagnostic variable (crosses in Fig. 2).

Dynamical characteristics of the model jet regimes
The differences in time-mean EKE and EKE spectra of the three jet regimes (plus the mixed case) reproduced by the model are related to their different spatio-temporal variabilities (see also Sect. 3.1, and Lachmy and Harnik, 2020). This, in turn, suggests that their dynamical systems characteristics should be different. Wave driving is the dominant mechanism affecting the jet and its variability on a wide range of time scales. A strong EKE and a wide EKE spectrum should favour nonlinearities 195 in the flow and allow for a rich set of possible evolutions. The converse holds for weak EKE and a narrow EKE spectrum. We would thus expect the jet stream to be more persistent (smaller θ) and display a lower local dimension d when the EKE is weak and the spectrum is narrower (see also Sect. 5).
We analyse the results for d and θ computed on both barotropic zonal-mean zonal wind and barotropic wave vorticity ( (turquoise dots), to the mid-range EKE merged jet (amber dots), to the strong EKE eddy-driven jet regimes (black dots), via the mixed cases (blue dots). Although the spread in d and θ within each simulation is large, the centroids of the different simulations are organized in specific regions in the d-θ phase-space, according to the corresponding jet regime (Fig. 3d).
The picture from the barotropic wave vorticity is more nuanced (Fig. 3a, c). The subtropical jet regime (turquoise dots) typically displays low d and θ, matching the expected weak wave activity and the dominance of thermal maintenance mech-205 anisms. The eddy-driven jet regime (black dots) shows higher d and θ, reflecting the dominant role of eddies and the large meridional excursions in jet location. The merged jet regime (amber dots) and mixed cases (blue dots) display on average a lower d yet a higher θ than the eddy-driven jet regime, although the mixed and eddy-driven cases show a substantial overlap in both indicators. A low local dimension coupled with a low persistence is a somewhat unusual combination, both in terms of our a priori expectations on the dynamical characteristics of the different jet regimes and in terms of the positive correlation 210 between d and θ displayed by most atmospheric variables (e.g., Faranda et al., 2017b, a;Messori et al., 2017). This result may reflect the narrow wave spectrum driving the merged jet regime, leading to a quasi-periodic behaviour in which wave trains with a dominant zonal wavenumber develop and circle the globe with a relatively high phase speed (due to the relatively strong jet) before decaying. The wave field evolution is thus a highly predictable eastward translation with a low windspeed variability and a weak meridional meandering of the jet, as reflected in the low local dimension. At the same time, the wave field displays 215 a low persistence due to its relatively fast phase propagation. Consistently, the quasi-linear nature of this regime, in which the waves weaken the jet in place while they decay, results in d and θ values of the barotropic zonal-mean zonal flow which lie in-between those of the linear subtropical jet regime and of the turbulent eddy-driven jet regime (Fig. 3b, d).
A more detailed picture of the evolution of the flow's dynamical characteristics can be obtained by ranking the different simulations by decreasing EKE, which closely mirrors the division between the different jet regimes (Fig. 4). The barotropic 220 wave vorticity displays an increase in θ at the transition between the eddy driven and merged jet regimes, and a large discontinuous decrease in θ at the transition between the merged and subtropical jet regimes around simulation 23 (Fig. 4a). The latter simulation is discussed in further detail below. As discussed above, we interpret the large values of θ in the merged jet regime, which indicate low persistence, as a result of the regime's narrow spectrum with a single dominant propagating wave mode.
The barotropic wave vorticity d shows, on the other hand, a decrease as a function of EKE, albeit with significant variabil-225 ity within the individual jet regimes. The highest values of d appear in the most energetic simulations of the eddy driven-jet regime, consistent with their turbulent nature as seen in the wave spectrum (Fig. 2b) and also discussed in Lachmy and Harnik (2016). Both dynamical metrics for the barotropic zonal-mean zonal wind (Fig. 4b, d) display a largely monotonic decrease with increasing simulation index. Again, the transitions between different jet regimes, and especially that between the merged and subtropical jet regimes (around simulation 23), emerge as discontinuities. Fig. 4 also highlights the range of variability 230 of the dynamical systems metrics within each simulation. Some distributions, such as the local dimension of the subtropical jet regime (Fig. 4c), are very peaked, indicating homogeneous characteristics of the flow throughout the simulations. Others, such as the eddy-driven regime in the same panel, display a broader peak and a wider range of variability, pointing to the rich dynamical structure of the flow.
We conclude the overview of the QG model's jet regimes by focusing on simulation no. 23, which serves well to illustrate 235 the sensitivity of the dynamical systems metrics to the instantaneous characteristics of the atmospheric flow. This simulation falls within the subtropical jet regime (Table A1 and Fig. 2a), but displays some anomalous characteristics during a period of approximately 1200 days in the middle of the simulation (Fig. 5c). These are clearly visible in Fig. 4, but also emerge in Fig.   3a as a cloud of turquoise points extending towards high θ values. During this period, the flow is characterised by vacillations between a high-EKE, southerly merged jet regime and a subtropical jet. Unlike the mixed jet simulations, however, this is an 240 isolated occurrence, explaining why the simulation falls into the subtropical jet classification. The dynamical systems metrics enable us to identify clearly this anomalous transition period (red dots in Fig. 5a, b), which covers a separate region of the d-θ space compared to the rest of the simulation. Indeed, the outcropping region in Fig. 3a can be almost exclusively ascribed to the anomalous period (cf. Figs. 3a, 5a). The persistence of this period for the barotropic wave vorticity is intermediate between the merged jet and subtropical jet clusters, while the local dimension is close to that of the subtropical jet cluster (turquoise 245 cross in Fig. 3c). The barotropic zonal-mean zonal wind instead displays anomalously low local dimension and persistence, with the former being much lower than even those of the linear subtropical jet regime (turquoise crosse in Fig. 3d). Thus, in simulation no. 23 the dynamical systems metrics reflect the transition of the jet to a state with peculiar dynamical properties.

Dynamical characteristics of jet regimes in reanalysis data
As discussed in the introduction, pure thermally-driven or eddy-driven jets are largely theoretical constructs. Moreover, for a 250 given region and season, the jet typically displays large intraseasonal variability and a range of different flow characteristics (e.g., Bals-Elsholz et al., 2001;Nakamura and Shimpo, 2004;Woollings et al., 2010;Harnik et al., 2014;Messori and Caballero, 2015;Madonna et al., 2017). On these bases, we expect the jets observed in the real atmosphere to blend the dynamical characteristics of the different jet regimes reproduced in the QG model.
To verify whether our dynamical systems approach can distinguish between different jet regimes in the real atmosphere, 255 we therefore take a slightly different angle from that used in our idealised model, although the analysis tools are identical to those used above. We compute d and θ on barotropic zonal wind (see Sect. 2.2), but consider flow anomalies associated with concurrent low or high values of d and θ. We define these as values beyond the 10 th or 90 th percentiles of the respective distributions. The choice of joint high or low percentiles is dictated by the fact that the different jet regimes in Fig. 3b, d align chiefly along a d-θ diagonal. The analysis focuses on the JJA months over the Pacific domain shown in Fig. 1, chosen because   260 it displays a variety of jet regimes. During the Austral winter, the signature of the subtropical -and predominantly thermallydriven -jet is evident at upper levels around 30 • S, while that of the polar -and predominantly eddy-driven -jet emerges most clearly at lower levels around 60 • S (Fig. 1). Periods when these flows are well-separated alternate with periods dominated by a single near-barotropic flow, akin to a merged jet.
When both d and θ are anomalously low, the upper-level zonal wind (300 hPa) displays a clear maximum around 30 • S, 265 coincident with the climatological location of the subtropical jet and characterised by above-average speeds (Fig. 6a). The signature of this jet extends to the lower troposphere, where a low-level (850 hPa) zonal wind maximum located slightly poleward of the 300 hPa maximum, is clearly visible on the same days, especially in the eastern part of the domain (Fig. 6c).
This inference is confirmed by analysing the vertical cross-section of the zonal flow, which evidences a poleward displacement of the location of the zonal wind maximum at 850 hPa compared to 300 hPa (Fig. 6e). The EKE shows widespread negative 270 anomalies across the central and southern portions of the domain, and localised positive anomalies at 300 hPa matching the latitude of the jet (Fig. 7a, c). We interpret the above characteristics as the signature of a merged jet.
When both d and θ are anomalously high, the upper-level zonal flow evidences two maxima: a weaker-than-climatology jet around 30 • S and an anomalous secondary maximum at around 55 • S, associated with large positive zonal windspeed anomalies ( Fig. 6b). At low levels, there is a single jet located just to the North of 60 • S and associated with large positive zonal flow 275 anomalies (Fig. 6d). The southern jet therefore has a more pronounced barotropic structure than its northern counterpart, as also seen in the vertical cross-section of the zonal flow (Fig. 6f). The EKE anomalies are predominantly positive across the central and southern portions of the domain (Fig. 7b, d). The above points to a primarily eddy-driven jet in the South and a primarily thermally-driven jet further North, in agreement with the theoretical framework of Lee and Kim (2003) and Son and Lee (2005). It also suggests that the persistence of the double-jet configuration, with a strong eddy-driven jet, is lower (high θ) 280 than that of the single merged jet (low θ).

Discussion and conclusions
We have analysed different jet regimes in a set of idealised simulations with a QG model and reanalysis data. The QG model reproduces the full range of theoretically-derived jet regimes (e.g., Lee and Kim, 2003;Son and Lee, 2005;Lachmy and Harnik, 2016). These are: an eddy-driven jet, a merged jet and a subtropical jet. It additionally reproduces transition states between 285 a merged and eddy-driven jet, which we term mixed jets. Pure eddy-driven or thermally-driven jets are largely theoretical constructs, and in the real atmosphere the separation between them is often blurred. Nonetheless, it is possible to identify primarily eddy-driven, primarily thermally-driven and merged jets (e.g., Koch et al., 2006;Eichelberger and Hartmann, 2007;Li and Wettstein, 2012). The South Pacific sector is of particular interest in this context, since it is a region where the three main jet regimes may be observed at the same longitude (e.g., Bals-Elsholz et al., 2001;Nakamura and Shimpo, 2004). 290 Relating the different jet regimes identified in idealised models to those identified in reanalysis data is far from immediate.
Here, we have proposed an analysis approach which may be applied to both datasets, and which provides a direct link between the characteristics of the jets in the QG model and those seen in the ERA-Interim dataset. Such approach is grounded in dynamical systems theory and is based on two metrics, d and θ, which characterise the instantaneous (local in phase-space) dynamical characteristics of the jet. The local dimension d is a proxy for the number of active degrees of freedom of the system.

295
The persistence θ −1 is a measure of the typical residence time of the system in the neighbourhood of a given state. Their computation issues from an analysis of recurrences of the system, and both may be related to the concept of predictability.
Indeed, one may expect low d, high θ −1 states to be more predictable than high d, low θ −1 situations (see Sect.

2.3). This interpretation is confirmed by ongoing work which finds a direct link between precipitation forecast skills at single stations in
France and the values of d and θ computed for 500 hPa geopotential height fields in the preceding days. Unlike conventional 300 approaches which diagnose the driving mechanisms of the jet in terms of complex physical processes, such as convection or eddy momentum flux convergence, the dynamical systems metrics allow to delineate the dynamical properties of the system uniquely from computing d and θ on the wind itself. Differently from previous studies using this approach, here we have thus attempted to link the dynamical systems metrics directly to the physical mechanisms of jet maintenance. As a caveat, we note that the absolute values of these metrics should be interpreted in a relative sense within each dataset they are computed on, 305 and that direct comparison of their magnitudes across datasets should be treated with caution. From a theoretical standpoint, there are cogent arguments supporting the use of recurrences of observables of a system to investigate the properties of the system's underlying phase-space, even though all the variables defining said phase-space cannot be considered (e.g., Faranda et al., 2017a, c;Barros et al., 2019). This makes our analysis affordable in terms of computational and data requirementsindeed d and θ for a single variable from a reanalysis dataset may be computed by a modern laptop computer in a matter of 310 hours -and versatile in terms of applicability to very different datasets (cf. Rodrigues et al., 2018;Faranda et al., 2019c;Pons et al., 2020). Moreover, d and θ do not require prior knowledge of the exact location of the jet. In the QG model, d and θ allow to characterise the different jet regimes. Specifically, when computed on barotropic zonalmean zonal wind they highlight the eddy-driven jet as being a low-persistence, high-dimensional regime, the subtropical jet as being a high-persistence, low-dimensional regime and the merged jet as having intermediate dynamical characteristics. This 315 reflects the wide wave spectrum and large latitudinal fluctuations in the eddy-driven jet, versus the lower latitudinal variability of the subtropical jet. Indeed, in the former regime the energetic wave activity and the comparatively broad wave spectrum and range of wave phase speeds offers a variety of different temporal evolutions of the atmospheric flow. In the latter regime, the flow instead displays a weak eddy kinetic energy and a narrow wave spectrum, with a corresponding narrow range of wave phase speeds and meridional excursions. Moreover, the typically similar wave length and phase speed of the waves implies that 320 the temporal evolution of the atmospheric flow is relatively repetitive. Harnik (2016, 2020) further analyzed the EKE budgets of different jet regimes in terms of the contributions by linear versus eddy-eddy interactions, and found that the degree of linearity of the dynamics changes. The subtropical regime is comparatively linear, the merged regime's variability is dominated by wave-mean flow interactions, and the eddy-driven regime is nonlinear, with a turbulent upscale energy cascade.
Intuitively, we expect the predictability of the flow to decrease, and the jet structure to vary more strongly in time, as the flow 325 becomes more nonlinear. The dynamical systems metrics thus reflect the theoretical expectations on the role of increased EKE in decreasing predictability (Leith, 1971). In contrast with the jet variability, the model's jet speed (Fig. A1) is not related to the degree of nonlinearity of the flow. Because of this, changes in jet speed may not map directly onto differences in d and θ.
A natural question to ask is whether d and θ also reflect variability within the individual simulations in the QG model.
A preliminary analysis on the barotropic zonal-mean zonal wind for one of the simulations classified as eddy-driven regime 330 (simulation no. 2, see Table A1) shows that, when d and θ are very low, the jet has a large variability (red line in Fig. A2c) and a narrow profile reminiscent of the merged jet regime (cf. the red line in Fig. A2b with the amber line in Fig. A1). This may be contrasted with the average jet structure for this simulation which is quite wide, with an inflection suggestive of a separation between the eddy-driven and subtropical jets (gray line, Fig. A2b). During low d and θ timesteps, the jet also displays a much sharper EKE spectrum (not shown), and a more equatorward lower level jet (red line, Fig. A2d) than climatology. On the other 335 hand, when d and θ are anomalously large, the lower-level jet is slightly more poleward and intense than usual (cf. the grey and blue lines in Fig. A2b,d), and the spread in jet profiles is unusually low (blue line in Fig. A2c). The resemblance of low d and θ timesteps in an eddy-driven jet simulation to the flow characteristics seen in merged jet regime simulations is consistent with inter-regime differences in the two dynamical systems indicators. Indeed, simulations classified as merged jet regimes generally display lower d and θ values than simulations within the eddy-driven regime (Fig. 3b). Lachmy and Harnik (2020) 340 also found that the temporal variability of the eddy-driven jet also spans a merged jet-like state.
The results for reanalysis data largely mirror those found for the idealised simulations. Focusing on the Pacific sector of the Southern Ocean, d and θ computed on the barotropic zonal wind again discriminate between different jet regimes. Specifically, high d, high θ days display both a thermally-driven, subtropical jet and a strong eddy-driven, polar-front jet, while low d, low θ days display a merged jet. A key difference from the results in the QG model is that, in the reanalysis data, the south Pacific 345 sector we consider does not display a purely thermally-driven nor a purely eddy-driven jet. Indeed, both jet regimes co-occur at different latitudes during winter. The upper-level jet is stronger in the single-jet merged configuration compared to the double jet configuration (cf. Fig. 6a and b), consistent with the idealized results of Son and Lee (2005).
Our intention is not to provide a systematic analysis of atmospheric jet characteristics in different geographical regions.
Rather, by building upon both idealised simulations and reanalysis data, we illustrate the potential of dynamical systems 350 indicators to diagnose jet regimes in a conceptually intuitive fashion.

Appendix A: Description of the 2-layer QG model
The 2-layer QG model, described in section 2.1, solves the equations outline below.
The barotropic (vertical mean) and baroclinic (half the difference between the upper and lower layers) components of the zonal-mean zonal momentum equation: and the diagnostic equation for (V a ) T : where U ≡ u · cos φ, V ≡ v · cos φ and µ ≡ sin φ. The variables u and v are the zonal and meridional winds, respectively, 360 and φ is the latitude. Subscripts M and T denote the barotropic and baroclinic components, respectively. Overbars denote zonal means and primes denote deviations from the zonal mean. The constants Ω, a, τ f , τ r and ν are Earth's rotation rate and radius, the surface friction time scale, the radiative damping time scale and the numerical diffusion coefficient, respectively.
The non-dimensional parameter is defined as 8 aΩ N H 2 , where N is the Brunt-Väisälä frequency.
The barotropic and baroclinic components of the wave potential vorticity (PV) equation are: where λ is the longitude in radians. q and ψ are the PV and stream function, respectively, which satisfy the following relation: The components of the mean PV gradient, ∂q M ∂µ and ∂q T ∂µ , which appear in Eqs. (A3a) and (A3b) are: The model includes several parameters which were fixed for all the experiments in this study. The surface friction and radiative relaxation time scales are: τ f = 3.9 days and τ r = 11.6 days, respectively. The radiative equilibrium profile of the thermal wind, (u T ) E (φ), is the following profile smoothed by a running average filter: where ( the wave spectrum (Lachmy and Harnik, 2016). The diagnosed regimes for each of the simulations are listed in Table A1.
The different structures of the zonal-mean zonal wind are shown in Fig. A1 for three simulations, one from each regime. The upper layer zonal wind (thick solid lines) represents the upper tropospheric jet stream. The lower layer zonal wind (thin dashed lines) indicates the structure of the mean meridional circulation, since according to the zonal-mean momentum balance, it is positive in the Ferrel cell and negative in the Hadley and polar cells (Lachmy and Harnik, 2014). We identify the mechanism 390 maintaining the jet according to the relative location of the upper and lower layer zonal wind maxima. The subtropical jet (turquoise lines) is thermally driven, since its upper-layer zonal wind maximum is at the Hadley cell edge, where the lower layer zonal wind is zero. The merged jet (amber lines) and eddy-driven jet (black lines) are located inside the Ferrel cell, where the lower layer zonal wind is maximal, indicating that they are driven by eddy momentum flux convergence. In the merged jet regime the jet inside the Ferrel cell is collocated with the maximum vertical shear of the zonal wind, which indicates that 395 it represents a merging of the subtropical and eddy-driven jets. In the eddy-driven jet regime the two maxima are separated, indicating that the upper-layer zonal wind maximum represents a purely eddy-driven jet (Lachmy and Harnik, 2016).
The attractor of a dynamical system is a geometric object defined in the space hosting all the possible states of the system (phase-space). For Each point ζ on the attractor, two dynamical indicators can be computed to characterise the recurrence 400 properties of ζ. The local dimension d indicates the number of degrees of freedom active locally around ζ; the persistence θ −1 is a measure of the mean residence time of the system around ζ (Faranda et al., 2017b). To determine d, a key observation is the connection between extreme value theory and the Poincaré recurrences of a state ζ in a chaotic dynamical system. We consider long trajectories of a dynamical system -in our case successions of daily atmospheric latitude-longitude maps -as a sequence of states on the attractor. This equates to approximating the attractor with its Sinai-Ruelle-Bowen (SRB) measure 405 (Young, 2002), and establishes a connection between the attractor and the histogram of visits in phase-space. For a given point ζ in phase-space (e.g., a given atmospheric latitude-longitude map), one may then compute the probability that the system returns within a ball of radius centered on ζ, which follows an extreme value distribution. Indeed, the Freitas-Freitas-Todd theorem (Freitas et al., 2010), modified by Lucarini et al. (2012), states that logarithmic returns: yield a probability distribution such that: This is the exponential member of the Generalized Pareto Distribution family, with z = g(x(t)) and s a high threshold associated to a quantile q of the series g(x(t)). Requiring that the orbit falls within a ball of radius around the point ζ is equivalent to asking that the series g(x(t)) exceeds the threshold s; therefore, the ball radius is simply e −s(q) . The parameters 415 µ and σ, namely the location and the scale parameters of the distribution, depend on the point ζ in phase space. µ(ζ) corresponds to the threshold s(q), while the local dimension d(ζ) is obtained as: σ = 1/d(ζ).
When x(t) contains all the variables of the system, the estimation of d based on extreme value theory has a number of advantages over traditional methods (e.g. the box counting algorithm (Liebovitch and Toth, 1989;Sarkar and Chaudhuri, 1994)). First, it does not require to estimate the volume of different sets in scale-space: the selection of s(q) based on the 420 quantile provides a selection of different scales s which depend on the recurrence rate around the point ζ. Moreover, it does not require the a priori selection of the maximum embedding dimension, as the observable g is always a univariate time-series.
The persistence of the state ζ is measured via the extremal index 0 ≤ ϑ(ζ) ≤ 1, an adimensional parameter which measures the inverse of cluster lengths for consecutive exceedances of s or, in other terms, recurrences. From this we extract θ(ζ) = ϑ(ζ)/∆t, where ∆t is the timestep of the data. θ(ζ) is therefore the inverse of the average residence time of trajectories around 425 ζ and is in units of frequency (in this study days −1 ). If ζ is a fixed point of the attractor θ(ζ) = 0. For a trajectory that leaves the neighborhood of ζ at the next time iteration, θ = 1. To estimate ϑ, we adopt the Süveges maximum likelihood estimator (Süveges, 2007). This approximates the value of ϑ for processes with asymptotically independent inter-cluster times, with the threshold s and the length of the timeseries as asymptotes. The estimator performs a first guess of the value of ϑ by using the exponential nature of the inter-cluster times. Then, a second estimate is given by extending the definition of clusters to 430 include sequences of exceedances interrupted by short "gaps" which are not distributed exponentially, and may therefore be viewed as single clusters. The first estimate provides a upper bound on ϑ; the second estimate, a lower bound. The algorithm is then iterated until it reaches a stable point. Heuristically, clusters should be much shorter than inter-cluster times. The other assumptions underlying this calculation largely match those made to compute d, discussed below. For its ease of computation, we prefer the Süveges estimator to the exact formulas provided by Caby et al. (2019), that are anyways conceived for periodic The above framework, and in particular the Freitas-Freitas-Todd theorem (Freitas et al., 2010), holds for stationary, Axiom A systems in the limit of infinite timeseries. In this respect climate data, and specifically reanalysis data, poses a number of challenges. The dynamics of the system are non Axiom A and non-stationary, which means that different areas of the attractor  This points to the applicability of the extreme value framework to real-world data. Empirical tests support this conclusion, and indeed Buschow and Friederichs (2018) have shown that d successfully reflects the dynamical characteristics of the atmosphere even when applied to datasets that deviate from the theoretical case, and where universal convergence to the exponential member of the GPD is not achieved.
Code availability. The code to compute the dynamical systems indicators used in this study is made freely available through the cloud   Table A1). Units are m 2 s −2 for both panels. The markers on both panels indicate the jet regime of each simulation. Open circles, black circles, gray circles and ×'s indicate the different jet regimes as in the legend. The larger markers indicate the simulations shown in Fig. A1.  the different jet regimes. Note that higher (lower) θ values corresponds to lower (higher) persistence, and that the ordinates' ranges differ across panels. (see Table A1).     Figure A1. zonal-mean zonal wind (m s −1 ) in the upper (uu, thick solid lines) and lower (u l , thin dashed lines) layers of the 2-layer QG model, for simulation number 1 from the eddy-driven jet regime, simulation number 20 from the merged jet regime and simulation number 27 from the subtropical jet regime. Colours indicate the different jet regimes as in the legend. See Table A1 for the parameters of the simulations. Figure A2. d-θ diagram of barotropic zonal-mean zonal wind (a) for simulation no. 2 (see Table A1) and composites over timesteps with low (< 10 th ) and high (> 90 th ) percentiles of both d and θ, of upper (b) and lower (d) level zonal mean zonal wind (m s −1 ) and (c) barotropic zonal mean zonal wind standard deviation. The vertical bars in (b) and (d) show the standard deviations of the relevant composites at a given latitude. The red and blue dots (a) and lines and bars (b-d) correspond to low and high d-θ timesteps, respectively, while gray marks all other times (a) or the simulation's climatology (b-d). The numbers in parentheses in (a) show the sample size of the each composite. Thick lines in (b), (d) denote latitudes at which the high or low d and θ composites are significantly different from the simulation's climatology at the 95% confidence level, based on a two-sided student t-test.