the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
On the interaction of stochastic forcing and regime dynamics
Joshua Dorrington
Tim Palmer
Stochastic forcing can, sometimes, stabilise atmospheric regime dynamics, increasing their persistence. This counterintuitive 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 sixmode 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 lowdimensional 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.
 Article
(3782 KB) 
Supplement
(1164 KB)  BibTeX
 EndNote
Rossby (1940) first suggested the atmosphere may possess multiple largescale quasiequilibrium 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 crisisinduced intermittency (Ott, 2002). The first theoretical study of multiple atmospheric equilibria, carried out by Charney and DeVore (1979) (hereafter CdV79), showed that the flow in a loworder midlatitude barotropic channel could equilibrate to either a zonally symmetric flow or a stationary blocking pattern. They suggested realworld 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 largescale weather patterns.
Similar results were obtained in other simple models (Charney and Straus, 1980; Kallen, 1981; De Swart, 1988), 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 Pierrehumbert, 1982; Cehelsky and Tung, 1987). Therefore later work moved away from analysis of fixed points, treating regimes as chaotic flows originating from attractor merging (Itoh and Kimoto, 1996, 1997) and exploring the role of homoclinic orbits in preferred transitions between regimes (Kimoto and Ghil, 1993; Kondrashov et al., 2004; Selten and Branstator, 2004).
Despite a deepening mathematical understanding of regimes, the Earth system appears to be highdimensional, 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. Earthsystem models often contain explicit stochastic parameterisation schemes and in addition a range of smallscale processes that can be thought of as randomly forcing the largescale 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 (Lorenz, 1963) 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 Palmer, 2015), as well as the representation of other slowvarying climate modes (Palmer and Weisheimer, 2011; 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él, 2011a). It has been shown for such systems how stochastic forcing can delay the escape of trajectories in the case of 1D (Reimann, 1996; Lai and Tél, 2011b) and 2D (Altmann and Endler, 2010) 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 sixdimensional, continuoustime, nonhyperbolic 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.
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:
The six variables {x_{i}} 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 nonlinearities. 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:
The system is integrated using the Euler–Maruyama difference scheme with a time step of $\mathrm{d}t=\mathrm{2}\times {\mathrm{10}}^{\mathrm{4}}$ model time units (MTU) and is sampled every 1 MTU. We compute integrations of length 200 000 MTU for a range of σ values.
To visualise the model phase space, we project the sixdimensional 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 lowdimensional 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 scrolllike 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 nearzero phase speed, corresponding to persistent “blocked” flow patterns, divided by highly chaotic, fastevolving 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 saddlenode bifurcation (a transition from one to multiple fixedpoint 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 steadystate 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 saddlenode bifurcations merge: one fixed point bifurcates into two fixed points and a limit cycle. This foldHopf point has codimension 2, meaning two model parameters must be tuned to obtain it. Highcodimension bifurcations seem to be of reduced relevance to the real world, representing exceptional states in parameter space. However, the foldHopf 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 Kirk, 2004). 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 (Crommelin, 2002) including the Lorenz 1984 model (Shilnikov et al., 1995) and in the dynamics of nonlinear circuits (Pusuluri and Shilnikov, 2018) and neurons (Izhikevich, 2006).
This type of regime behaviour is of a fundamentally different type to that of the archetypal Lorenz 1963 system. In Lorenz 1963, crisisinduced 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 reenter the slowevolving regime. For a recent discussion of the coexistence 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 sixdimensional state vector, rather than restricting it to the x_{1}–x_{4} 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 quasipredictable 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.
For nonzero 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 attractor^{1} 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 x_{1} mode, and the phase velocity ${v}_{t}=\frac{\mathrm{1}}{\mathrm{2}}\sqrt{{\sum}_{i=\mathrm{1}}^{\mathrm{6}}({x}_{i,t+\mathrm{1}}{x}_{i,t\mathrm{1}}{)}^{\mathrm{2}}}$, for different σ values, coloured by regime assignment. The characteristics of the regimes are broadly preserved even up to $\mathit{\sigma}=\mathrm{1.6}\times {\mathrm{10}}^{\mathrm{2}}$: the persistent blocking state has low x_{1} values and very low phase velocities, and the zonal state has high x_{1} values and relatively low phase velocities, while the transitional regime has high x_{1} values and high phase velocity. However, for the larger value of $\mathit{\sigma}=\mathrm{3.0}\times {\mathrm{10}}^{\mathrm{2}}$, it is apparent that there is no meaningful regime structure remaining: the phase velocity never drops below 0.2, indicating a loss of quasistationarity. At this point the deterministic dynamics play second fiddle to the stochastic term, and the assignment of regimes becomes meaningless.
Figure 5a shows Fourier transforms of the total mode amplitude $\sqrt{{\sum}_{i}{x}_{i}^{\mathrm{2}}}$ for three representative amplitudes of stochastic forcing. The deterministic σ=0 spectrum possesses a clear peak consistent with the quasiperiodic 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 lowfrequency region decreases, and the spectrum flattens. Figure 5b shows this nonmonotonicity clearly, with a peak in lowfrequency variability close to σ=0.01.
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 lowfrequency 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.
How can the ability of random forcing to stabilise regime dynamics be understood? What causes the nonmonotonicity between stochastic amplitude and persistence, and why is the response here asymmetric, primarily in the blocking regime? As mentioned, nonmonotonic changes of transient lifetime due to stochastic forcing are understood for lowdimensional 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 (Reimann, 1996; Lai and Tél, 2011b). 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 Endler, 2010; Faisst and Eckhardt, 2003).
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 sixdimensional 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 dynamicsdominated 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 fattailed, and so the mean diverges (Kallenberg, 2017). 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 sixdimensional 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 noisedominated neighbourhood of the fixed point. There is also an important interaction between stochasticity and the lowdimensional, 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 slowevolving dynamics, the system is at a point x=x_{i} and will stay in the quasistationary regime until it reaches x=x_{c}. 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:x\to t$, $g:t\to x$. If we let x=f(t) be any convex, monotonic function defined over the interval [x_{0},x_{c}], then $t=g\left(x\right)={f}^{\mathrm{1}}\left(x\right)$, 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 x_{p}, drawn from an even probability distribution P(Δx) with fastdecaying tails and support over [x_{0}−x_{p}, x_{c}−x_{p}]. Informally, this amounts to choosing a point x_{p} 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 x_{c}−x_{0}, 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:
However this is not in general the case for 〈Δt〉:
where G_{p}, as a translation of g, is also concave. Using the definition of the expectation we find
where the limits are justified by our requirement for vanishing tails for P and our choice of x_{p}.
Because P(Δx) is even, then if we Taylorexpand G_{p}, we will lose all odd powers on integration and be left with
And having taken Δx to be small, we can drop highorder terms:
Because G_{p} is concave, this derivative is negative, and we obtain
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)=x_{0}e^{λt}, and $P\left(\mathrm{\Delta}x\right)=\mathcal{N}(\mathrm{0},\mathit{\sigma})$, which is interesting for our system, we obtain analytically $\langle \mathrm{\Delta}t\rangle =\frac{{\mathit{\sigma}}^{\mathrm{2}}}{\mathrm{2}\mathit{\lambda}{x}_{\mathrm{p}}^{\mathrm{2}}}$ (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 noninvertible 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 x_{p} 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 lowdimensional, 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 noisedominated regime where the stochastic terms are much larger than the deterministic tendency.
In this paper we have investigated the interaction between stochasticity and strongly nonlinear dynamics in a sixequation 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 quasistationary 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 noisedominated, 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 Endler, 2010), 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 nonattracting 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 continuoustime fully chaotic conceptual models such as the archetypal Lorenz 1963 system (Lorenz, 1963), 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 discretetime 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 Gritsun, 2019) and regime systems (Maiocchi et al., 2022). Characterising the unstable periodic orbits associated with each regime could even allow the nonmonotonic 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 EuroAtlantic region (Dorrington et al., 2022). A natural next step is to ask how stochasticity impacts regimes in quasigeostrophic 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 nearcrisis 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 spatiotemporal 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.
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 massconservation equation.
The momentum equation in a rotating frame is given by
while the massconservation equation is
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 $(\mathbf{a}\cdot \mathbf{b})\mathbf{c}=\mathbf{a}\times (\mathbf{b}\times \mathbf{c})+(\mathbf{a}\cdot \mathbf{c})\mathbf{b}$ to Eq. (A1) gives us
and by introducing the absolute vorticity, $\mathit{\xi}=\mathrm{\nabla}\times \mathbf{v}$, and taking the curl of both sides, we get the shallow water vorticity equation:
Again, we apply the triple product identity and the chain rule and find
where $\mathrm{\nabla}\cdot \mathit{\xi}=\mathrm{\nabla}\cdot (\mathrm{\nabla}\times \mathbf{v})=\mathrm{0}$ by definition and $(\mathit{\xi}\cdot \mathrm{\nabla})\mathbf{v}=\mathrm{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:=f_{z} for convenience. The divergent velocity term can be rewritten with the aid of Eq. (A2) by rearranging
and substituting it into Eq. (A5). This leaves us with
which can be put into a clearer form using the chain rule:
That is, the potential vorticity $q=\frac{\mathit{\xi}+f}{H}$ 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=Dh(x,y)$. We assume shallow orography, so we may take the limit where $\frac{h}{D}\ll \mathrm{1}$ and Taylorexpand H, truncating at first order so that $\frac{\mathrm{1}}{H}\approx \frac{\mathrm{1}}{D}+\frac{h}{{D}^{\mathrm{2}}}$:
In the midlatitude 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 quasigeostrophic 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:
We define $\mathit{\gamma}=\frac{{f}_{\mathrm{0}}}{D}$ and introduce the streamfunction Ψ, given by $u=\frac{\mathrm{\Psi}}{y}$, $v=\frac{\mathrm{\Psi}}{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
or more concisely
where the Jacobian operator is given by
Taking our forcing to be a linear relaxation towards a background state, we now have our explicit equation for the evolution of Ψ:
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 nopenetration and noslip conditions in the meridional direction. Of course the real atmosphere has no hard meridional boundaries; however the impacts of the subtropical and midlatitude jets have been shown to produce a waveguiding effect that meridionally confines lowfrequency variability in a similar way, both in loworder models (Hoskins and Ambrizzi, 1993; Yang et al., 1997) and in observations (Branstator, 2002), which makes this a reasonable first approximation.
We consider a channel of dimensionless zonal length $x=[\mathrm{0},\mathrm{2}\mathit{\pi}]$ and meridional length $y=[\mathrm{0},b\mathit{\pi}]$, 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 nonlinearly 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:
The projection coefficient of the jth and kth modes on the ith mode, c_{ijk}, is computable simply through a direct pointwise product of the ith mode and the Jacobian, integrated over the domain and appropriately normalised:
Under these constraints, we can now rewrite Eq. (A14) projected onto the ith basis function as
where the double sum indexes over all basis functions. We make an explicit choice of basis functions now, taking a Fourier basis defined as
where the meridional mode n>0 and the zonal mode m can take any noninteger 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 $\leftm\right\le \mathrm{1}$ and a meridional threshold of n≤2, restricting us to just six basis functions, as in CdV79, C04, and K14.
Although up to this point our approximations have been motivated well physically, we cannot rigorously justify introducing such a drastic truncation without resorting to selfinterest: such a system is much easier to study and can be understood in more detail than a highdimensional 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 nonlinear 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 poletoequator 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)=\mathit{\gamma}\sqrt{\mathrm{2}}\mathrm{cos}\left(x\right)\mathrm{sin}\left(\frac{y}{b}\right)$ (this is shown in Fig. A2).
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 https://github.com/joshdorrington/Thesis_materials/blob/main/10d_mode_derivation.ipynb, 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 sixequation system:
where colours encode the linear relaxation (green), Coriolis forces (blue), orographic impacts (orange), and nonlinear advection (red). We see that there are far fewer than the 36 promised nonlinear terms per equation: most quadratic interactions integrated over the domain vanish due to symmetry considerations. $\mathit{\left\{}{\mathit{\alpha}}_{n}\mathit{\right\}},\mathit{\left\{}{\mathit{\beta}}_{n}\mathit{\right\}},\mathit{\left\{}{\mathit{\delta}}_{n}\mathit{\right\}},\mathit{\left\{}{\mathit{\gamma}}_{n}\mathit{\right\}},\mathit{\left\{}\stackrel{\mathrm{\u0303}}{{\mathit{\gamma}}_{n}}\mathit{\right\}}$, 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:

${x}_{i}^{*}$ is the forcing vector, representing thermal relaxation towards a background state. We use a zonally symmetric forcing profile with x_{1}=0.95 and ${x}_{\mathrm{4}}=\mathrm{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 halfwidth, set to b=0.5, which gives a channel of 6300 km × 1600 km.
In Sect. 3 we found that
If we consider x to be undergoing exponential growth and subject it to a Gaussian perturbation so that $x=f\left(t\right)={x}_{\mathrm{0}}{e}^{\mathit{\lambda}t}$, and $P\left(\mathrm{\Delta}x\right)=\mathcal{N}(\mathrm{0},\mathit{\sigma})$, then we can obtain an analytical form for Δt. Remembering that ${G}_{p}\left(\mathrm{\Delta}x\right)=g({x}_{\mathrm{p}}+\mathrm{\Delta}x)g\left({x}_{\mathrm{p}}\right)$, and seeing that $g\left(x\right)=\frac{\mathrm{1}}{\mathit{\lambda}}\mathrm{log}\left(\frac{x}{{x}_{\mathrm{0}}}\right)$, we have
We transform our variables to $w=\frac{\mathrm{\Delta}x}{\mathit{\sigma}}=\frac{x{x}_{\mathrm{p}}}{\mathit{\sigma}}$, which implies $\frac{x}{{x}_{\mathrm{p}}}=\frac{{x}_{\mathrm{p}}+\mathit{\sigma}w}{{x}_{\mathrm{p}}}=\mathrm{1}+\frac{\mathit{\sigma}}{{x}_{\mathrm{p}}}w$ and $dx=\mathit{\sigma}\cdot \mathrm{d}w$. Our limits of integration become ${w}_{\mathrm{0}}=\frac{{x}_{\mathrm{0}}{x}_{\mathrm{p}}}{\mathit{\sigma}}<\mathrm{0}$ and ${w}_{\mathrm{c}}=\frac{{x}_{\mathrm{c}}{x}_{\mathrm{p}}}{\mathit{\sigma}}>\mathrm{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 ${w}_{\mathrm{0}}\approx \mathrm{\infty}$, ${w}_{\mathrm{c}}\approx +\mathrm{\infty}$. This gives us
While not directly integrable, we can make progress by Taylorexpanding 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 $\frac{\mathit{\sigma}}{{x}_{\mathrm{p}}}w=\frac{\mathrm{\Delta}x}{{x}_{\mathrm{p}}}$ 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
which, from any book of identities, we can find gives
Python code for deriving Eq. (1) and Fortran code for integrating it are available at https://doi.org/10.5281/zenodo.7602855 (Dorrington, 2023).
No data sets were used in this article.
The supplement related to this article is available online at: https://doi.org/10.5194/npg30492023supplement.
JD performed all simulations and drafted the manuscript, under the supervision of TP.
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.
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.
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 openaccess publication were covered by the Karlsruhe Institute of Technology (KIT).
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.: Noiseenhanced trapping in chaotic scattering, Phys. Rev. Lett., 105, 244102, https://doi.org/10.1103/PhysRevLett.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, https://doi.org/10.1175/JCLID1100297.1, 2012. a
Branstator, G.: Circumglobal Teleconnections, the Jet Stream Waveguide, and the North Atlantic Oscillation, Tech. Rep., 14, https://doi.org/10.1175/15200442(2002)015<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 twolayer models, J. Atmos. Sci., 44, 3282–3303, https://doi.org/10.1175/15200469(1987)044<3282:TOMEAW>2.0.CO;2, 1987. a
Champneys, A. R. and Kirk, V.: The entwined wiggling of homoclinic curves emerging from saddlenode/Hopf instabilities, Physica D, 195, 77–105, https://doi.org/10.1016/j.physd.2004.03.004, 2004. a
Charney, J. G. and DeVore, J. G.: Multiple Flow Equilibria in the Atmosphere and Blocking, J. Atmos. Sci., 36, 1205–1216, https://doi.org/10.1175/15200469(1979)036<1205:MFEITA>2.0.CO;2, 1979. a
Charney, J. G. and Straus, D. M.: Formdrag instability, multiple equilibria and propagating planetary waves in baroclinic, orographically forced, planetary wave systems, J. Atmos. Sci., 37, 1157–1176, https://doi.org/10.1175/15200469(1980)037<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, https://doi.org/10.1007/s0038201422399, 2015. a
Crommelin, D. T.: Homoclinic Dynamics: A Scenario for Atmospheric UltralowFrequency Variability, J. Atmos. Sci., 59, 1533–1549, https://doi.org/10.1175/15200469(2002)059<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, https://doi.org/10.1175/15200469(2004)061<1406:amfarb>2.0.co;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, https://doi.org/10.1103/PhysRevE.60.3936, 1999. a
Cvitanović, P., Artuso, R., Mainieri, R., Tanner, G., and Vattay, G.: Chaos: Classical and Quantum, Niels Bohr Inst., http://chaosbook.org/ (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, https://doi.org/10.1007/s003820142238x, 2015. a
De Swart, H. E.: Loworder spectral models of the atmospheric circulation: A survey, Acta Applicandae Mathematicae, 11, 49–96, https://doi.org/10.1007/BF00047114, 1988. a, b
Dorrington, J.: Software for “On the interaction of stochastic forcing and regime dynamics”, Zenodo [code], https://doi.org/10.5281/zenodo.7602855, 2023. a
Dorrington, J., Strommen, K., and Fabiano, F.: Quantifying climate model representation of the wintertime EuroAtlantic circulation using geopotentialjet regimes, Weather Clim. Dynam., 3, 505–533, https://doi.org/10.5194/wcd35052022, 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, https://doi.org/10.1016/j.jcp.2013.10.042, 2014. a
Faisst, H. and Eckhardt, B.: Lifetimes of noisy repellors, Phys. Rev. E, 68, 026215, https://doi.org/10.1103/PhysRevE.68.026215, 2003. a
Franaszek, M. and Fronzoni, L.: Influence of noise on crisisinduced intermittency, Phys. Rev. E, 49, 3888, https://doi.org/10.1103/PhysRevE.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, https://doi.org/10.1175/15200469(1993)050<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 LowFrequency Variability, J. Atmos. Sci., 53, 2217–2231, https://doi.org/10.1175/15200469(1996)053<2217:maacii>2.0.co;2, 1996. a
Itoh, H. and Kimoto, M.: Chaotic itinerancy with preferred transition routes appearing in an atmospheric model, Physica D, 109, 274–292, https://doi.org/10.1016/S01672789(97)00064X, 1997. a
Izhikevich, E. M.: Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting, MIT Press, ISBN 9780262090438, 2006. a
Kallen, E.: The Nonlinear Effects of Orographic and Momentum Forcing in a LowOrder, 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, https://doi.org/10.1007/9783319415987, 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, https://doi.org/10.1175/15200469(1993)050<2625:mfritn>2.0.co;2, 1993. a
Kondrashov, D., Ide, K., and Ghil, M.: Weather Regimes and Preferred Transition Paths in a ThreeLevel Quasigeostrophic Model, J. Atmos. Sci., 61, 568–587, https://doi.org/10.1175/15200469(2004)061<0568:WRAPTP>2.0.CO;2, 2004. a
Kwasniok, F.: Enhanced regime predictability in atmospheric loworder models due to stochastic forcing, Philos. T. Roy. Soc. A, 372, 20130286–20130286, https://doi.org/10.1098/rsta.2013.0286, 2014. a, b
Lai, Y. C. and Tél, T.: Introduction to Transient Chaos, Appl. Math. Sci.Switzerland, 173, 3–35, https://doi.org/10.1007/9781441969873_1, 2011a. a
Lai, Y. C. and Tél, T.: Noise and Transient Chaos, Appl. Math. Sci.Switzerland, 173, 107–143, https://doi.org/10.1007/9781441969873_4, 2011b. a, b, c, d
Lorenz, E. N.: Deterministic nonperiodic flow, Universality in Chaos, 2nd edn., 20, 367–378, https://doi.org/10.1201/9780203734636, 1963. a, b
Lucarini, V. and Gritsun, A.: A new mathematical framework for atmospheric blocking events, Clim. Dynam., 54, 575–598, https://doi.org/10.1007/S00382019050182, 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 quasiinvariant sets, Chaos: An Interdisciplinary Journal of Nonlinear Science, 32, 033129, https://doi.org/10.1063/5.0067673, 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, https://doi.org/10.7717/peerjcs.103, 2017. a
Ott, E.: Chaotic transitions, in: Chaos in Dynamical Systems, edited by: Ott, E., Cambridge University Press, Cambridge, 2nd edn., 283–294, https://doi.org/10.1017/CBO9780511803260.010, 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, https://doi.org/10.1080/03091929.2010.547194, 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, https://doi.org/10.1002/QJ.4257, 2022. a
Pusuluri, K. and Shilnikov, A.: Homoclinic chaos and its organization in a nonlinear optics model, Phys. Rev. E, 98, 040202, https://doi.org/10.1103/PhysRevE.98.040202, 2018. a
Reimann, P.: Noisy onedimensional maps near a crisis. I. Weak Gaussian white and colored noise, J. Statist. Phys., 82, 1467–1501, https://doi.org/10.1007/BF02183392, 1996. a, b, c
Reinhold, B. B. and Pierrehumbert, R. T.: Dynamics of Weather Regimes: QuasiStationary Waves and Blocking, Mon. Weather Rev., 110, 1105–1145, https://doi.org/10.1175/15200493(1982)110<1105:dowrqs>2.0.co;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, https://doi.org/10.1103/PhysRevA.41.784, 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, https://doi.org/10.1175/15200469(2004)061<2267:prtrae>2.0.co;2, 2004. a
Shen, B. W., Pielke, R. A., Zeng, X., Baik, J. J., FaghihNaini, 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, https://doi.org/10.1175/BAMSD190165.1, 2021. a
Shilnikov, A., Nicolis, G., and Nicolis, C.: Bifurcation and predictability analysis of a loworder atmopsheric circulation model, Int. J. Bifurcation and Chaos, 05, 1701–1711, https://doi.org/10.1142/S0218127495001253, 1995. a, b
Strommen, K., Chantry, M., Dorrington, J., and Otter, N.: A topological perspective on weather regimes, Clim. Dynam., https://doi.org/10.1007/s0038202206395x, 2022. a
Wackerbauer, R. and Kobayashi, S.: Noise can delay and advance the collapse of spatiotemporal chaos, Phys. Rev. E, 75, 066209, https://doi.org/10.1103/PhysRevE.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, https://doi.org/10.1175/15200469(1997)054<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.
 Abstract
 Introduction
 Model formulation
 Stochastic persistence
 Conclusions
 Appendix A: Derivation of the Charney–deVore 1979 model
 Appendix B: Impact of a Gaussian perturbation on exponential growth
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Special issue statement
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
toychaotic system and theoretical arguments to explain why this strange effect occurs – at least in simple models.
 Abstract
 Introduction
 Model formulation
 Stochastic persistence
 Conclusions
 Appendix A: Derivation of the Charney–deVore 1979 model
 Appendix B: Impact of a Gaussian perturbation on exponential growth
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Special issue statement
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement