A nonlinear oscillator describing storm track variability

We construct a two‐variable model that describes the interaction between local baroclinicity and eddy heat flux in order to understand aspects of the variance in storm tracks. It is a heuristic model for diabatically forced baroclinic instability close to baroclinic neutrality. The two‐variable model has the structure of a nonlinear oscillator. It exhibits some realistic properties of observed storm track variability, most notably the intermittent nature of eddy activity. This suggests that apparent threshold behaviour can be more accurately and succinctly described by a simple nonlinearity. An analogy is drawn with the triggering of convective events.


Introduction
The nonlinear phase of a baroclinic lifecycle is usually understood as an extension of the linear instability problem, in that one is presented with a baroclinically unstable background state in which an instability evolves.The main difference is that, in a nonlinear life-cycle experiment, the instability is allowed to grow to finite amplitude and this will modify the background state.However, it can be argued that this is an unusual set-up to study naturally occurring baroclinic lifecycles, as very unstable basic states cannot be widespread by virtue of their instability.
By analogy, when studying the evolution of thermals, we would not normally start with a statically unstable state and let the instability evolve.We would start with a statically stable state and then slowly heat the bottom surface.As a result, the basic state would at some point cross the instability threshold and thermals would develop.Such a set-up is appropriate for studying a dry convecting boundary layer.The two key differences between this set-up and the set-up starting from a statically unstable state are that the system (i) is diabatically forced and (ii) mostly resides close to marginal stability.We argue that the nonlinear phase of baroclinic instability must also be understood in such a way: diabatically forced and close to baroclinic neutrality.In this article, we describe a heuristic model for baroclinic instability and storm track variability based on this viewpoint.
Storm tracks are maintained by a balance between the production of potential energy, which can be identified by baroclinicity, and the erosion of this potential energy by the eddies in the storm track, which can be identified by transient eddy heat fluxes (Chang and Orlanski, 1993).Maxima in either of these fields can be used to diagnose the location of storm tracks (Hoskins and Valdes, 1990).However, there is a potential ambiguity regarding the maxima of baroclinicity in relation to the maxima of the † The copyright for this article was changed on 9 April 2014 after original online publication.
eddy activity: a high baroclinicity presumably means a favourable environment for eddy growth, but a high number of eddies would, all else being equal, correspond to a low baroclinicity.The first aspect is, for example, reflected in studies of energy balance models (e.g.Sellers, 1969;Hwang and Frierson, 2010), where energy fluxes increase linearly with temperature gradients.The second aspect reflects the well-established view that eddies feed on the available potential energy stored in the temperature gradient (Holton and Hakim, 2012).Here we construct a simple two-variable model describing these conflicting aspects.We find that the model corresponds to a nonlinear oscillator, which exhibits some important properties of the observed variability of the storm track, such as the observed intermittent, bursting behaviour of the meridional heat flux in storm tracks (Swanson and Pierrehumbert, 1997;Messori and Czaja, 2013a).
Our model is perhaps the simplest possible model that describes the nonlinear growth of an instability in a background state, where the feedback of the instability on the background state is taken into account.It may appear surprising that this feedback in our model changes exponential growth to oscillatory behaviour, although this has been observed in previous studies of frontal instability (e.g.Spall, 1997), as well as in weakly nonlinear extensions of unforced barcolinic instability theory (Pedlosky, 1987, ch 7).In other words, the feedback to the background state becomes strong enough to dominate, and in our case indeed reverse, the exponential growth of the instability.
Such behaviour is normally associated with states near marginal instability.Indeed, the two competing interactions between baroclinicity and eddies admit an equilibrium state close to neutrality, where the local diabatic production of instability is consumed rapidly by the instability itself.Such a near-neutral state is thought to reflect the mean state of the atmosphere (Stone, 1978;Hall and Sardeshmukh, 1998) and is consistent with a predictive scaling relation for the height of the tropopause (Lindzen, 1993).However, the latter arguments do not describe the variance of the state.
One way to understand the variance induced by the competing effects of eddy growth on high baroclinicity and decay of baroclinicity due to those eddies is to think of a gradual diabatic build-up of instability which, above a certain threshold, will be released by the eddies.Such models are perhaps reminiscent of archetypal models of collapsing sand-piles, earthquakes, percolation networks and many more model systems, associated with phase transitions, scaling behaviour and self-organized criticality (Burridge and Knopoff, 1967;Bak et al., 1987).Here, we do not explore any formal associations between such models and our proposed model, except to indicate that our model, like models of self-organized criticality, also includes an external supply of instability that is released by eddy growth.In contrast to models of self-organized criticality, our release mechanism is much more simple and does not admit any scaling behaviour.Ideas of threshold behaviour have also been successfully applied to the release of convective instability (Peters and Neelin, 2006;Yano et al., 2012).
Such models of near-critical variance are complex and can lead to rich behaviour.Here we simplify this picture drastically by going to first principles on what we know about the expected interaction between baroclinicity and eddies.It turns out that this naturally leads to an oscillator equation with a simple nonlinearity.The model is derived and discussed in the next section.In section 3, we present a first comparison of our nonlinear oscillator model with observed data.Section 4 contains some concluding remarks and an analogy with convective instability.

A nonlinear oscillator model
Here, we derive and discuss a two-variable model describing the interaction between baroclincity and baroclinic eddies.We take the baroclinicity s to be proportional to the horizontal temperature gradient near the surface, s = −kT y . (1) The baroclinicity represents the maximum growth rate of baroclinic eddies, such as in the archetypal Eady model for baroclinic instability.The fastest growing mode will likely dominate the baroclinic activity locally, so that the baroclinicity will provide the correct dominant time-scale for the growth of the eddies.Alternatively, the growth could be locally dominated by non-modal growth, but in this case the baroclinicity also provides the relevant bulk time-scale for the baroclinic growth.This is mainly an expression of the fact that all baroclinic growth extracts its energy from the vertical wind shear in the background state and this vertical shear then provides the physically relevant time-scale.
The evolution of the baroclinicity can be described in terms of a mean imposed forcing and an erosion by eddy heat fluxes.The imposed forcing of the baroclinicity represents the diabatic tendency to increase the temperature gradient between Pole and Equator and independent dissipative tendencies due to eddies generated upstream of the high-baroclinicity area (e.g.lee cyclogenesis over North America for the North Atlantic storm track), as well as the tendencies due to the mean winds, such as those described by a frontal deformation field perhaps imposed by planetary stationary waves.The effect of an imposed deformation field on linear instability of baroclinic (Spall, 1997), barotropic (Dritschel et al., 1991) and surface Rossby waves (Harvey and Ambaum, 2010) can have a complicated interaction with the linear instability.A deformation field typically increases the potential for instability but at the same time kinematically suppresses the growing wave.For example, in the case described in Harvey and Ambaum (2010) these two competing effects become important at different times in the evolution.In the present set-up, we simplify the combined effect of diabatic heating and imposed large-scale flow by a simple, in our case constant, positive tendency of the baroclinicity, denoted F: where in the last equality we assumed that the meridional extent of the eddy heat flux is dominated by a meridional scale of 1/l.Hereafter we will define the scaled eddy heat flux as The eddy heat flux itself depends on the precise structure of the eddies, but it scales with the square of amplitude of eddies.The baroclinic instability therefore leads to a growth rate of 2s for the eddy heat flux.We also expect dissipative processes to occur and we will model these as a fixed linear decay rate of 2s 0 (see Hall and Sardeshmukh, 1998).Thus, we find the following set of two equations for the evolution of the baroclinicity and the eddy heat flux: The second equation introduces a nonlinearity, because the growth rate of the heat flux varies with the heat flux itself.The dissipative processes described by s 0 can be absorbed in a constant offset of the baroclinicity.It turns out that the solution for the excess baroclinicity, s − s 0 , is a symmetric oscillation, which means that s itself oscillates around a positive value of s 0 .
This set of equations cannot be solved analytically in terms of tabulated functions, but we can transform the set such that their behaviour becomes clear.This is achieved by transforming the heat flux f to a new flux variable, y, as (5) a transformation that is well-defined because f can be assumed to be positive-definite.With the new heat flux variable, the above set of equations becomes This set of first-order equations corresponds to the following nonlinear oscillator equation in y: For small y, i.e. situations where the eddy heat flux f closely compensates for the imposed forcing F, we find a linear oscillation equation with frequency √ 2F and from Eq. (6b) it follows that the baroclinicity s oscillates about s 0 with the same frequency and a phase of π/2 ahead of the heat flux, as expected from causality arguments.
To understand the behaviour of this equation for larger deviations, we can write the above nonlinear oscillation in terms of a potential V as We include an offset in the potential such that the zero point of the potential corresponds to y = 0. Figure 1 shows the potential as a function of y.The potential is positive-definite, everywhere concave and asymmetric.
From the potential form of the equation it becomes clear that the nonlinear oscillator has a Lyapunov function, or conserved 'energy' E, defined by The different behaviours of the nonlinear oscillator are parametrized by the value of this energy.For example, the extreme values of y are defined by the implicit equation V(y) = E.This is a transcendental equation, which has trivial asymptotic solutions for small and large E. For small E and y, the potential is approximated by V(y) ≈ Fy 2 , so that the extremes of the oscillation correspond to y = ± √ E/F.In other words, y oscillates for the nonlinear oscillator in Eq. ( 8), scaled with 2F, as a function of the transformed heat flux y.The dashed lines correspond to the two asymptotic forms for small y, V(y)/2F ≈ y 2 /2, and for negative y, V(y)/2F ≈ −y − 1.
symmetrically around y = 0 with a frequency of √ 2F.Such an oscillation corresponds to the untransformed heat flux f varying between F ± √ EF with a frequency of √ 2F.Equation (6b) then indicates that s also oscillates with a frequency of √ 2F, between values of s 0 ± √ E/2.The behaviour for large E is more interesting, with different asymptotic behaviour for positive and negative values of y.Using the language of a mass oscillating on a nonlinear spring, we find that for positive y the spring is hardening and for negative y the spring is softening and asymptoting to a constant-force spring.The extremes of the oscillation in y are approximated by [−E/2F − 1, ln(1 + E/2F)], so the eddy heat flux f varies in the range [Fe −E/2F−1 , F + E/2].(The top value underestimates the real maximum by a sizeable fraction.)For large E, therefore, the heat flux remains low (negative y) for a substantial period of time and then increases rapidly to a value of F + E/2, before collapsing to low values again.According to Eq. (4a), the baroclinicity builds up nearly linearly during the periods of low f and then collapses to a value below s 0 before the quasilinear build-up starts again.This type of behaviour resembles a relaxation oscillation, although in most models of relaxation oscillators explicit threshold properties are used.
We can derive several more analytical properties of the system in the various asymptotic regimes, but these may be considered less relevant for the application of this system, given the simplifications used to derive the original set of equations.Instead, we illustrate numerical solutions for various values of E in Figure 2. The figure illustrates the transition from a simple oscillation for low E to a relaxation oscillation for larger E. Also notable is the shift to a larger period for larger E, corresponding to the increasing time spent at low values of f .In other words, the nonlinearity decouples the time-scale between two bursts and the time-scale of the burst itself.
Despite the obvious simplifications in this model, the ensuing relaxation oscillation appears to have several relevant implications for the observed variance in storm tracks.We will discuss these in more detail in the next section and in the concluding section.

Observational evidence
We have analysed the variations in baroclinicity and heat flux on the upstream side of the North Atlantic storm track using the December-January-February (DJF) data of the ERA-40 reanalysis dataset, spanning 1957-2002.Climatological results are illustrated in Figure 3.As in Hoskins and Valdes (1990), the local baroclinicity at 775 hPa was calculated as where f = 2 sin φ and N = g d ln θ/dZ.The height derivatives of u and θ were calculated using a second-order centred finite-difference approximation between 700 and 850 hPa, using  geopotential height (Z) as the vertical coordinate.Hoskins and Valdes (1990) note that the baroclinicity at these levels is dominated by wind-shear variations rather than stability variations and we confirmed this to be true for the ERA-40 reanalysis date set as well.By the thermal wind relation, we then find that the above expression for baroclinicity is essentially the same as the expression used in Eq. ( 1).In order to analyze the typical baroclinicity related to the storm track, a sector spanning 30-50 • N and 30-80 • W was selected for spatial averaging to obtain a time series for baroclinicity, as illustrated in Figure 3.
The perturbations used to calculate the instantaneous transient heat flux between 700 and 925 hPa were computed by subtracting 10 day running means (which were calculated using a Lanczos filter) of instantaneous values for the vertically averaged v and T. This 10 day cut-off filter is suitable for representing highfrequency eddies (Athanasiadis and Ambaum, 2009;Lorenz and Hartmann, 2001).Again, a specific sector of the North Atlantic, spanning 30-60 • N and 30-80 • W, was selected for spatial averaging to obtain a time series for the eddy heat flux.
The high heat-flux regions in Figure 3 are not precisely colocated with the high baroclinicity regions, as also emphasized in Chang and Orlanski (1993).Indeed, we expect the heat-flux maxima to be located somewhat downstream of the baroclinicity maxima, although for the Atlantic region this offset is fairly small, possibly due to lee cyclogenesis over North America.Figure 4 shows a selection of time series for the winters 1995-1999.These were chosen for no particular reason other than that they reflect the type of variance within a winter season and variance across different years.It is clear from these time series that the heat flux tends to come in bursts, reminiscent of the nonlinear oscillator model at higher energies.The observed baroclinicity does not vary as drastically as the heat flux, but it can be seen that maxima of the heat flux tend to coincide with periods where the baroclinicity decreases.Additionally, the strictly periodic behaviour of the oscillator model is not really reflected in the observations, although periods of quasi-periodic behaviour appear to occur, for example, in the first half of the 1998-1999 winter season.
Figure 5 shows a composite (superposed epoch) centred at maxima of the heat flux, where only maxima larger than 30 K m s −1 were selected.This corresponds to an average of five events per winter season.The composite is only weakly dependent on the chosen threshold value.The composite heat flux shows a clear peak and a uniform, symmetric decay away from the central time.This is not an entirely conclusive test, as any signal will show such behaviour under composites whenever the temporal autocorrelation decays as it does in the figure.However, the time series in Figure 4 seem to indicate that the composite is representative of the behaviour of the observed bursts in the heat flux.The periodic behaviour of the oscillator model is not evident in the composite, indicative of a lack of periodicity or a weak quasi-periodicity in the data.However, the heat flux beyond 2 day lags does not decay to very low values in the composite, contrary to the quiescent periods in Figure 4.This indicates the presence of further bursts in heat flux at those lags in the full data set.
Figure 5 further shows that the composite excess baroclinicity (here calculated as the excess over a running mean which, according to our model, is expected to coincide with the mean linear dissipation rate s 0 ) is also consistent with the nonlinear oscillator model.The peak in composite heat flux coincides with a reduction in the composite baroclinicity.The maximum of decrease in baroclinicity is in fact lagged slightly ahead of the peak in heat flux, perhaps indicating that the eddy heat flux is not completely coincident with the eddy processes that consume the available potential energy.It is possible that other measures of eddy growth rate and eddy activity show a phase relation that is more aligned with that expected from the oscillator model.The build-up of baroclinicity between bursts of heat flux, as evident in the oscillator model, is not so clear in the composite because of the regression to the mean away from the central time.Furthermore, the interquartile range of the baroclinicity is large, so any trends will be hard to observe.Nonetheless, the time series in Figure 4 show anecdotal evidence of a build-up of baroclinicity ahead of a burst in heat flux and the fact that the composite baroclinicity shows a clear decay around the peak of the heat flux implies that, away from the peaks, the baroclinicity has to build up on average.
The lack of periodicity in the observations is evident.The waiting time distribution shows a broad distribution (perhaps consistent with Messori and Czaja (2013b)) with, however, a possible peak at short waiting times.The statistics are not strong enough to favour either possibility, without independent evidence for the underlying models.It is hard to make further quantitative links between our observations and the oscillator model.For example, choosing a lower threshold to define the heat flux maxima will result in more events and, correspondingly, a shift to shorter waiting times.In the context of our nonlinear oscillator model, this could be related to higher values of F or smaller values of E, or both.

Concluding remarks
The nonlinear oscillator model described in section 2 displays some salient features that can be seen in observations, here illustrated with data from the North Atlantic storm track.
Firstly, the heat flux does not appear to be uniform in time, or even uniformly random, but comes in bursts of activity.This aspect of the heat flux has been observed before (Swanson and Pierrehumbert, 1997;Messori and Czaja, 2013a).In the full three-dimensional atmosphere, this burst corresponds to the development of individual systems or periods of high storm activity.In the model, the time-scale between systems is set by the forcing of the mean temperature gradient F, which restores the baroclinicity after its collapse.However, the time-scale between individual bursts is non-trivially associated with the time-scale of the forcing because of the nonlinear nature of the system.Nonlinear oscillations operating at higher energies will have longer time-scales for the same forcing time-scale.
Secondly, the baroclinicity seems to build up to a maximum value before collapsing to a stable state in which the dissipative processes dominate the baroclinic instability and the environment is not conducive to system growth.This appears to correspond to the increasing northward tilt of the storm track before it returns to a more zonal regime, following a mechanism for storm-track tilt proposed by Orlanski (1998).This latter aspect is currently being investigated in detail and will be reported on in a separate article.Indications of this behaviour in time can be gleaned, for example, from Frame et al. (2011Frame et al. ( , 2013) ) and Franzke et al. (2011), where the North Atlantic eddy-driven jet is seen to tilt progressively northeastward before collapsing back to a more zonal direction.
Thirdly, on the time-scales of this storm track variability, there is a non-trivial, out-of-phase relationship between baroclinicity and heat flux.Although the spatial maxima of the two quantities have both been used to diagnose the locations of storm tracks, on synoptic time-scales they do not vary together.This reconciles the seemingly contradictory observation that high baroclinicity should be conducive to high heat fluxes, but at the same time high heat fluxes should lead to low baroclinicity.
There are some obvious caveats to our nonlinear oscillator model.The most serious is possibly the lack of a clear periodicity in the observed data.Indeed, recent evidence (Messori and Czaja, 2013b) indicates that meridional heat flux covers a very broad range of frequencies, from synoptic to planetary-wave scales.In response, we would argue that the strictly periodic behaviour of the model depends on the external forcing F being constant in time.Given that the forcing describes both radiative and kinematic effects, it would follow that the forcing in fact would vary over time.If the time-scale of this variation is longer than the time-scale of the forcing itself, the relevant behaviour of our model does not change, but the strict periodicity of the model will be broken.
There appears to be a strong analogy in the formation of vertical convection.In fact, Yano and Plant (2012) derive a formally equivalent set of equations based on mass-flux schemes (their Eqs (19a) and (19b)).In convective events, we generally observe a build-up of convective available potential energy (CAPE), which may then be released by a convective event (Yano et al., 2012).The typical interpretation would be that the convection acts on a shorter time-scale than the build-up of CAPE, which is implicit in the quasi-equilibrium assumption for convective parametrizations.Convective parametrizations contain triggers where a convective event is initiated when a certain threshold of instability is achieved.The present nonlinear oscillator model gives a different perspective, in that we do not need a trigger to achieve the same effect.The nonlinearity gives the impression of trigger behaviour: a sudden onset of a perturbation, which stabilizes the flow on a short time-scale.The nonlinearity also gives the impression of time-scale separation, where the time-scale of an individual perturbation is much shorter than the time-scale between perturbations.The difference between these time-scales is regulated by the energy parameter (the amplitude) in our oscillator model.If such a simple nonlinearity can be included in convective parametrizations, then we do not need to invoke the quasi-equilibrium assumption and we can remove switches from convective parametrizations.This has implications for the construction of adjoint numerical models.

Figure 2 .
Figure 2. Time series of heat flux f (solid lines), rescaled with F, and excess baroclinicity s − s 0 (dashed lines), rescaled with √ F. The time is rescaled with the natural frequency √ 2F of the system and each tick mark corresponds to one period of this natural frequency.The three panels correspond to E/2F = 0.1, E/2F = 1 and E/2F = 10.

Figure 3 .
Figure 3. Polar stereographic view of the baroclinicity (solid contours, displaying values of 0.5 and 0.6 day −1 ) and heat flux (dashed contours, displaying values of 10 and 20 K m s −1 ) averaged over the 1957-2002 DJF winters.The sector used for spatial averaging is outlined by the thin solid line.

Figure 4 .
Figure 4. Time series of heat flux (solid lines, left axes, units k m s −1 ) and baroclinicity (dashed lines, right axes, units day −1 ) for the DJF winter, for years 1995-1999 ((a)-(e)).The time is in days since 1 December of the indicated year.The heat flux and the baroclinicity are evaluated at the upstream side of the North Atlantic storm track, as indicated in Figure 3.

Figure 5 .
Figure 5. Composite of heat flux and baroclinicity for the winters of 1957-2001, centred around the maxima of the heat flux.The solid line is the median value of the heat flux and the dashed line is the median value of the baroclinicity.The shading corresponds to the interquartile range of each quantity.The anomalous (excess) baroclinicity has been plotted; the mean offset in the baroclinicity is 0.46 day −1 .