Articles | Volume 30, issue 1
Research article
07 Feb 2023
Research article |  | 07 Feb 2023

On the interaction of stochastic forcing and regime dynamics

Joshua Dorrington and Tim Palmer

Stochastic forcing can, sometimes, stabilise atmospheric regime dynamics, increasing their persistence. This counter-intuitive effect has been observed in geophysical models of varying complexity, and here we investigate the mechanisms underlying stochastic regime dynamics in a conceptual model. We use a six-mode truncation of a barotropic β-plane model, featuring transitions between analogues of zonal and blocked flow conditions, and identify mechanisms similar to those seen previously in work on low-dimensional random maps. Namely, we show that a geometric mechanism, here relating to monotonic instability growth, allows for asymmetric action of symmetric perturbations on regime lifetime and that random scattering can “trap” the flow in more stable regions of phase space. We comment on the implications for understanding more complex atmospheric systems.

1 Introduction

Rossby (1940) first suggested the atmosphere may possess multiple large-scale quasi-equilibrium states to describe the empirically observed tendency for European weather patterns to take on a finite range of frequently recurrent configurations. These are now known in meteorology as atmospheric regimes, and in dynamical systems, such behaviour is known as crisis-induced intermittency (Ott2002). The first theoretical study of multiple atmospheric equilibria, carried out by Charney and DeVore (1979) (hereafter CdV79), showed that the flow in a low-order mid-latitude barotropic channel could equilibrate to either a zonally symmetric flow or a stationary blocking pattern. They suggested real-world equivalents of these equilibria states might be destabilised by perturbations from unresolved modes (in this idealised model, this included spatial scales of ∼𝒪(1000) km), causing shifts from one state to another and introducing a vacillation between large-scale weather patterns.

Similar results were obtained in other simple models (Charney and Straus1980; Kallen1981; De Swart1988), but these equilibria were often absent in less truncated systems, rendered “unrealisable” by shifts in the model attractor or collapsing completely due to structural instability (Reinhold and Pierrehumbert1982; Cehelsky and Tung1987). Therefore later work moved away from analysis of fixed points, treating regimes as chaotic flows originating from attractor merging (Itoh and Kimoto1996, 1997) and exploring the role of homoclinic orbits in preferred transitions between regimes (Kimoto and Ghil1993; Kondrashov et al.2004; Selten and Branstator2004).

Despite a deepening mathematical understanding of regimes, the Earth system appears to be high-dimensional, and as a result linking the regime behaviour of simple models to that of realistic models has proved challenging, especially when further complicated by the introduction of stochastic forcing. Earth-system models often contain explicit stochastic parameterisation schemes and in addition a range of small-scale processes that can be thought of as randomly forcing the large-scale flow. Contrary to linear intuition, there is evidence that this fast variability can stabilise regimes: Kwasniok (2014) (hereafter K14) found additive noise increases regime lifetime in both the Lorenz 1963 system (Lorenz1963) and a chaotic reformulation of the CdV79 model found by Crommelin et al. (2004) (hereafter C04). Further, stochastic forcing has been shown to improve regime predictability by increasing persistence in truncations of the Lorenz 1996 model (Düben et al.2014; Christensen et al.2015) and climate model regime representation (Dawson and Palmer2015), as well as the representation of other slow-varying climate modes (Palmer and Weisheimer2011; Berner et al.2012).

Outside of geoscience, similar stochastic persistence effects have been observed and explained for transient chaotic maps. For such open systems, chaotic dynamics only persist for finite time, eventually “escaping”, often to a steady state (Lai and Tél2011a). It has been shown for such systems how stochastic forcing can delay the escape of trajectories in the case of 1D (Reimann1996; Lai and Tél2011b) and 2D (Altmann and Endler2010) maps via a “fattening” of the fractal structure. This is of interest for fully chaotic geophysical regime systems as the transition between regimes can be conceived of as an escape from one regime to another, and so increases in regime persistence correspond to decreases in escape rate. This approach was introduced in Franaszek and Fronzoni (1994) to understand the stochastically forced Duffing oscillator.

We will study the stochastically forced CdV79 model, as retuned by C04, and show that similar mechanisms explain stochastic persistence in this six-dimensional, continuous-time, non-hyperbolic regime system as in simpler transiently chaotic maps. We introduce the model and its regime structure in Sect. 2. We then discuss the impact of stochasticity on that regime structure in Sect. 3 and demonstrate that this can be understood as a result of the theory of hitting time problems in combination with a geometric mechanism. We summarise our results and discuss the implications for weather and climate in Sect. 4.

2 Model formulation

The deterministic Charney–deVore model is obtained by spectral truncation of a barotropic flow in a β-plane channel geometry. It can be expressed as a system of six ordinary differential equations, derived in full in Appendix A:

(1) x ˙ 1 = - C ( x 1 - x 1 * ) + γ ̃ 1 x 3 x ˙ 2 = - C x 2 + β 1 x 3 - α 1 x 1 x 3 - δ 1 x 4 x 6 x ˙ 3 = - C x 3 - β 1 x 2 + α 1 x 1 x 2 + δ 1 x 4 x 5 - γ 1 x 1 x ˙ 4 = - C ( x 4 - x 4 * ) + ϵ ( x 2 x 6 - x 3 x 5 ) + γ ̃ 2 x 6 x ˙ 5 = - C x 5 + β 2 x 6 - α 2 x 1 x 6 - δ 2 x 4 x 3 x ˙ 6 = - C x 6 - β 2 x 5 + α 2 x 1 x 5 + δ 2 x 4 x 2 - γ 2 x 4 .

The six variables {xi} represent spectral mode amplitudes of a streamfunction field, subject to a linear relaxation to a zonally symmetric background state, Coriolis and orographic forces, and quadratic advective non-linearities. To explore stochastic effects, we introduce an additive white noise vector to the tendency, given by a normal distribution ξi with standard deviation σ, as in K14:

(2) x ˙ i = x ˙ i , deterministic + σ ξ i .

The system is integrated using the Euler–Maruyama difference scheme with a time step of dt=2×10-4 model time units (MTU) and is sampled every 1 MTU. We compute integrations of length 200 000 MTU for a range of σ values.

Figure 1(a) Projection of the deterministic attractor into the space of the leading EOFs, capturing 92 % of the total variance. (b, c) The flow patterns corresponding to the first and second EOFs of the system. Blue and red dots show the positions of fixed points corresponding to blocked and zonal flows respectively.


Figure 2A 2000 MTU integration of the CdV79 system, showing the evolution of each mode amplitude. The dynamics are chaotic but weakly so, showing clear quasi-periodic behaviour. Two regimes of behaviour can be seen: fast chaotic oscillations, divided by periods of slow evolution with irregular duration. The mode values at the blocked and zonal fixed points are shown with dashed blue and red lines.


To visualise the model phase space, we project the six-dimensional attractor into the space of the two leading empirical orthogonal functions (EOFs) of the deterministic system, as shown in Fig. 1, capturing 94 % of the system's variability (Fig. S1 in the Supplement shows the various projections of the six x modes). The attractor is highly structured, with a number of strongly preferred trajectories through phase space (note the logarithmic colour scheme in Fig. 1). This low-dimensional attractor shows no obviously separable phase space regions such as the twin lobes of the Lorenz 1963 system, and we see that the trajectories are highly twisted together, producing a scroll-like geometry.

A representative sample of the model's temporal evolution is shown in Fig. 2. Intermittency is readily apparent, with predictable slowly evolving periods, with near-zero phase speed, corresponding to persistent “blocked” flow patterns, divided by highly chaotic, fast-evolving flows that contain both zonally symmetric states and transient wave activity. The regime structure is therefore quite asymmetrical, in contrast to the Lorenz 1963 or 1996 systems whose regimes are essentially equivalent to each other. Instabilities in the CdV system arise from two basic mechanisms: an orographically induced saddle-node bifurcation (a transition from one to multiple fixed-point equilibria), as identified in the original CdV79 paper, and a barotropic Hopf bifurcation (a transition from a fixed point to a periodic orbit) associated with the onset of a zonally propagating waves in the model domain. Both bifurcations are codimension 1, existing along a surface with one fewer dimension than the parameter space, and so can in general be found by tuning one model parameter. Stable steady-state and periodic wave solutions to Eq. (1) are common for arbitrary sets of model parameters (not shown), while chaotic – and especially intermittently chaotic – solutions are harder to find. In the original formulation of CdV79, two stable fixed points were found, corresponding to a zonally symmetric state, and a “blocked” flow, with strong wave activity. The model parameters used here, and in C04 and K14, are close to a parameter set where the barotropic Hopf and orographic saddle-node bifurcations merge: one fixed point bifurcates into two fixed points and a limit cycle. This fold-Hopf point has codimension 2, meaning two model parameters must be tuned to obtain it. High-codimension bifurcations seem to be of reduced relevance to the real world, representing exceptional states in parameter space. However, the fold-Hopf bifurcation is known to produce a range of lower codimension and so more relevant bifurcations at nearby parameter values, including paths to chaotic dynamics (Champneys and Kirk2004). Most relevant to us is the generation of homoclinic orbits, a type of bifurcation where orbits stretch from a fixed point back to itself, with an infinite period. These orbits can lead to fully fledged chaotic dynamics and the creation of a strange attractor through intersection with an unstable toroidal limit cycle. At the boundary crisis this creates a “whirlpool repeller” (Shilnikov et al.1995) and exactly the scrolling structure seen in Fig. 1. The two fixed points, analogues of those in CdV79, remain but are now unstable, and at the parameter values chosen here the strange attractor moves between the vicinity of these two fixed points, as seen in Figs. 1 and 2, shadowing collapsed homoclinic orbits. The flow patterns corresponding to these fixed points are shown in Fig. S2. See De Swart (1988) and C04 for more details on the bifurcation structure of the model. As a result the regimes of the CdV system are rooted in the dynamics of unstable homoclinic orbits, with close approaches to the unstable fixed point leading to periods of slow evolution and then rapid “bursting” to chaotic, more turbulent flow. Such behaviour is seen in other simple atmospheric models (Crommelin2002) including the Lorenz 1984 model (Shilnikov et al.1995) and in the dynamics of non-linear circuits (Pusuluri and Shilnikov2018) and neurons (Izhikevich2006).

This type of regime behaviour is of a fundamentally different type to that of the archetypal Lorenz 1963 system. In Lorenz 1963, crisis-induced intermittency creates regime dynamics when the basins of attraction of two attractors merge together, triggering movement between two distinct regimes of fully chaotic behaviour. The Charney–deVore system, in contrast, experiences almost deterministic dynamics before drifting away and entering a transient chaotic regime. At an unpredictable later time, the system will once again return to the neighbourhood of the periodic orbit and re-enter the slow-evolving regime. For a recent discussion of the co-existence of predictable and chaotic flow states in atmospheric systems see Shen et al. (2021).

Regimes in CdV

When studying a multimodal system, bulk metrics such as the mean or variance can obscure the underlying dynamical structure. Therefore it is natural to use a regime framework as a way of partitioning phase space into dynamically distinct regions of interest. Following the approach of K14, we use a hidden Markov model (HMM) to identify qualitatively different dynamical regimes, with the minor methodological change of fitting the HMM to the full six-dimensional state vector, rather than restricting it to the x1x4 subspace as in K14. As for most clustering methods, we must make a subjective choice of regime number when fitting the HMM. We use three regimes, which were found in K14 to separate the qualitatively distinct model dynamics well, capturing both the persistent quasi-predictable blocking regime and the highly chaotic zonal regime and the transitions between them. Figure 3a shows the regions of the deterministic attractor assigned to each regime, and Fig. S3 shows the corresponding flow patterns. The blocking regime is small in phase space extent, focused on the tightly spiralled trajectories that start close to the blocking fixed point. It clearly separates from the swirling trajectories assigned to the transitional regime. The zonal state contains the very chaotic “branching” of the attractor into different orbits and contains the early parts of these orbits until their transition back to blocking becomes unavoidable. Interestingly, these regimes broadly correspond to regions of the attractor with different numbers of unstable dimensions (Fig. 3b), supporting our choice of three regimes as dynamically appropriate. Our regime assignment also has some similarities to the regimes found in Strommen et al. (2022), which used persistent homology to identify regimes in the system studied here and found a persistent component similar to our blocking state.

Figure 3(a) The assignment of points on the deterministic model attractor to the three hidden Markov model regimes projected into the EOF1–EOF2 plane. (b) The number of unstable dimensions as identified by the positive local Lyapunov exponents, calculated at each point along the deterministic model attractor, projected into the EOF1–EOF2 plane.


Figure 4Left: time evolution of the x1 mode in deterministic and stochastic cases, with points coloured by their regime assignment. Right: the corresponding phase velocities.


For non-zero noise amplitude, σ≠0, we identify regimes in the same way, refitting the HMM for each σ value. This is necessary as while the qualitative shape of the attractor1 is robust over a range of noise amplitudes, the fine structure seen in Fig. 1 breaks down, and the stochastic simulations experience shifts of probability density closer to the location of the blocking fixed point (shown in Fig. S4).

Figure 4 shows the temporal evolution of the zonally symmetric x1 mode, and the phase velocity vt=12i=16(xi,t+1-xi,t-1)2, for different σ values, coloured by regime assignment. The characteristics of the regimes are broadly preserved even up to σ=1.6×10-2: the persistent blocking state has low x1 values and very low phase velocities, and the zonal state has high x1 values and relatively low phase velocities, while the transitional regime has high x1 values and high phase velocity. However, for the larger value of σ=3.0×10-2, it is apparent that there is no meaningful regime structure remaining: the phase velocity never drops below 0.2, indicating a loss of quasi-stationarity. At this point the deterministic dynamics play second fiddle to the stochastic term, and the assignment of regimes becomes meaningless.

3 Stochastic persistence

Figure 5a shows Fourier transforms of the total mode amplitude ixi2 for three representative amplitudes of stochastic forcing. The deterministic σ=0 spectrum possesses a clear peak consistent with the quasi-periodic nature of the dynamics. For the σ=0.01 case this peak has collapsed, producing a flatter spectrum with power shifted towards low frequencies. This is not a monotonic effect: for σ=0.04, the total power in the low-frequency region decreases, and the spectrum flattens. Figure 5b shows this non-monotonicity clearly, with a peak in low-frequency variability close to σ=0.01.

Figure 5(a) The power spectral density of the total mode amplitude, ixi2, for three different strengths of stochastic forcing, σ. (b) The spectral power averaged over frequencies from 1/500 to 1/10000 MTU−1 (corresponding to the shaded region of panel a) as a function of σ.


Figure 6The distribution of blocking regime lifetimes, as a function of stochastic forcing amplitude. Shaded regions mark the interquartile range, with the thick black line showing the median and the red line the mean. Dotted black lines show the 0th and 95th percentiles of the distribution.


This is in excellent agreement with K14, which showed increasing regime persistence with noise amplitude that peaked at σ=0.008 and then decreased. We show here that this is associated with a genuine slowing of the system variables visible in the Fourier spectrum and is observable independent of adopting a regime perspective. With this reassurance, and given the HMM maintains the same qualitative characteristics for σ values ≤0.016, we judge that the regime assignment provides a sensible way to analyse the amplified low-frequency variability seen in Fig. 5 in more detail. Figure 6 shows the distribution of lifetimes for each of the three regimes, as a function of increasing σ. While all regimes show some increase in regime persistence for moderate noise amplitudes, the effect is by far the most significant in the blocking regime, showing a near trebling of mean regime lifetime that peaks at σ=0.008. As shown in Fig. S5, this peak is associated with a 150 % increase in the number of blocks lasting at least 100 MTU and a 600 % increase in those lasting longer than 175 MTU. Changes in the median blocking lifetime are smaller, and the distribution is clearly skewed, with the 95th percentile of blocking lifetime exceeding 600 MTU when persistence is maximum.

Figure 7(a) The mean (solid lines) and median (dashed lines) Euclidean distance between points in the stochastic simulation and their closest neighbour on the deterministic attractor, shown for the blocking and non-blocking regimes. (b, c) As in (a) but showing the distance between points in each simulation and the blocking (b) and zonal (c) equilibria respectively.


Figure 8A schematic showing one variable, x, of a chaotic process, which moves between periods of chaotic variability and periods of predictable, monotonic growth (shaded in grey). The predictable regime is entered at a random position xi, in the neighbourhood of an unstable fixed point x0, and once it reaches an amplitude xc, it will re-enter the chaotic regime. At a position xp, we consider the impact of applying a random perturbation Δx to the x variable, drawn from some symmetric probability distribution, with vanishing support outside the interval [x0xp,xcxp]. We see that, due to the convex growth of x, a perturbation with negative sign will increase the duration of the predictable period more than a corresponding perturbation with positive sign will shorten it.


How can the ability of random forcing to stabilise regime dynamics be understood? What causes the non-monotonicity between stochastic amplitude and persistence, and why is the response here asymmetric, primarily in the blocking regime? As mentioned, non-monotonic changes of transient lifetime due to stochastic forcing are understood for low-dimensional transiently chaotic maps. Symmetric perturbations can decrease on average the region of phase space in which chaotic trajectories can “escape” (i.e. decay) for a large class of 1D concave maps, lengthening their lifetime (Reimann1996; Lai and Tél2011b). Additionally, stochasticity smooths the probability density of the deterministic system. If the escape region of the system is of higher than average density due to low dimensionality, then stochastic forcing can on average knock more trajectories out of this region than in to it, again increasing persistence (Altmann and Endler2010; Faisst and Eckhardt2003).

In our system we first consider the mean state changes. The average distance between points in the stochastic simulations and their closest neighbour on the deterministic attractor increases, approximately linearly, as the amplitude of the stochastic forcing increases (Fig. 7a). Once again there is asymmetry between regimes, with the minimum distance during blocking growing far faster than during the other regimes. Figure 7b and c show that this is associated with a drift towards the neighbourhood of the unstable blocking fixed point, with no equivalent result seen for the zonal fixed point (also visible in Fig. S4). If we consider the dynamics in the neighbourhood of the blocking fixed point, we realise the deterministic tendency there is near zero, and so the additive stochastic terms dominate entirely: the dynamics approach six-dimensional Brownian motion. How long, in expectation, will it take for this Brownian motion to move the system away from the fixed point once again and return it to a dynamics-dominated flow state? This question is simply a reformulation of the classic hitting time problem: at what time will a random process first reach a distance a from its initial condition? For Brownian motion in three or more dimensions it is known that in fact the average time taken is infinite: the distribution of hitting times is fat-tailed, and so the mean diverges (Kallenberg2017). Therefore the skewed distribution of blocking lifetimes seen in Fig. 6 can be understood in terms of the amount of time spent near the fixed point, as shown in Fig. 7b. As σ increases past 0.008, the average distance from both the attractor and the blocking fixed point increases. As the stochastic forcing becomes stronger than the deterministic damping along the stable manifold, the system increasingly explores the full six-dimensional phase space, and so the directionality of the forcing vanishes. There is a conceptual similarity here to the second stochastic persistence mechanism of Altmann and Endler (2010) found for mixed phase space systems, whereby trajectories in a 1D transiently chaotic map can scatter into islands of stability, which introduces a fat tail to the transient lifetime. However, the increased persistence we find is not entirely a result of excursions to the noise-dominated neighbourhood of the fixed point. There is also an important interaction between stochasticity and the low-dimensional, ordered flow of the blocking regime. We present a heuristic mechanism similar in character to that of Lai and Tél (2011b) but now considering a 1D flow as in the schematic shown in Fig. 8. We can think of this as a parameterisation along the 1D unstable manifold of our CdV79 system moving between a highly chaotic state and a slowly evolving, predictable state. Immediately upon entering the regime of slow-evolving dynamics, the system is at a point x=xi and will stay in the quasi-stationary regime until it reaches x=xc. We assume that x grows monotonically with the time t in this regime, so we can parameterise x in terms of t or vice versa: f:xt, g:tx. If we let x=f(t) be any convex, monotonic function defined over the interval [x0,xc], then t=g(x)=f-1(x), and so g is therefore a concave, monotonic function with negative curvature. We consider applying a single perturbation Δx to the system at a point xp, drawn from an even probability distribution Px) with fast-decaying tails and support over [x0xp, xcxp]. Informally, this amounts to choosing a point xp far from the edge of the predictable regime, and so there is only a vanishing probability of Δx directly triggering a regime transition. When the standard deviation of P is small compared to xcx0, this will apply to most values of x in the interval. From the evenness of P, we see that the impact of the perturbation on the variable x, in expectation, is zero:

(3) Δ x = 0 .

However this is not in general the case for 〈Δt:

(4) Δ t = g ( x p + Δ x ) - g ( x p ) = G p ( Δ x ) ,

where Gp, as a translation of g, is also concave. Using the definition of the expectation we find

(5) Δ t = - G p ( Δ x ) P ( Δ x ) d Δ x ,

where the limits are justified by our requirement for vanishing tails for P and our choice of xp.

Because Px) is even, then if we Taylor-expand Gp, we will lose all odd powers on integration and be left with

(6) Δ t = - n = 1 1 ( 2 n ) ! d 2 n G p d ( Δ x ) 2 n P ( Δ x ) d Δ x .

And having taken Δx to be small, we can drop high-order terms:

(7) Δ t = - 1 2 d 2 G p d ( Δ x ) 2 P ( Δ x ) d Δ x .

Because Gp is concave, this derivative is negative, and we obtain

(8) Δ t = - - 1 2 d 2 G p d ( Δ x ) 2 P ( Δ x ) d Δ x = - I ,

where I is positive. Therefore any convex perturbation growth will be delayed (i.e. intermittency persisted) by small symmetric random perturbations, and the average lifetime of the stable regime will increase. The case for a negatively oriented variable with concave curvature is precisely equivalent, while for a positively oriented variable with concave curvature or a negatively oriented variable with convex curvature, decreased persistence would be seen. For the linear case with zero curvature, there is of course no impact of the perturbation on persistence.

When considering the particular case of f(t)=x0eλt, and P(Δx)=N(0,σ), which is interesting for our system, we obtain analytically Δt=-σ22λxp2 (see Appendix B). This quadratic dependence on σ2 is exactly what is found in Lai and Tél (2011b) and Reimann (1996) for the escape rate for 1D maps. However it is not clear that the two derivations are totally equivalent; we consider a convex monotonic instability and they a concave non-invertible map. Higher σ corresponds to longer persistence (remembering our overall weak noise approximation), as does small λ, indicating that nearly stable intermittent dynamics can be most substantially persisted. We also see in this case that smaller xp leads to larger drift, indicating that perturbations soon after the system has entered the intermittent state are most impactful. While we have considered only a single perturbation, the reasoning generalises to a continuous stochastic process, where the 〈Δt term we have calculated is now a spatially dependent drift term. The crux of the argument can be understood solely from examination of Fig. 8. Any low-dimensional, convexly growing instability in the climate system should be susceptible to the same stochastic suppression. We can also see why only the blocking regime showed dramatic persistence increases, as it is the only regime where the variables evolve with a constant sign of curvature and where there is only one unstable dimension. The decrease of lifetime for larger noise values does not emerge naturally from this framework but can be understood as entering a noise-dominated regime where the stochastic terms are much larger than the deterministic tendency.

4 Conclusions

In this paper we have investigated the interaction between stochasticity and strongly non-linear dynamics in a six-equation conceptual model of atmospheric flow. We consider the case of additive stochastic forcing applied to a chaotic flow which transitions intermittently between chaotic zonal and quasi-stationary blocked regimes, which Kwasniok (2014) showed exhibited stochastically reinforced regime persistence for a range of noise amplitudes. We have investigated the mechanisms driving this behaviour and provide heuristic arguments based on the stochastic suppression of instability growth and random scattering close to an unstable fixed point. The low dimensionality and ordered nature of the blocking regime leads to a slow monotonically growing instability, which will on average be persisted by a random perturbation rather than accelerated. This effect is amplified, and a heavy tail introduced into the lifetime of blocking events, through chance excursions to the neighbourhood of a deterministically inaccessible fixed point, where the overall deterministic tendency is negligible. Here, the dynamics are almost entirely noise-dominated, and the system will remain close to the fixed point until it scatters away. Ultimately, if stochastic forcing is increased past a certain point, persistence decreases. This is not a resonance effect but is instead a distinct feature exhibited by transient systems and systems near boundary crises (Altmann and Endler2010), where a balance must be struck between providing enough diffusion to escape the deterministic transition trajectories and providing so much diffusion that new purely stochastic routes out of the regime are created. The mechanisms driving stochastic regime dynamics have many similarities to the dynamics of stochastic transiently chaotic maps. Understanding regime transitions as escapes from the neighbourhood of one non-attracting set to the neighbourhood of another, which will itself eventually escape, provides a different perspective that makes the connection to the theory of stochastically perturbed escape rates clear. In atmospheric physics more attention is generally paid to continuous-time fully chaotic conceptual models such as the archetypal Lorenz 1963 system (Lorenz1963), but the study of regime systems may well benefit from a deeper consideration of the framework of transient chaos and the wealth of work done on discrete-time maps. Indeed, all chaotic systems can be understood in terms of transient shadowing of unstable periodic orbits (Cvitanović et al.2016), and such ideas have been recently applied to blocking flows (Lucarini and Gritsun2019) and regime systems (Maiocchi et al.2022). Characterising the unstable periodic orbits associated with each regime could even allow the non-monotonic dependence of regime lifetime on stochasticity to be analytically computed using perturbative methods (Cvitanović et al.1999). Understanding the impacts of stochastic parameterisation on truly transient processes, such as mesoscale convective features, is also of great importance (Pickl et al.2022). It remains to be seen however how relevant these simple mechanisms are in more complex atmospheric models, and indeed the impact of stochastic physics on climate model regime structure is still unclear, complicated by high sampling variability in the Euro-Atlantic region (Dorrington et al.2022). A natural next step is to ask how stochasticity impacts regimes in quasi-geostrophic settings and to explore how the mechanisms uncovered here evolve (or, indeed, are rendered extinct) as the complexity of the flow increases. Finally, intermittent dynamical systems such as those possessing circulation regimes are typically close to a bifurcation. Such near-crisis dynamics are particularly sensitive to stochastic forcing, and as we have seen, the response of the system to forcing amplitude can be complex. As a further example, Wackerbauer and Kobayashi (2007) demonstrated high sensitivity of the escape rate in a spatio-temporal system to the spatial correlation of the noise process. This highlights the potential importance of carefully selecting the amplitude and kind of stochastic parameterisations in complex models: not all representations of stochasticity in a system can be expected to be equivalently effective.

Appendix A: Derivation of the Charney–deVore 1979 model

Here we derived the CdV79 model introduced in Sect. 2 from the shallow water equations, marking simplifying approximations in boldface. The shallow water system is a simple but highly useful model of fluid flow. Consider a single layer of fluid, with constant density ρ, negligible viscosity, and a vertical pressure gradient deriving solely from gravity (i.e. the fluid is in hydrostatic equilibrium) and which, as the name implies, has a representative horizontal length scale L much larger than the vertical length scale H. The dynamics of this system are completely determined by the momentum equation and the mass-conservation equation.

The momentum equation in a rotating frame is given by

(A1) v t + ( v ) v + f × v = - g H ,

while the mass-conservation equation is

(A2) D H D t + H v = 0 ,

where v is velocity, f represents Coriolis forces, g is simply the gravitational constant (assumed constant), and H is the depth of the fluid column. Here the mass conservation simply expresses that if velocities converge at a point, the height of the fluid column must rise in response. Application of the identity (ab)c=-a×(b×c)+(ac)b to Eq. (A1) gives us

(A3) v t + ( × v ) × v + 1 2 v 2 + f × v = - g H ,

and by introducing the absolute vorticity, ξ=×v, and taking the curl of both sides, we get the shallow water vorticity equation:

(A4) ξ t + × ( ( ξ + f ) × v ) = 0 .

Again, we apply the triple product identity and the chain rule and find

(A5) ξ t + ( v ) ( ξ + f ) + ( ξ + f ) ( v ) = 0 ,

where ξ=(×v)=0 by definition and (ξ)v=0 because ξ, as a curl of horizontal wind, does not project onto the variability of v. Indeed, at this point we note that both ξ and f have no horizontal components, and so we focus only on the scalar equation for the z component, letting ξ:=ξz and f:=fz for convenience. The divergent velocity term can be rewritten with the aid of Eq. (A2) by rearranging

(A6) v = - 1 H D H D t

and substituting it into Eq. (A5). This leaves us with

(A7) D ( ξ + f ) D t - ( ξ + f ) H D H D t = 0 ,

which can be put into a clearer form using the chain rule:

(A8) D D t ( ξ + f ) H = 0 .

That is, the potential vorticity q=ξ+fH is materially conserved. If the thickness of the fluid increases, then the absolute vorticity must likewise increase, due to the stretching of the flow.

We divide our fluid height H into some average fluid depth D and an anomaly component h(x,y) describing the height of surface orography, so that H=D-h(x,y). We assume shallow orography, so we may take the limit where hD1 and Taylor-expand H, truncating at first order so that 1H1D+hD2:

(A9) D D t ξ + f + ( ξ + f ) h ( x , y ) D = 0 .

In the mid-latitude atmosphere, this is a reasonable approximation at the largest scales. For example, if we let D be the height of the tropopause – around 8 km – and we let h be the average elevation of Europe or North America – between 500 and 800 m – then the limit holds. We also choose to place ourselves in the quasi-geostrophic limit of low Rossby number, where ξ<f. Finally, we take a β-plane approximation to the Coriolis force, and this brings us to the barotropic vorticity equation in the presence of shallow orography:

(A10) D D t ξ + β y + f 0 D h ( x , y ) = 0 = D q D t .

We define γ=f0D and introduce the streamfunction Ψ, given by u=-Ψy, v=Ψx, and ξ=∇2Ψ. We also generalise slightly by introducing an external source (or sink) of vorticity F. Writing our material derivative out explicitly, we now have

(A11) 2 Ψ t - Ψ y x [ 2 Ψ + γ h + β y ] + Ψ x y [ 2 Ψ + γ h + β y ] = F

or more concisely

(A12) 2 Ψ t + J ( Ψ , 2 Ψ + γ h ) - β Ψ x = F ,

where the Jacobian operator is given by

(A13) J ( A , B ) := A x B y - B x A y .

Taking our forcing to be a linear relaxation towards a background state, we now have our explicit equation for the evolution of Ψ:

(A14) t 2 Ψ = - J ( Ψ , 2 Ψ + γ h ) - β Ψ x - C ( Ψ - Ψ * ) .

This is the CdV79 equation, before spectral truncation, where Ψ* is our background streamfunction. We consider the flow in a channel geometry, with periodic boundary conditions in the zonal direction and with no-penetration and no-slip conditions in the meridional direction. Of course the real atmosphere has no hard meridional boundaries; however the impacts of the subtropical and mid-latitude jets have been shown to produce a waveguiding effect that meridionally confines low-frequency variability in a similar way, both in low-order models (Hoskins and Ambrizzi1993; Yang et al.1997) and in observations (Branstator2002), which makes this a reasonable first approximation.

We consider a channel of dimensionless zonal length x=[0,2π] and meridional length y=[0,bπ], so that b then is a free parameter controlling the shape of the channel.

Discretising the model

To convert this infinite dimensional partial differential equation into a system of ordinary differential equations, we can use a Galerkin projection. This involves expanding our spatially dependent variables (Ψ, Ψ*, and h) in terms of an infinite linear sum of complete spatial basis functions, labelled {Φ(x,y)i}:


We start by defining this basis abstractly, requiring the basis functions to be eigenfunctions of the Laplacian operator and the zonal derivative:


The Jacobian term is the thorniest point in the spectral analysis, as it non-linearly mixes spectral modes together. It is the impact of this term that prevents an arbitrary spectral truncation of Eq. (A14) being a solution of the full equation. In full generality, we simply consider that the Jacobian of each pair of basis functions has some projection on every other basis function:

(A20) J ( Φ j , Φ k ) = c i j k Φ i .

The projection coefficient of the jth and kth modes on the ith mode, cijk, is computable simply through a direct pointwise product of the ith mode and the Jacobian, integrated over the domain and appropriately normalised:

(A21) c i j k = 0 2 π 0 b π Φ i J ( Φ j , Φ k ) d y d x 0 2 π 0 b π d y d x .

Under these constraints, we can now rewrite Eq. (A14) projected onto the ith basis function as

(A22) t Ψ i = β a i - 2 b i Ψ i - C ( Ψ i - Ψ i * ) - j , k a i - 2 c i j k ( a k 2 Ψ j Ψ k - γ Ψ j h k ) ,

where the double sum indexes over all basis functions. We make an explicit choice of basis functions now, taking a Fourier basis defined as

2cos(nyb)if m=02eimxsin(nyb)else,

where the meridional mode n>0 and the zonal mode m can take any non-integer value. It is easy to verify that these basis functions satisfy the boundary conditions and eigenfunction constraints defined above. To actually solve these spectrally transformed equations, it will be necessary to introduce a truncation, disregarding wave modes exceeding a cutoff wavenumber. Here we choose to introduce a zonal threshold of |m|1 and a meridional threshold of n≤2, restricting us to just six basis functions, as in CdV79, C04, and K14.

Figure A1The real-valued basis functions used in the spectrally truncated CdV79 system, plotted as streamfunctions.


Although up to this point our approximations have been motivated well physically, we cannot rigorously justify introducing such a drastic truncation without resorting to self-interest: such a system is much easier to study and can be understood in more detail than a high-dimensional model. Nonetheless, we may reason that any phenomena observed in this truncated model, which will have to derive solely from the dynamics of the largest scales, may have analogues in the true atmospheric circulation, although there they will be shaped and influenced by coupling to smaller scales. The resulting basis functions are shown in Fig. A1.

Concretely, while the true spectral equation for the largest wave modes, Eq. (A22), contains an infinite number of non-linear coupling terms, after truncation we are reduced to only 36 such terms, representing only interactions between the resolved modes: the ability for small scales to feed vorticity into the untruncated large scales is neglected.

We choose our linear relaxation to be zonally symmetric, representing the forcing effect of the pole-to-equator temperature gradient. We choose our orography to include a single peak straddling the periodic boundary, with a trough in the centre of the domain. h(x,y)=γ2cos(x)sinyb (this is shown in Fig. A2).

Figure A2The orographic profile used in the CdV79 system, with a trough in the centre of the domain and a ridge over the periodic boundary. Sub-panels show cross sections across the domain, following the dotted black lines.


Figure A3The background state in the Crommelin 1904 formulation of the CdV79 system. It is zonally symmetric and represents the impact of the meridional temperature gradient. Again, sub-panels show cross sections across the domain, following the dotted black lines.


Obtaining the coefficients in Eq. (A22) is now simply a matter of calculating integrals but is very involved and is best performed using symbolic computation software (an implementation using SymPy (Meurer et al.2017) can be found at, last access: 3 February 2023). After computing coefficients, we now diverge slightly from CdV79, transforming our variables to a real basis, as in C04, given by


This transformation is purely cosmetic; the behaviour of the system is unchanged. At the end of this process we obtain an explicit six-equation system:

(A23) ,

where colours encode the linear relaxation (green), Coriolis forces (blue), orographic impacts (orange), and non-linear advection (red). We see that there are far fewer than the 36 promised non-linear terms per equation: most quadratic interactions integrated over the domain vanish due to symmetry considerations. {αn},{βn},{δn},{γn},{γñ}, and ϵ are simple geometric constants of 𝒪(1) given in C04.

There are now a number of dimensionless free parameters in the model, which we set equal to the values used in K14 and C04:

  • xi* is the forcing vector, representing thermal relaxation towards a background state. We use a zonally symmetric forcing profile with x1=0.95 and x4=-0.76095; all other terms are zero (shown in Fig. A3).

  • C is the thermal relaxation timescale, set to C=0.1, describing a damping time of ∼10 d.

  • γ is the orographic amplitude, set to γ=0.2, corresponding to a 200 m amplitude.

  • β is the Coriolis parameter, set to β=1.25, defining a central latitude of 45.

  • b defines the channel half-width, set to b=0.5, which gives a channel of 6300 km × 1600 km.

Appendix B: Impact of a Gaussian perturbation on exponential growth

In Sect. 3 we found that

(B1) Δ t = - - 1 2 d 2 G d ( Δ x ) 2 P ( Δ x ) d Δ x = - I .

If we consider x to be undergoing exponential growth and subject it to a Gaussian perturbation so that x=f(t)=x0eλt, and P(Δx)=N(0,σ), then we can obtain an analytical form for Δt. Remembering that Gp(Δx)=g(xp+Δx)-g(xp), and seeing that g(x)=1λlogxx0, we have

(B2) Δ t = x 0 x c 1 λ log x x p 1 σ 2 π e - 1 2 x - x p σ 2 d x .

We transform our variables to w=Δxσ=x-xpσ, which implies xxp=xp+σwxp=1+σxpw and dx=σdw. Our limits of integration become w0=x0-xpσ<0 and wc=xc-xpσ>0.

If, as in the main text, we assume σ is small relative to a typical length scale of x (that is the perturbation is vanishingly unlikely to force the system into the chaotic regime directly), then these limits will be very large in absolute magnitude. Given this and the exponentially decaying tails of the Gaussian distribution, we approximate w0-, wc+. This gives us

(B3) Δ t = 1 λ 2 π - + log 1 + σ x p w e - w 2 2 d w .

While not directly integrable, we can make progress by Taylor-expanding the logarithm. Firstly we recognise that as the Gaussian distribution is even in w, then all odd terms in the Taylor expansion will vanish. Secondly we note that σxpw=Δxxp will be small if, as above, σ is small, so we can drop terms of fourth order and higher in w. Thus we end up with

(B4) Δ t = - σ 2 2 λ x p 2 2 π - + w 2 e - w 2 2 d w ,

which, from any book of identities, we can find gives

(B5) Δ t = - σ 2 2 λ x p 2 2 π 2 π = - σ 2 2 λ x p 2 .
Code availability

Python code for deriving Eq. (1) and Fortran code for integrating it are available at (Dorrington2023).

Data availability

No data sets were used in this article.


The supplement related to this article is available online at:

Author contributions

JD performed all simulations and drafted the manuscript, under the supervision of TP.

Competing interests

The contact author has declared that neither of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Interdisciplinary perspectives on climate sciences – highlighting past and current scientific achievements”. It is not associated with a conference.


The authors would like to thank Glenn Shutts, Daan Crommelin, and Frank Kwasniok for enlightening discussions. We would like to thank one anonymous reviewer and Tamás Bódai for valuable feedback that has helped to strengthen the manuscript considerably.

Financial support

This research has been supported by the Natural Environment Research Council, National Centre for Earth Observation (grant no. NE/L002612/1).

The article processing charges for this open-access publication were covered by the Karlsruhe Institute of Technology (KIT).

Review statement

This paper was edited by Lesley De Cruz and reviewed by Tamás Bódai and one anonymous referee.


Altmann, E. G. and Endler, A.: Noise-enhanced trapping in chaotic scattering, Phys. Rev. Lett., 105, 244102,, 2010. a, b, c, d

Berner, J., Jung, T., and Palmer, T. N.: Systematic Model Error: The Impact of Increased Horizontal Resolution versus Improved Stochastic and Deterministic Parameterizations, J. Climate, 25, 4946–4962,, 2012. a

Branstator, G.: Circumglobal Teleconnections, the Jet Stream Waveguide, and the North Atlantic Oscillation, Tech. Rep., 14,<1893:CTTJSW>2.0.CO;2, 2002. a

Cehelsky, P. and Tung, K. K.: Theories of multiple equilibria and weather regimes – a critical reexamination. Part II: baroclinic two-layer models, J. Atmos. Sci., 44, 3282–3303,<3282:TOMEAW>2.0.CO;2, 1987. a

Champneys, A. R. and Kirk, V.: The entwined wiggling of homoclinic curves emerging from saddle-node/Hopf instabilities, Physica D, 195, 77–105,, 2004. a

Charney, J. G. and DeVore, J. G.: Multiple Flow Equilibria in the Atmosphere and Blocking, J. Atmos. Sci., 36, 1205–1216,<1205:MFEITA>2.0.CO;2, 1979. a

Charney, J. G. and Straus, D. M.: Form-drag instability, multiple equilibria and propagating planetary waves in baroclinic, orographically forced, planetary wave systems, J. Atmos. Sci., 37, 1157–1176,<1157:FDIMEA>2.0.CO;2, 1980. a

Christensen, H. M., Moroz, I. M., and Palmer, T. N.: Simulating weather regimes: impact of stochastic and perturbed parameter schemes in a simple atmospheric model, Clim. Dynam., 44, 2195–2214,, 2015. a

Crommelin, D. T.: Homoclinic Dynamics: A Scenario for Atmospheric Ultralow-Frequency Variability, J. Atmos. Sci., 59, 1533–1549,<1533:HDASFA>2.0.CO;2, 2002. a

Crommelin, D. T., Opsteegh, J. D., and Verhulst, F.: A Mechanism for Atmospheric Regime Behavior, J. Atmos. Sci., 61, 1406–1419,<1406:amfarb>;2, 2004. a

Cvitanović, P., Søndergaard, N., Palla, G., Vattay, G., and Dettmann, C. P.: Spectrum of stochastic evolution operators: Local matrix representation approach, Phys. Rev. E, 60, 3936,, 1999. a

Cvitanović, P., Artuso, R., Mainieri, R., Tanner, G., and Vattay, G.: Chaos: Classical and Quantum, Niels Bohr Inst., (last access: 3 February 2023), 2016. a

Dawson, A. and Palmer, T. N.: Simulating weather regimes: impact of model resolution and stochastic parameterization, Clim. Dynam., 44, 2177–2193,, 2015. a

De Swart, H. E.: Low-order spectral models of the atmospheric circulation: A survey, Acta Applicandae Mathematicae, 11, 49–96,, 1988. a, b

Dorrington, J.: Software for “On the interaction of stochastic forcing and regime dynamics”, Zenodo [code],, 2023. a

Dorrington, J., Strommen, K., and Fabiano, F.: Quantifying climate model representation of the wintertime Euro-Atlantic circulation using geopotential-jet regimes, Weather Clim. Dynam., 3, 505–533,, 2022. a

Düben, P. D., McNamara, H., and Palmer, T. N.: The use of imprecise processing to improve accuracy in weather & climate prediction, J. Comput. Phys., 271, 2–18,, 2014. a

Faisst, H. and Eckhardt, B.: Lifetimes of noisy repellors, Phys. Rev. E, 68, 026215,, 2003. a

Franaszek, M. and Fronzoni, L.: Influence of noise on crisis-induced intermittency, Phys. Rev. E, 49, 3888,, 1994. a

Hoskins, B. J. and Ambrizzi, T.: Rossby Wave Propagation on a Realistic Longitudinally Varying Flow, J. Atmos. Sci., 50, 1661–1671,<1661:RWPOAR>2.0.CO;2, 1993. a

Itoh, H. and Kimoto, M.: Multiple Attractors and Chaotic Itinerancy in a Quasigeostrophic Model with Realistic Topography: Implications for Weather Regimes and Low-Frequency Variability, J. Atmos. Sci., 53, 2217–2231,<2217:maacii>;2, 1996. a

Itoh, H. and Kimoto, M.: Chaotic itinerancy with preferred transition routes appearing in an atmospheric model, Physica D, 109, 274–292,, 1997. a

Izhikevich, E. M.: Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting, MIT Press, ISBN 978-0-262-09043-8, 2006. a

Kallen, E.: The Nonlinear Effects of Orographic and Momentum Forcing in a Low-Order, Barotropic Model, J. Atmos. Sci., 38, 2150–2163, 1981. a

Kallenberg, O.: Random Measures, Theory and Applications, vol. 77 of Probability Theory and Stochastic Modelling, Springer International Publishing, Cham,, 2017. a

Kimoto, M. and Ghil, M.: Multiple Flow Regimes in the Northern Hemisphere Winter. Part I: Methodology and Hemispheric Regimes, J. Atmos. Sci., 50, 2625–2644,<2625:mfritn>;2, 1993. a

Kondrashov, D., Ide, K., and Ghil, M.: Weather Regimes and Preferred Transition Paths in a Three-Level Quasigeostrophic Model, J. Atmos. Sci., 61, 568–587,<0568:WRAPTP>2.0.CO;2, 2004. a

Kwasniok, F.: Enhanced regime predictability in atmospheric low-order models due to stochastic forcing, Philos. T. Roy. Soc. A, 372, 20130286–20130286,, 2014. a, b

Lai, Y. C. and Tél, T.: Introduction to Transient Chaos, Appl. Math. Sci.-Switzerland, 173, 3–35,, 2011a. a

Lai, Y. C. and Tél, T.: Noise and Transient Chaos, Appl. Math. Sci.-Switzerland, 173, 107–143,, 2011b. a, b, c, d

Lorenz, E. N.: Deterministic nonperiodic flow, Universality in Chaos, 2nd edn., 20, 367–378,, 1963. a, b

Lucarini, V. and Gritsun, A.: A new mathematical framework for atmospheric blocking events, Clim. Dynam., 54, 575–598,, 2019. a

Maiocchi, C. C., Lucarini, V., and Gritsun, A.: Decomposing the dynamics of the Lorenz 1963 model using unstable periodic orbits: Averages, transitions, and quasi-invariant sets, Chaos: An Interdisciplinary Journal of Nonlinear Science, 32, 033129,, 2022. a

Meurer, A., Smith, C. P., Paprocki, M., Čertík, O., Kirpichev, S. B., Rocklin, M., Kumar, A. T., Ivanov, S., Moore, J. K., Singh, S., Rathnayake, T., Vig, S., Granger, B. E., Muller, R. P., Bonazzi, F., Gupta, H., Vats, S., Johansson, F., Pedregosa, F., Curry, M. J., Terrel, A. R., Roučka, Š., Saboo, A., Fernando, I., Kulal, S., Cimrman, R., and Scopatz, A.: SymPy: Symbolic computing in python, PeerJ Comput. Sci., 2017, e103,, 2017.  a

Ott, E.: Chaotic transitions, in: Chaos in Dynamical Systems, edited by: Ott, E., Cambridge University Press, Cambridge, 2nd edn., 283–294,, 2002. a

Palmer, T. N. and Weisheimer, A.: Diagnosing the causes of bias in climate models – why is it so hard?, Geophys. Astrophys. Fluid Dynam., 105, 351–365,, 2011. a

Pickl, M., Lang, S. T., Leutbecher, M., and Grams, C. M.: The effect of stochastically perturbed parametrisation tendencies (SPPT) on rapidly ascending air streams, Q. J. Roy. Meteor. Soc., 148, 1242–1261,, 2022. a

Pusuluri, K. and Shilnikov, A.: Homoclinic chaos and its organization in a nonlinear optics model, Phys. Rev. E, 98, 040202,, 2018. a

Reimann, P.: Noisy one-dimensional maps near a crisis. I. Weak Gaussian white and colored noise, J. Statist. Phys., 82, 1467–1501,, 1996. a, b, c

Reinhold, B. B. and Pierrehumbert, R. T.: Dynamics of Weather Regimes: Quasi-Stationary Waves and Blocking, Mon. Weather Rev., 110, 1105–1145,<1105:dowrqs>;2, 1982. a

Romeiras, F. J., Grebogi, C., and Ott, E.: Multifractal properties of snapshot attractors of random maps, Phys. Rev. A, 41, 784–799,, 1990. a

Rossby, C. G.: Planetary flow patterns in the atmosphere, Q. J. Roy. Meteor. Soc., 66, 68–87, 1940. a

Selten, F. M. and Branstator, G.: Preferred Regime Transition Routes and Evidence for an Unstable Periodic Orbit in a Baroclinic Model, J. Atmos. Sci., 61, 2267–2282,<2267:prtrae>;2, 2004. a

Shen, B. W., Pielke, R. A., Zeng, X., Baik, J. J., Faghih-Naini, S., Cui, J., and Atlas, R.: Is weather chaotic? Coexistence of chaos and order within a generalized lorenz model, B. Am. Meteorol. Soc., 102, E148–E158,, 2021. a

Shilnikov, A., Nicolis, G., and Nicolis, C.: Bifurcation and predictability analysis of a low-order atmopsheric circulation model, Int. J. Bifurcation and Chaos, 05, 1701–1711,, 1995. a, b

Strommen, K., Chantry, M., Dorrington, J., and Otter, N.: A topological perspective on weather regimes, Clim. Dynam.,, 2022. a

Wackerbauer, R. and Kobayashi, S.: Noise can delay and advance the collapse of spatiotemporal chaos, Phys. Rev. E, 75, 066209,, 2007. a

Yang, S., Reinhold, B., and Källé, E.: Multiple Weather Regimes and Baroclinically Forced Spherical Resonance, J. Atmos. Sci., 54, 1397–1409,<1397:MWRABF>2.0.CO;2, 1997. a


The “pullback attractor” of Romeiras et al. (1990) gives formal meaning to the concept of an attractor for stochastically forced chaotic systems.

Short summary
Atmospheric models often include random forcings, which aim to replicate the impact of processes too small to be resolved. Recent results in simple atmospheric models suggest that this random forcing can actually stabilise certain slow-varying aspects of the system, which could provide a path for resolving known errors in our models. We use randomly forced simulations of a toy chaotic system and theoretical arguments to explain why this strange effect occurs – at least in simple models.