Marginal stability and predator–prey behaviour within storm tracks

A predator–prey relationship between storm track intensity and growth rate is revealed in reanalysis data of the North Atlantic and North Pacific, as well as in an idealized global circulation model with a zonally asymmetric heating dipole. Averaging in the phase space of these two quantities reveals that both quantities oscillate on approximately monthly time‐scales. These oscillations occur due to quasi‐periodic bursts in storm track activity that reduce excess baroclinicity and bring the flow back towards a state that is marginally stable to those bursts. Many detailed properties of these oscillations are reproduced well by a two‐dimensional dynamical system, especially in respect of the North Atlantic storm track which is more zonally constrained than that in the North Pacific. It is predicted and observed that on average stronger storm events occur less frequently but grow on a shorter time‐scale. The results suggest that nonlinearly oscillating behaviour around a state of baroclinic neutrality is a general feature of localized storm tracks, and they offer a new perspective on the study of baroclinic instability.


Introduction
In the steady state of the Northern Hemisphere, maxima in storm track activity and baroclinicity are co-located with areas of intensified meridional temperature gradients on the eastern coasts of Eurasia and North America (Hoskins and Valdes, 1990).One may expect this co-location, as baroclinicity can be defined by the sharpness of these gradients and is proportional to the growth rate of baroclinic eddies that form the storm tracks.Conversely, when these eddies are sufficiently intense, they reduce the temperature gradients and thus the baroclinicity, an idea which has led to low-order parametrizations of eddy heat flux (which is a measure of storm track activity) as eddy diffusion (e.g.Larichev and Held, 1995;Griffies, 1998;Barry et al., 2002).This nonlinear relationship between baroclinicity and storm track activity complicates the temporal variability of these two quantities.As a result, the traditional linear theory of baroclinic instability (e.g.Charney, 1947;Eady, 1949), where a constant shear induces exponentially growing perturbations, is only valid for the initial stage of the baroclinic eddy development.
There is substantial theoretical evidence that this nonlinear relationship manifests itself as oscillations of the mean vertical wind shear (which is related to baroclinicity, by thermal wind balance) and baroclinic eddies in the governing equations.For example, using the two-layer quasi-geostrophic (QG) model with a small departure from a critical shear (beyond which instability occurs), Pedlosky (1970) studied the weakly nonlinear interaction between the mean flow and the primary wave.This nonlinearity was shown to yield an oscillation of the perturbation amplitude due to it alternately changing the sign of its vertical phase.Such phase changes are associated with reducing and increasing the overall wind shear, resulting in transitions between marginally stable and marginally unstable states, respectively.In the inviscid case, the system can oscillate reversibly around a state of baroclinic neutrality.Because this behaviour arises from the inertial properties of the interaction between the mean flow and the primary wave alone, these oscillations do not require any forcing or interacting harmonics to maintain them.Pedlosky (1982) later found that, for situations at the minimum critical shear, the oscillatory behaviour is actually more complex due to non-negligible resonant behaviour.Nevertheless, beyond the minimum critical shear his theory still provides a good example of the mutual interaction between perturbations and the mean flow.Thompson (1987) also used the two-layer QG model and found an oscillation between the mean flow and the primary wave in a truncated five-dimensional closed system.One mode of this system represents the mean thermal wind and the remaining four modes represent the perturbation.However, in this case the period of this oscillation was found to be directly proportional to the ratio between the square roots of heating and friction (Thompson, 1988), implying that diabatic processes are essential for the existence of these oscillations.Realistic values of this diabatic forcing yielded a period of approximately 23 days.Similar oscillations in the two-layer channel model were recovered by Lorenz (1963b), who used a less severe truncation of the model to obtain a 14-dimensional system.Lorenz additionally found that transitions of the system into different more stable or more chaotic regimes were possible, depending on the scale of the eddies and the magnitude of the imposed diabatic forcing.Smaller eddies and larger heating led to a less stable system and eventually chaos.Fleming (2014) found similar characteristics in the Thompson (1987) model, while observing that, when chaos is approached, baroclinic instability becomes very strong and reminiscent of explosive baroclinic development observed in the atmosphere.
In light of the theories above, it is apparent that some periodic or quasi-periodic temporal interaction between the mean shear and eddies is expected in the atmosphere.However, although laboratory experiments can replicate some of the theoretically predicted oscillatory behaviour (e.g.Hide, 1958;Hart, 1972), the baroclinic behaviour of the atmosphere is often deemed to be too irregular and characteristic of non-periodic chaotic behaviour (e.g.Lorenz, 1963a).
Nevertheless, some recent observational evidence suggests a coherent mutual interaction between storm track activity and baroclinicity on various time-scales in both hemispheres (Ambaum and Novak, 2014;Thompson and Woodworth, 2014;Novak et al., 2015;Thompson and Li, 2015).These studies heuristically indicate that a predator-prey relationship may exist between baroclinicity (prey) and storm track activity (predator), as follows.The build-up of baroclinicity due to a diabatic forcing, such as land-sea contrast, orography and oceanic boundary currents (Brayshaw et al., 2011), enhances eddy growth until eddy mixing becomes sufficient to reduce baroclinicity, overriding the effect of the diabatic forcing.Once the baroclinicity is low enough to inhibit any further eddy generation, the storm track activity falls and allows the diabatic forcing to replenish the baroclinicity and the cycle repeats.Such cycling additionally agrees with the quasi-oscillatory behaviour between the mean end eddy energies observed in the Southern Hemisphere (Randel and Stanford, 1985).Despite the qualitative agreement of this predator-prey mechanism with observations, little attempt has yet been made to quantitatively compare the observed variability with existing theoretical predictions.
This article analyses ERA-Interim December-February (DJF) data to uncover an explicit mean oscillatory relationship between the storm track activity and baroclinicity in the storm tracks of the North Atlantic and the North Pacific.Additionally, results from an idealized global circulation model (GCM) with no orography indicate that this behaviour may apply to all localized storm tracks.The detailed properties of the observed oscillatory behaviour are then analysed and compared to those of the Ambaum and Novak (2014) nonlinear oscillator model.
The Ambaum-Novak model is a two-dimensional dynamical system and therefore lacks the potential to become chaotic (Hirsch and Smale, 1974).Despite its simplicity, the authors found a good qualitative agreement between the model behaviour and that of the North Atlantic storm track, using time-lag composites centred around extreme bursts of storm track activity.This article uses more quantitative methods and several storm tracks to produce a more thorough comparison with the Ambaum-Novak model, and to provide more insight into the detailed properties of the observed cycling behaviour.
Section 2 briefly describes the Ambaum-Novak model and some insightful properties of its oscillations.Section 3 discusses the data used and section 4 outlines the construction of a phase space plot from noisy observations.Sections 5 and 6 compare the qualitative and quantitative properties of the Ambaum-Novak model to those of the observed phase space plots in the North Atlantic.Section 7 then investigates whether the oscillating properties also apply to other localized storm tracks and other measures of storm track activity.The remaining section discusses the results and their wider implications.Ambaum and Novak (2014) proposed a two-dimensional model to describe the nonlinear oscillator relationship between baroclinicity and eddy meridional heat flux (a measure of storm track activity) as follows:

The Ambaum-Novak model
where s = −kdT/dy is baroclinicity which is proportional to the negative meridional temperature gradient, and f = kl 2 v T is the scaled heat flux assumed to scale with the squared eddy amplitude.k is a constant and l the meridional wavenumber.s 0 represents a constant eddy dissipation rate and F is a constant diabatic forcing that restores baroclinicity.It is apparent that the constant diabatic forcing is the steady-state value of heat flux and the eddy dissipation rate is the steady-state value of baroclinicity.Investigation of the steady state of the model and its application to the Earth's atmosphere is the subject of a separate study and will not be discussed here further.
Although the Ambaum-  Thompson (1987) and Lorenz (1984), both of which are based on simplifications of the primitive equations (Thompson, 1987;van Veen, 2003).In the Thompson (1987) model, the equation for the mean shear flow is identical to Eq. ( 1), except for an additional linear damping term.The four perturbation equations of the Thompson (1987) model can be written in the form of a matrix which, under vanishing β, is analogous in structure to Eq. (2).Similarly, the equations of the squared perturbations in Lorenz (1984) model can be combined to form an analogous equation system to the Ambaum-Novak model.Although the Lorenz (1984) system includes additional terms of asymmetric heating and linear damping of the mean flow, both of these effects could be argued to be small, based on global observations (Oort, 1964).Ambaum and Novak (2014) show that a useful transformation of the heat flux as allows Eqs ( 1) and ( 2) to be combined to form This equation is a nonlinear oscillator.For small y it linearizes to d 2 y/dt 2 = −2Fy, with a natural frequency of ω 0 = √ 2F.The intensity of the diabatic forcing therefore determines the natural frequency of the oscillator, as was the case for the Thompson (1987) model.At higher amplitudes the nonlinear terms can no longer be neglected.Equation (4) can be simplified using a successive approximation method (Landau and Lifshitz, 1976) that yields a correction to the natural frequency, which for the second-order approximation is where y max is the amplitude of oscillations of the transformed heat flux, y.The model therefore predicts that the period of oscillation increases with the amplitude of heat flux events due to the nonlinearity.Ambaum and Novak (2014) further demonstrate that the symmetry of the model oscillations changes with amplitude.At higher amplitudes, heat flux peaks become more pronounced, more abrupt and less frequent, and baroclinicity exhibits a saw-tooth pattern, falling abruptly during heat flux events but recovering relatively slowly.Such behaviour is more consistent with a relaxation oscillation rather than a simple linear oscillation for which the shape and frequency are constant with amplitude.
Another property of the Ambaum-Novak model is that the amplitude of the perturbations is constant in time, as is evidenced by the zero time derivative of the Lyapunov function: where dE/dt = 0, and V(y) = 2F(e y − y − 1) is the 'potential' of the system, so the system oscillates around its equilibrium in neutral limit cycles.
A further insight into the cycling behaviour can be gained by studying the baroclinicity-heat flux phase space.One can do this by plotting baroclinicity against heat flux, using timeseries from Eqs (1) and (2) obtained for several different amplitudes (i.e.different values of E).An example of such a phase space plot for the two non-dimensionalized variables is shown in Figure 1(a).Time is in the clockwise direction and the thickness of the contours is proportional to the speed of propagation along the trajectories in the non-dimensional phase space.The aforementioned predator-prey relationship emerges whereby baroclinicity increases at times of low heat flux and the subsequent Baroclinicity Heat flux bursts in heat flux erode the baroclinicity to values at which heat flux production is inhibited.The low heat flux then allows the baroclinicity to be slowly replenished again.At high amplitudes, the burst-like nature of the heat flux is additionally evident from the skewness of the contours along the horizontal axis (similar to Fleming, 2014), as well as from the faster speed of oscillation when the heat flux is high.
According to Eq. ( 3), the forcing parameter F is used to compute the transformed heat flux, y, which produces the phase space plot in Figure 1(b).This transformation is applied because it is the y variable that facilitates the mathematical description of the system as an oscillator in Eq. (4).It is this phase space that will be predominantly used to compare the oscillatory properties of the Ambaum-Novak model with those of both the ERA-Interim reanalysis and the GCM data below.Due to the logarithm in Eq. (3), the negative part of y is stretched compared to the untransformed heat flux.The central trajectories are almost circular, similar to those of a linear oscillator.The trajectories then become more skewed and variable in the speed of propagation as the amplitude increases and the system becomes more nonlinear.
The following properties predicted by the Ambaum-Novak model will be tested against atmospheric observations in sections 5 and 6: (i) the spiky nature of heat flux and the slow recovery of baroclinicity, implying two time-scales; (ii) the clockwise oscillation in the baroclinicity-heat flux phase space which becomes more skewed with amplitude; (iii) the similarity of the system to a linear oscillator at small amplitudes; (iv) the decrease of frequency with amplitude in a quadratic fashion; and (v) the increase of speed of the onset and termination of heat flux events with increasing amplitude.

Reanalysis
Observations and 700 hPa and based on perturbations from the 10-day Lanczos filter, in order to isolate transient eddy variability within the storm tracks.Following Hoskins and Valdes (1990), baroclinicity is measured by the maximum Eady growth rate at 775 hPa, s = 0.31f 0 N −1 dU/dZ, where f 0 = 2 sin φ, N = g d(ln θ )/dZ and Z is the geopotential height.
To produce one-dimensional time series of each variable, spatial averaging over the beginning of each storm track was applied.In the North Atlantic the averaging sector spans 35-50 • N and 290-320 • E for heat flux, and 30-50 • N and 270-320 • E for baroclinicity.In the North Pacific the averaging spans 30-50 • N  and 150-195 • E for heat flux, and 25-45 • N and 120-170 • E for baroclinicity.These sectors were chosen to comprise the respective maxima of the two variables (Figure 2), the latitudes of which are fairly constant in time (Novak et al., 2015).
In the North Atlantic, other measures of storm track activity have also been tested.These include eddy kinetic energy {1/2(u 2 + v 2 )} at 250 and 850 hPa, meridional wind variance at 250 hPa, potential vorticity variance at 315 K, geopotential height variance at 250 and 500 hPa, and the mean sea-level pressure variance.The upper-level eddy kinetic energy, which is discussed in more detail below, was spatially averaged over 40-55 • N and 290-330 • E, over its North Atlantic maximum.

Aquaplanet GCM: PUMA
The Portable University Model of Atmosphere (PUMA; Fraedrich et al., 1998) was used here to test whether the predator-prey relationship also occurs in storm tracks of medium complexity models.The model is adapted from the Hoskins and Simmons (1975) spectral GCM with dry dynamics and a linear parametrization of the diabatic forcing and dissipation (Held and Suarez, 1994).The frictional time-scale increases exponentially from 0.5 days at the lowermost level to infinity at the 0.8 sigma level and above, and the thermal damping time-scale ranges from 1.2 days at the lowermost level to 30 days in the uppermost level, following a similar set-up to Fraedrich et al. (2005).
The T42 horizontal resolution (corresponding approximately to a 2.8 • ×2.8 • grid) and 10 vertical sigma levels were used to simulate a single localized storm track without orography.This storm track was produced by imposing a midlatitude heating-cooling dipole with a southwest to northeast orientation (following Frisius et al., 1998).
Daily data of eddy heat flux and baroclinicity (both calculated using the above method) were extracted for 25 winters of the experiment, with the first year being discarded as a spin-up year.The time-mean view of these variables relative to the heating-cooling dipole is displayed in Figure 3.The figure also shows the spatial averaging sector spanning 30-60 • N and 100-180 • E for both variables.

Construction of phase plots from data
This section outlines the method employed to construct baroclinicity-heat flux phase space plots from the relatively noisy time series of the reanalysis data (example time series can be found in Ambaum and Novak, 2014), to extract the average behaviour of the two variables.
Baroclinicity was first scaled as σ = 2s/ √ 2F (according to the Ambaum-Novak model), so as to make the system nondimensional.The value of F/kl 2 is estimated directly from the phase space plots by locating the central value of heat flux which the trajectories encircle.The long-time mean of heat flux is very close to this central value.kl 2 is assumed to be a constant factor, and the choice of its value is discussed in the Appendix.
Next the tendencies of each variable were obtained using the centred-difference method.These were then averaged on a 30 × 30 grid of the phase space, using two-dimensional Gaussian kernel averaging in the phase space (i.e. the closest points to a grid Gaussian weighted average to all data points (of all DJF time series) in this phase space to obtain one average tendency vector per grid point, as in (d).The last step is to plot streamlines using these average tendencies.
point were allocated the highest weights in the average for that grid point).The coarseness of this averaging can be quantified by the standard deviation of the filtering kernel.Using the resulting gridded phase-space vectors streamlines were plotted to display the overall behaviour while ignoring the grid boxes that only contained less than three data points of the original timeseries.
The main steps of this method are depicted in Figure 4.

North Atlantic: qualitative comparison with the model
Examples of the North Atlantic phase space plots are shown in Figure 5.The black arrows are the trajectory streamlines and their thickness is proportional to the speed of oscillation in the scaled non-dimensional phase space.The shading is a two-dimensional histogram, with the darker shading representing a higher data point density extracted from the raw time series.The phase space plots for fine (small filter size), medium and coarse (large filter  constant baroclinicity offset (s 0 ) was chosen to be the centre of the oscillations (0.56 day −1 ), which is close to the time-mean value of baroclinicity (0.57 day −1 ).The clockwise direction of propagation around a single point is reminiscent of the phase space plot of the Ambaum-Novak model (Figure 1(a)).The skew towards high heat flux also concurs with the model and it reflects the spiky nature of the observed heat flux that was previously noted by Swanson and Pierrehumbert (1997) and Messori and Czaja (2013).Furthermore, it is apparent from the thickness of the observed trajectories that the speed of oscillation during the bursts in heat flux does increase with amplitude, as predicted above.
However, there are also some marked differences between the phase space plots of the data and the Ambaum-Novak model.For example, the model does not contain any negative heat flux values and thus its phase space plot can be very compressed near zero for high amplitudes.On the other hand, the reanalysis data show some negative values and less skewness along the heat flux axis, especially for coarser smoothing.Although fine smoothing produces a skewed shape that is more consistent with the model phase space plot, it additionally reveals a more fine-scale structure which is uncharacteristic of the model.Some of this fine-scale structure most likely arises from sampling issues.As is shown in the Appendix, too little or too much filtering seems to cause a departure from the model, and there is an optimal size of the smoothing filter, which produces consistent results.
In order to quantify whether the oscillatory properties of the North Atlantic storm track are consistent with the predicted ones, heat flux was transformed into the y variable.The corresponding phase space plots equivalent to those discussed above are displayed in Figures 5(d)-(f), respectively.It should be noted that the sparse negative values of heat flux were neglected due to the natural logarithm in Eq. ( 3).The phase space plots exhibit similar characteristics to that of the Ambaum-Novak model (Figure 1(b)), though again more structure is apparent for the finest smoothing.The system oscillates around a single point where the data density peaks.The smallest amplitude oscillations are more regular and thus are closer to being linear, concurring with the prediction above.
It is noted that the filtered tendencies are composed of raw tendencies with relatively large standard deviations of sizes comparable to the filtered tendency values (not shown).The reader is thus reminded that the filtered oscillatory behaviour should be interpreted as an average property of the system, rather than representing the system's instantaneous behaviour.
Nevertheless, it is apparent that despite its simplified nature, the Ambaum-Novak model can qualitatively replicate the general features of the average observed behaviour.It is stressed here that this is the case for any (visually) reasonable size and shape of the kernel chosen to smooth the tendencies over the phase space.The next section demonstrates that additional optimization of the filter shape can result in a good quantitative agreement between the model and the observations.

North Atlantic: quantitative comparison with the model
To compare the model and observed data more quantitatively, this section focuses on the case of the medium-sized filter (with standard deviation of 0.65 units of the scaled phase space) and the scaling factor of kl 2 = 3.1×10 −13 K −1 m −1 s −1 .The choice of these parameters yielded appropriate scaling of the baroclinicity (Appendix).Figure 5(e) is the associated phase space plot.
In order to estimate the frequency of the North Atlantic oscillations for trajectories of different amplitudes, the tendencies of all grid points were first interpolated across the phase space using a cubic spline interpolation.The system was then numerically integrated in time (using the fourth-order Runge-Kutta method) from different starting amplitudes of y at s = 0, until a full cycle was completed.The time period (T = 2π/ω) was obtained by summing the number of the 6hourly intervals it took for each trajectory to complete a cycle, and was then converted to angular frequency.Because the observed trajectories are not quite closed, error bars were calculated using the same procedure, but starting the integration in time from  different points along the semicircle of each trajectory in the region of positive y (between the highest and lowest points of baroclinicity).The interquartile range and median of the resulting frequency values were then calculated and displayed in Figure 6(a).Note that due to the finite temporal resolution of the data, sampling at the smallest amplitudes was omitted.
The exact numerically integrated solution of the Ambaum-Novak model and its second-order approximation (from Eq. ( 5)) are also displayed.The model was parametrized using F = 1/2ω 2 0 where ω 0 was estimated from the observed frequency of the smallest amplitude (assumed to be approximately equal to the natural frequency).The tendency of y at y = 0 with increasing amplitude was also plotted (Figure 6(b)), in order to quantify the speed of the onset and termination of heat flux events of different magnitudes.
The frequency of oscillation is almost constant for the smallest amplitudes, approximately marking the natural frequency, ω 0 , of the oscillator.As the amplitude increases further, the frequency decreases and drops off sharply after approximately y max = 1.This pattern of change agrees very well with the idealized model even at high amplitudes.In fact, the agreement is better between the data and the exact model solution than the agreement between the exact model solution and its second-order approximation.The strong agreement between the model and reanalysis data also appears in the plot of y-tendency at y = 0.As the amplitude of oscillations increases, the tendencies of both branches increase in magnitude.The agreement is less good for the positive branch at large amplitudes, most likely because at that point the numerically integrated trajectories, along which the speed was estimated, had already travelled through the region of negative y of the phase space (with the trajectories having been initiated in the region of positive y).Tendencies at negative y values are less well defined due to the logarithm in Eq. ( 3).This behaviour can be contrasted with that of a linear oscillator in which case the frequency would be constant and the speed would increase linearly with amplitude.
The results above and the discussion in the Appendix imply that there is a preferred (natural) frequency of the observed system, though the averaging introduces some dependency.It is therefore of interest to see whether this natural frequency can be detected in the power spectra of the two variables.These were constructed in Figure 7 using the raw November-March reanalysis time series of heat flux and baroclinicity, and the Welsch spectrum estimation method, with a Hanning window and no overlap (similar to Thompson and Li, 2015).The power spectra are the thick solid lines.In order to determine whether the observed variability is significantly different from red noise, red noise timeseries were calculated using lag-1 autocorrelations to produce red-noise spectra (thin solid lines) fitted to the baroclinicity and heat flux spectra.The 99% confidence level (dashed lines) was calculated using the chi-squared statistic.Following Thompson and Li (2015), the red noise fit to heat flux is based on high-pass (30 day) filtered timeseries to remove the dominant low-frequency variability, thus enhancing the fit to the observed power spectrum.However, the choice of the width of this high-pass filter is not unique and the level of significance is sensitive to this choice.
The baroclinicity spectrum behaves very much like the red noise spectrum, but peaks rather more prominently at the very low frequencies, starting from 0.2 rad day −1 and lower, with the higher end of these values being comparable to the natural frequency values found in the above phase space plots.While the heat flux spectrum also peaks at these low frequencies, the spectrum differs from the red noise by a more pronounced higher frequency variability.It would be difficult to conclude from these spectra alone that the timeseries oscillate (as was revealed by the phase space plots above and by the time-lag composites in Ambaum and Novak, 2014).The quasi-periodic nature of the oscillations means that dominant sharp peaks are not necessarily apparent in the power spectra.However, the similarity of the two spectra at the low frequencies (close to the natural frequency value) is an indication that these two timeseries vary most prominently at similar timescales.

Applicability of the Ambaum-Novak model to other eddy variables and storm tracks
In order to investigate whether the above oscillating relationship in the North Atlantic applies independently to other measures of storm track activity, the identical analysis to that above was performed for eddy kinetic energy (EKE= 1/2{u 2 + v 2 }) at 250  and 850 hPa, meridional wind variance (v 2 ) at 250 hPa, potential vorticity variance (PV 2 ) at 315 K, geopotential height variance (Z 2 ) at 250 and 500 hPa, and the mean sea-level pressure variance (MSLP 2 ).Note that, as before, these variables are the 'unfiltered variances', where the products of the perturbations were not subsequently averaged (as is the case in traditional covariance calculations).This last step was omitted in order to investigate the instantaneous variability rather than the 10-day filtered signal.Because all the above storm track activity measures yielded similar (or better) results, only the results from the upper-level EKE analysis and its logarithmic transform, y (EKE) = ln (EKE/F (EKE) ), are presented here.Again, F (EKE) is the steady-state value of EKE and corresponds to the diabatic forcing in the model.The baroclinicity-y (EKE) phase space plot is shown in Figure 8(a).Qualitatively, the data concur with the model, whereby the phase space again exhibits a general oscillating behaviour in the clockwise direction, and the speed of oscillation increases with amplitude during the onset and termination of high storminess events.However, this time the trajectories are less regular and are being suppressed in the negative y territory by an apparent second attracting node, most likely resulting from sampling issues.Although the grey-scale shading shows that most data points are within the primary oscillation, the quantitative comparison between the data and the model is affected by this secondary node.The frequency change with amplitude (Figure 8(b)) is sharper than that of the model and the natural frequency is lower than that estimated from the heat flux in the previous section.Additionally, the y tendency at y = 0 contains substantially larger error bars (Figure 8(c)).
Inspecting the EKE power spectrum in Figure 8(d) reveals that the very low frequency peak that was found in the heat flux case is not apparent.Instead, two other large peaks appear around 0.4 and 0.8 rad day −1 .Since the EKE spectrum does not exhibit dominant variability at the same frequencies as the baroclinicity spectrum, it is unsurprising that the phase space plot is less neat than the heat flux case.Due to the resulting additional higherfrequency structure in the phase space, coarser filtering was applied leading to an underestimation of the natural frequency.
It was found that the agreement with the model improves if the sector used for spatial averaging of EKE is moved to a more upstream location, rather than being centred over the EKE maximum (which is more downstream than the maxima in heat flux and baroclinicity).Furthermore, if the low-level EKE is used, results that are very similar to those obtained when using heat flux can be recovered, with a natural frequency of ω 0 = 0.11 rad day −1 (not shown).Similar results apply for the MSLP variance analysis with ω 0 = 0.13 rad day −1 (not shown).
Having found that the North Atlantic exhibits an oscillatory behaviour with a period of around 30-50 days, it is of interest to investigate whether such behaviour exists in other storm tracks.The present analysis is restricted to localized storm tracks, and so the North Pacific storm track was additionally investigated.Its phase space plot for baroclinicity and the transformed heat flux is shown in Figure 9(a).The general clockwise major circulation with speed of oscillation increasing with amplitude agrees with the picture observed in the North Atlantic.However, the oscillations are less regular and not quite centred around the highest data point density.Additionally, a secondary node appears in the region of negative y and positive excess baroclinicity.The natural frequency (Figure 9(b)) is much lower than in the North Atlantic, even if the filtering kernel is decreased to the smallest value that allows a coherent oscillation around the phase space (not shown).Despite these inconsistencies, the quantitative agreement between the model and data is still reasonable for both the frequency and the negative y-tendency changes with amplitude (Figures 9(b) and  (c)).The positive branch of y-tendencies is evidently suppressed by the additional structure in that region of the phase space.The secondary circulation forces the trajectories to deform and enter the more central region of the phase space plot, where lower y-tendencies persist (i.e.where the energy of the oscillator is lower).Investigating the power spectra (Figure 9(d)) reveals that the most dominant heat flux modes have an angular frequency of 0.55 rad day −1 .In contrast, the Pacific baroclinicity spectrum peaks at frequencies that are more than five times smaller, similar to the baroclinicity frequency peaks found in the North Atlantic.As a result, the oscillating behaviour is less well defined in the phase space.The maximum of the storm track activity in the North Pacific is situated further downstream from the maximum of the baroclinicity (Figure 2).As in the Atlantic upper-level EKE case, it is plausible that the lowered spatial co-location between baroclinicity and heat flux in the Pacific degrades their mutual interaction and allows other (higher-frequency) variability to dominate the system.This can be confirmed by conducting the same analysis but averaging both baroclinicity and heat flux over the same geographical sector in the North Pacific.This yields a very good quantitative agreement with the Ambaum-Novak model (not shown).
Finally, the heat flux-baroclinicity phase space of the storm track of the 25 winters of the PUMA experiment is studied.Again, the phase space plot in Figure 10(a) exhibits a clockwise circulation around the maximum of the data point density, with increasing speed of oscillation with amplitude.However, the circulation trajectories are somewhat tilted in this case, implying that the heat flux peak occurs slightly earlier in relation to the dip in baroclinicity compared to the two real storm tracks.The more quantitative oscillatory properties, namely the changes with amplitude in frequency and y tendency at y = 0, are displayed in Figures 10(b) and (c).Again, the agreement with the Ambaum-Novak model is very good for both properties, and the natural frequency is consistent with that observed in the North Atlantic, despite the apparently large size of the kernel used to average over the phase space.Studying the power spectrum (Figure 10(d)) , the highest peaks are near 0.9 rad day −1 , though both variables also exhibit strong peaks at higher frequencies, which no doubt contribute substantially to the average behaviour in the phase space.Qualitatively similar behaviour can also be found when using the zonal mean to define the eddies, and also for perpetual equinox conditions (not shown).
It is noted that a better quantitative agreement with the Ambaum-Novak model can be obtained by increasing the temporal and spatial resolution of the model as well as by creating a more locally confined storm track (as will be shown as part of a forthcoming study).However, the model set-up here was intentionally chosen to be highly simplified in order to indicate that the predator-prey behaviour is likely present in any localized storm track.

Summary and discussion
The localized storm tracks of the North Atlantic, the North Pacific and the idealised GCM exhibit a predator-prey relationship between storm track activity (predator) and baroclinicity (prey).This relationship is based on the competing effects of high storm track activity reducing meridional temperature gradients (and thus baroclinicity) and a restoring forcing that enhances the baroclinicity when storm track activity is low.In time, this relationship leads to burst-like oscillations in storm track activity and saw tooth-like oscillations in baroclinicity (characterized by a fast reduction and a slow restoration).The observation of this behaviour in all of the studied storm tracks, as well as for several different measures of storm track activity, suggests that this is a feature of all localized storm tracks.This nonlinearly oscillating relationship was observed by Ambaum and Novak (2014) for approximately weekly timescales, though the present analysis of the system's phase space indicates that the dominant scales are closer to 30-50 days, which is also reflected in the similarity of low-frequency spectral peaks of baroclinicity and storm track activity.These time-scales are more akin to those observed by Randel and Stanford (1985) for oscillations between mean and eddy energies in the Southern Hemisphere, as well as to those associated with the Baroclinic Annular Modes in both hemispheres (Thompson and Woodworth, 2014;Thompson and Li, 2015).All of this observational evidence further supports the existence of the predator-prey relationship between storm track activity and its growth rate at the longer time-scales.As a result of this relationship, the system seems to undergo transitions between a baroclinicity-eroding state of intense storm track activity and a more stable baroclinicity-replenishing state, as suggested by Stone (1978).These arguments are also reminiscent of the index cycle   (Namias, 1950) which is characterized by transitions between regimes of zonal and wavy jet structures.Despite the cycling behaviour of the system, power spectra of both baroclinicity and storm track activity do not show unambiguous evidence of pure periodicity in the winter season (agreeing with Thompson and Li, 2015), but instead the analysis above implies a more quasi-periodic behaviour.Therefore care should be taken when computing composites of different large storm events, as the periodicity may be masked by the different time-scales of the events (e.g.Ambaum and Novak, 2014).For this reason, the properties of the observed oscillations were studied here in the baroclinicity-storm track activity phase space.The phase space analysis summarizes the oscillatory properties of the system at different amplitudes, without having to assume that different-amplitude events occur on similar time-scales.
A more quantitative investigation revealed some sensitivity to the size and shape of the filter that was used to smooth the raw timeseries over the baroclinicity-heat flux phase space.Since the raw data series are relatively noisy, averaging is necessary to recover their overall behaviour, though excessive averaging can damp the desired structures.An optimal size was therefore selected by mutual tuning of the detailed properties of the data and those of the Ambaum-Novak model.However, even in the absence of such tuning, any reasonable size and shape of the filter produces the same qualitative features.
It was found that the average time period of oscillation and speed of onset and termination of storm events both increase with the amplitude of those events, implying that more extreme events occur more suddenly but less frequently.This may be of considerable use in the context of climate change studies.For example, modelling studies (e.g.Pinto et al., 2009) have found that in warmer climates North Atlantic storms will become less frequent but retain their current amplitude.This suggests that the system would enter a new regime in the phase space with a lower diabatic forcing (F).This is consistent with the idea that diabatic processes, such as boundary currents and land-sea contrast, would be less efficient at sharpening the lowlevel meridional temperature gradient in a climate with a warmer Arctic (e.g.Barnes and Screen, 2015).Preliminary results show that the predator-prey behaviour exists in the Coupled Model Intercomparison Project Phase 5 (CMIP5) climate models and a forthcoming study aims to present the response of this behaviour to the climate change forcing more explicitly.
It should be noted that, although the oscillation in the phase space was observed in all storm tracks studied here, it was somewhat less defined in the North Pacific.Nakamura and Sampe (2002) have found that, as the East Asian winter monsoon intensifies, the enhanced westerlies trap the Pacific synopticscale eddies within the subtropical jet core and guide them southward and downstream, away from the maximum low-level baroclinicity.In doing so, they reduce the interaction between the upper and lower levels which is essential for baroclinic growth, as is consistent with the barotropic governor mechanism (James and Gray, 1986).Consequently, the Pacific storm track peaks in intensity further downstream of the baroclinicity maximum.This would move a substantial part of the eddy erosion of temperature gradients outside the region of maximum baroclinicity, where (instead of the diabatic heating sharpening temperature gradients) downstream fluxing of recycled energy from previous flow structures would dominate the eddy growth (Orlanski and Katzfey, 1991;Chang and Orlanski, 1993).The worse co-location of highest storm track activity and maximum baroclinicity would make the interaction between baroclinicity and eddies weaker and the predator-prey relationship would be less clear.In addition, some measures of storm track activity, such as PV variance, only exhibit a maximum over the North Atlantic in the winter time mean, because of the Pacific midwinter minimum in storm track activity.It was therefore more difficult to identify the correct spatial averaging domain over the North Pacific storm track.
In contrast, the stationary forcing and weaker winds in the North Atlantic allow the storm track to be more longitudinally confined and the heat flux (and thus the region of maximum mixing of temperature) is more co-located with the maximum of baroclinicity.Any changes in one variable can therefore be immediately felt by the other.As a result, the theory of baroclinic instability is more applicable here and the predator-prey behaviour is more detectable.It is noted that in general the barotropic   (1970) and the diabatically forced system of Ambaum and Novak (2014), showing (a) the non-dimensional potentials of the two systems, with the natural frequency of the Ambaum-Novak model being tuned to the frequency of the Pedlosky (1970) model near the potential minimum; and (b) the corresponding nondimenisonal phase space plot of the Pedlosky (1970) system.A represents the (third-order) amplitude of the baroclinic wave and s represents the dA/dt which relates to the vertical wind shear of the mean flow.The shading in (b) emphasizes different initial amplitudes, with darker shading marking higher values.Note that the potential of the Ambaum-Novak model in (a) was shifted to coincide with one of the potential minima of the Pedlosky (1970) model for better comparison.interaction between eddy momentum fluxes and the mean flow occurs further downstream where eddies begin to decay and where baroclinicity is relatively low.This 'barotropic' eddy-mean flow interaction should be distinguished from the interaction between mean baroclinicity and eddy heat fluxes, which occurs mostly during the growing stage of eddies, more upstream and closer to the maximum baroclinicity (e.g.James, 1994).While eddies may propagate downstream beyond the averaging region, this advection can be modelled as part of eddy dissipation.
The other focus of this article was to compare the above observed behaviour to the low-complexity yet nonlinear oscillator model of Ambaum and Novak (2014).Even though the model lacks a large part of the complexity of the real system (such as the slight convergence or divergence of the circulating trajectories of the phase space), it still proved capable of reproducing the general observed behaviour.For example, it exhibits the clockwise circulation in the baroclinicity-heat flux phase space.It also predicts the observed approximately quadratic increase of time period with amplitude and the approximately linear increase of the speed of onset and termination of heat flux events as the amplitude increases.For an appropriate filter, these properties show a strikingly good agreement between the model and the observed system.In some cases, the agreement between the observed data and the exact numerical model solution was better than the agreement between the exact solution and its second-order approximation.Even though the Ambaum-Novak model was parametrized with the observed natural frequency, such changes of the system with amplitude are not prescribed in the model and its agreement with the real system is not trivial.
It was also found that in general low-level measures of storm track activity produce a neater comparison with the Ambaum-Novak model.The upper levels may have more memory of pre-existing flow structures that propagate from the upstream regions, which may introduce additional variability into the phase space.The clearest picture was found using lowlevel heat flux, low-level EKE and MSLP variance, averaged over the region of maximum baroclinicity.
Furthermore, the phase space of the reanalysis data in conjunction with the Ambaum-Novak model imply that the recovery of the baroclinicity occurs slowly (with a time-scale 2π/ √ F estimated between 30 and 50 days, corresponding to the range of natural frequencies found for different coarseness of averaging) relative to the eddy dissipation (with a time-scale 1/s 0 of around 2 days).The difference in the time-scales can also be observed from the burst-like nature of storm track activity (Messori and Czaja, 2013), as well as from the speeds of oscillating trajectories across the phase space.These speeds are very slow during low storm track activity and fast during high activity.This supports Stone's (1978) baroclinic adjustment theory, whereby the system fluctuates around a state of baroclinic neutrality, between two regimes when eddies are either efficient or inefficient at reducing meridional temperature, so as to prevent the system from becoming too stable or unstable.Storm tracks can be seen to be baroclinically neutral in the steady state by definition.However, because the efficient-eddying regime occurs in short-term bursts whilst the baroclinicity-recovering regime occurs slowly, the storm tracks are marginally stable (with respect to the vigorous eddy mixing) for the majority of the time.
The Ambaum-Novak model further suggests that diabatic forcing is crucial in determining the frequency of the oscillations, as well as setting the steady-state value of the storm track activity.On the other hand, the steady-state baroclinicity is predicted to be set by eddy dissipation and remain insensitive to the diabatic forcing.This insensitivity of baroclinicity to diabatic forcing has been observed previously in the atmosphere and in idealized simulations, though the precise mechanisms are still unclear (Stone, 1978;Schneider and Walker, 2006;Zurita-Gotor, 2008;Jansen and Ferrari, 2013).While investigation of the explicit response of steady-state storm tracks to diabatic forcing and eddy dissipation is beyond the scope of the present study, it is suggested here that the Ambaum-Novak model may provide some useful insight into this behaviour.
It was additionally indicated that most of the reanalysis data lie within the near-linear regime, where the exact numerical and second-order approximated solutions only just start to diverge from each other.This observation strengthens the scope for using the weakly-nonlinear theory for the baroclinic instability of localized storm tracks.
However, it should be noted that the mechanism discussed here is different to that of Pedlosky (1970), who used the weakly nonlinear theory to predict oscillatory behaviour between the perturbation amplitude and the underlying baroclinicity.In this study the system is dependent on the overshooting of the eddy mixing effect (and thus phase angle reversal) which would allow eddies to return energy to the mean flow.Consequently, the two quantities oscillate inertially whilst propagating downstream with the perturbation.On the other hand, the model used in this study (as well as that of Thompson, 1987) does not include such a mutual feedback.Once the eddies decay to almost zero activity, their phase angle never attains negative values, and so the eddies do not have the ability to return available potential energy back to the flow.This is where the necessity of the diabatic heating becomes apparent for this model, as without it the baroclinicity would not be replenished.In other words, Pedlosky (1970)  dual role: they erode baroclinicity when it is high, and replenish it when it is below the critical value.Ambaum and Novak (2014) eddies are only able to erode it.The different potentials of the two systems can be compared in Figure 11(a).It is apparent that the two systems exhibit a similar behaviour near their steady-state equilibria, where their potentials are at their minima and their oscillations are nearly linear (not shown).For larger values, the two systems differ significantly, as can be observed by studying their respective phase space plots (Figures 1(b) and 11(b)).
It can be concluded that the predator-prey relationship between baroclinicity and storm track activity appears to be a general feature of all localized storm tracks.This relationship is characterized by both quantities oscillating quasi-periodically around their steady states due to the baroclinic adjustment of storm tracks.This agrees with the weakly nonlinear theory that the assumption of a constant background wind shear is insufficient to adequately describe storm track dynamics.sensitive: the size and the shape of the Gaussian smoothing filter that was used to average over the phase space.To demonstrate the effect of changing filter size, Figure A1(a) shows the change in frequency with amplitude for cases when fine, coarse and medium smoothing of the phase space was applied.It is apparent that excessive smoothing makes the decrease in frequency with amplitude very slow and it yields a lower natural frequency compared to the cases when smaller-sized filters are used.However, the smallest-sized filter exhibits a decrease in frequency that is no longer monotonous, reflecting the additional finescale structure of the phase space plot (Figures 5(a  filter.It should be noted that, while the analysis of the frequency change with amplitude is rather sensitive to the filter size, the analysis of the y-tendency change with amplitude during the onset and termination of heat flux events is less sensitive (Figure A1(b)).
The shape of the smoothing filter is closely related to the value of the scaling factor kl 2 , which is used to scale baroclinicity to its non-dimensional equivalent (σ = 2s/ √ 2F).A proper scaling of the system should allow the use of a circular filter in the nondimensional phase space.However, if the scaling factor was too small, then the baroclinicity would be too stretched relative to the y variable (and the filter too compressed along the baroclinicity axis).In the case of an excessively stretched baroclinicity, the Ambaum-Novak model predicts that the trajectories would travel too far along the baroclinicity axis in the phase space into regions of higher energy (or Lyapunov function, E) where the speeds are higher, which could make the overall oscillation period shorter, and vice versa.In the data, this effect of baroclinicity stretching on the time period was found to be true for all reasonable filter sizes and amplitudes tested (Figure A2).
As a further example, Figure A3 shows how the natural frequency, obtained using the medium-sized filter, decreases as the scaling factor increases.This figure also shows the analytical relationship that ω 0 = √ 2F = 2ckl 2 v T predicted by the Ambaum-Novak model, where c is an empirically determined constant which scales the time-mean heat flux to make it equal to the central point in the phase space plot.In this case c = 1.03 since the time-mean heat flux lies very close to the central value.It is reasonable to expect that the convergence of the predicted analytical relationship with the observed (filter-dependent) values would indicate the appropriate value for the scaling factor.For the medium-sized filter, the convergence is at ckl 2 = 0.9×10 −13 K −1 m −1 s −1 .The corresponding variability of the Lyapunov function of a selection of trajectories of different amplitudes is shown in Figure A4(a).Each trajectory represents integrations in time (using the fourth-order Runge-Kutta method), starting from the maximum heat flux when s = 0, until an entire cycle has been completed.It is apparent that when the trajectories reach the heat flux maxima and minima, the Lyapunov function for this particular scaling factor decreases substantially compared to the times when the trajectories are near the baroclinicity maxima and minima.This implies that the baroclinicity is still a little overstretched.The system was therefore further tested for slightly higher values of the scaling factor, with kl 2 = 3.1×10 −13 K −1 m −1 s −1 yielding the least pronounced fluctuations in the Lyapunov function, as shown in Figure A4(b).Evidently, the Lyapunov function also varies at times when the trajectories are not at the maxima and minima of heat flux or baroclinicity, reflecting the additional structure in the data that is not characteristic of the Ambaum-Novak model.The overall Lyapunov function also decreases slightly in time mirroring the trajectory convergence towards the centre, as is shown in the corresponding phase space plot (Figure 5(e)).As noted in the main text, low-amplitude oscillations are almost linear and so they exhibit smaller fluctuations in the Lyapunov function, as well as a reduced overall decay of the latter.
It should be pointed out that, if one explores the scaling of the kl 2 factor, it can be found empirically that the value of the Eady scale k = 0.31g/(NT 0 ) ≈ 0.78 m K −1 s −1 and c = 1.03.For kl 2 = 3.1×10 −13 K −1 m −1 s −1 , this would therefore yield the meridional extent of the storm track (or half of the wavelength) of 5000 km.This seems excessive, but is the right order of magnitude.It should be noted that the maximum Eady growth rate is known to overestimate the length-scales of cyclones (James, 1994), meaning that it may require larger than observed eddies to produce a more realistic baroclinicity.In addition, it is apparent that the data analysis includes a substantial amount of averaging of the phase space due to the additional noise of the data.It was shown that averaging reduces the natural frequencies and would thus increase the spatial scales of the disturbances estimated from the above scaling.
It was shown above that the filter can be tuned to adhere to the Ambaum-Novak model and the remnants of the data variability can be classed as noise.However, the model itself can be parametrized to adhere to the data.Thus it could be argued that there is no perfect precise size and shape of the filter.The results here merely state that a consistency between the Ambaum-Novak model and the observed data can be achieved by their mutual tuning.This consistency is by no means obvious and not always perfectly achievable, as seen in the North Pacific storm track case in the main text.Additionally, a moderate consistency between the qualitative behaviour of the data and the model can be achieved without the mutual tuning, as long as a reasonable filter size and shape are selected (i.e. a filter that reflects the data point distribution adequately along both axes of the phase space).This consistency of the idealized Ambaum-Novak model with the observed storm tracks suggests that similar dynamics underline them.For more unequivocal results, it may be useful to analyse the remaining data variability (here viewed as noise) and isolate it from the phase space oscillations in medium-complexity models.This is beyond the scope of this article.

Figure 1 .
Figure 1.Examples of phase space plots for the numerically integrated Ambaum-Novak model for (a) the instantaneous heat flux, f , and (b) the corresponding transformed heat flux, y = ln(f /F), where F is the diabatic forcing and s is the excess baroclinicity.The system was scaled to be non-dimensional.The line thickness is proportional to the relative speed of the oscillations in the non-dimensional phase space and time is in the clockwise direction.

Figure 2 .
Figure 2. Time-mean DJF plots of heat flux (black, 10 and 20 K m s −1 contours) and baroclinicity (grey, 0.5 and 0.6 day −1 contours) calculated from the ERA-Interim reanalysis.The averaging sectors for heat flux (thin black boxes) and baroclinicity (thin grey boxes) are also shown.

Figure 3 .Figure 4 .
Figure 3. Time-mean DJF plots of heat flux (black; 10, 15 and 20 K m s −1 contours) and baroclinicity (grey; 0.6, 0.7, 0.8 and 0.9 day −1 contours) calculated from the PUMA experiment.The averaging sector is shown by the thin black box and the location of the heating-cooling dipole is shown in light grey.
size) smoothing are shown in Figures 5(a)-(c), respectively.The c 2017 The Authors.Quarterly Journal of the Royal Meteorological Society published by John Wiley & Sons Ltd on behalf of the Royal Meteorological Society.

Figure 5 .
Figure 5. Phase space plots for (a)-(c) the observed instantaneous heat flux, fv T and (d)-(f) the transformed heat flux, y = ln(f /F),where F is the diabatic forcing (in this case around F/kl 2 = 28 K m s −1 ) and s is the excess baroclinicity.Line thickness is proportional to the speed of oscillation in the scaled non-dimensional phase space and the shading is the 2D histogram (showing values of 10, 50 and 100+ data points).The plots in (c,f) are averaged using a coarse filter (standard deviation: 1.1), while finer filters have been used for (b, e) and (a,d) (0.65 and 0.2 respectively).The radii of the grey ellipses at the bottom right corner mark the standard deviations of the averaging Gaussian filter.Only the filtered tendencies that were surrounded (within a distance of one grid point) by two or more raw data points are plotted.

Figure 6 .
Figure6.The change with increasing amplitude in (a) the angular frequency of oscillations (ω) and (b) the tendency of y at y = 0, showing the exact numerical model solution (solid line), the approximated second-order model solution (dashed line) and the reanalysis data filtered as in Figure5(e) (points).Note that the sampling of the smallest amplitudes below y max = 0.2 is omitted due to the data resolution being too coarse to resolve such small changes there.

Figure 7 .
Figure 7. Power spectra of the North Atlantic heat flux (black thick solid line) and baroclinicity (grey thick solid line) time series from the raw reanalysis timeseries, with their respective fitted red noise spectra (thin solid lines) and the 99% confidence interval (dashed lines).

c
2017 The Authors.Quarterly Journal of the Royal Meteorological Society published by John Wiley & Sons Ltd on behalf of the Royal Meteorological Society.

Figure 8 .
Figure 8. Analysis of the variability of y (EKE) (calculated at 250 hPa and averaged between 40-55 • N and 290-330 • E) and baroclinicity in the Atlantic storm track, showing (a) the phase space plot with the 2D histogram (grey shading, with contour intervals of 10, 50 and 100 data points) and the filter standard deviation (grey ellipse); (b) the angular frequency change with amplitude, as in Figure 6(a); (c) the y (EKE) -tendency change with amplitude at y (EKE) = 0, as in Figure 6(b); and (d) the power spectra of the raw baroclinicity (grey) and EKE (black; constructed as in Figure 7) time series.

Figure 9 .
Figure 9.As Figure 8, but for baroclinicity and heat flux in the North Pacific storm track.
c 2017 The Authors.Quarterly Journal of the Royal Meteorological Society published by John Wiley & Sons Ltd on behalf of the Royal Meteorological Society.

Figure 10 .
Figure 10.As Figure 8, but for baroclinicity and heat flux in the storm track of the PUMA experiment.
c 2017 The Authors.Quarterly Journal of the Royal Meteorological Society published by John Wiley & Sons Ltd on behalf of the Royal Meteorological Society.

Figure 11 .
Figure11.Comparison of the unforced inertial oscillatory system ofPedlosky (1970) and the diabatically forced system ofAmbaum and Novak (2014), showing (a) the non-dimensional potentials of the two systems, with the natural frequency of the Ambaum-Novak model being tuned to the frequency of thePedlosky (1970) model near the potential minimum; and (b) the corresponding nondimenisonal phase space plot of the Pedlosky (1970) system.A represents the (third-order) amplitude of the baroclinic wave and s represents the dA/dt which relates to the vertical wind shear of the mean flow.The shading in (b) emphasizes different initial amplitudes, with darker shading marking higher values.Note that the potential of the Ambaum-Novak model in (a) was shifted to coincide with one of the potential minima of thePedlosky (1970) model for better comparison.

Figure A1 .
Figure A1.Comparison of the effects of largest, medium and finest filtering on (a) the frequency change with amplitude and (b) the change in the y tendency at y = 0 with amplitude.This comparison is based on the North Atlantic heat flux and baroclinicity, and the scaling factor kl 2 = 3.1×10 −13 K −1 m −1 s −1 was selected.The numbers indicate the standard deviations of the filters.The medium filter changes are equivalent to Figure 6.

Figure A3 .
Figure A3.Natural frequency change with the scaling factor, kl 2 , tested for the North Atlantic heat flux and baroclinicity phase space (points), using the medium-sized filter.The analytical prediction of the relationship based on the Ambaum-Novak model is also shown by the solid line (see text for more detail).
) and (d)).As mentioned earlier, one must compromise between losing the desired structure and gaining additional fine-scale noise.The focus in the main text is therefore on the optimal medium-sized c 2017 The Authors.Quarterly Journal of the Royal Meteorological Society published by John Wiley & Sons Ltd on behalf of the Royal Meteorological Society.

Figure A4 .
Figure A4.Change of the Lyapunov function along the trajectories of one cycle in the phase space plot for different starting amplitudes of y max for (a) kl2 = 0.9×10 −13 K −1 m −1 s −1 and (b) kl 2 = 3.1×10 −13 K −1 m −1 s −1 .The legend indicates the initial amplitude (y max ) of each trajectory.The points at which the trajectories reached the positive and negative heat flux at s = 0 (squares and crosses, respectively), and positive and negative baroclinicity at y = 0 (diamonds and circles, respectively) are also shown.
Quarterly Journal of the Royal Meteorological Society published by John Wiley & Sons Ltd on behalf of the Royal Meteorological Society.
eddies have a c 2017 The Authors.Quarterly Journal of the Royal Meteorological Society published by John Wiley & Sons Ltd on behalf of the Royal Meteorological Society.