Many physical systems are forced by external inputs, which can sometimes take the form of chaotic variation. A particular example is found in applications related to weather and climate, where chaotic variation is prevalent across various timescales. If the system in question has multiple attracting solutions for a given range of forcing, rate-induced tipping can be triggered by the chaotic forcing, with the difference in timescales between the forcing and the system acting as a `rate' parameter. In this paper, we explore the interplay between these two timescales in a low-order model of ice age dynamics. The model exhibits bistability between two equilibria in one region of the parameter space and between an equilibrium and a periodic orbit in another region. When chaotic variation of the parameters is allowed within these bistable regions, the solutions of the forced system undergo rate-induced tipping from one attractor to another. Simulations of the forced system show that the timescale of the chaotic forcing induces a resonance-like behaviour, with an optimal timescale at which the likelihood of rate-induced tipping is at its maximum. We combine basin instability theory, finite-time Lyapunov exponents, and linear resonance analysis under periodic forcing to explain this resonance effect.

We consider a model thermokinetic reaction that is both autocatalytic and temperature sensitive. The mechanism consists of two reactions. The first is a reversible reaction with a temperature-dependent rate, and the second involves a cubically autocatalytic process. We demonstrate the presence of up to five steady states in the system, dependent on parameters, and we illustrate the fact that the system can exhibit tipping-point behaviour, in which small changes to initial conditions in an experiment could lead to radically different long-term outcomes. We prove that oscillatory behaviour is not possible over large regions in the parameter space; outside these regions, however, oscillations have been found, and we present a simple and robust method for computing them. Small-amplitude oscillations can arise through Hopf bifurcations at an equilibrium point. In addition, large-amplitude oscillations can bifurcate directly from global structures that are associated with equilibria behaving as saddles. This is illustrated with some numerical solutions of these highly nonlinear equations.

We consider a coupling of the Stommel box model and the Lorenz model, with the goal of investigating the so-called crises that are known to occur given sufficient forcing. In this context, a crisis is characterized as the destruction of a chaotic attractor under a critical forcing strength. We document the variety of chaotic attractors and crises possible in our model, focusing on the parameter region where the Lorenz model is always chaotic and where bistability exists in the Stommel box model. The chaotic saddle collisions that occur in a boundary crisis are visualized, with the chaotic saddle computed using the Saddle-Straddle Algorithm. We identify a novel sub-type of boundary crisis, namely a vanishing basin crisis. For forcing strength beyond the crisis, we demonstrate the possibility of a merging between the persisting chaotic attractor and either a chaotic transient or a ghost attractor depending on the type of boundary crisis. An investigation of the finite-time Lyapunov exponents around crisis levels of forcing reveals a convergence between two near-neutral exponents, particularly at points of a trajectory most sensitive to divergence. This points to loss of hyperbolicity associated with crisis occurrence. Finally, we generalize our findings by coupling the Stommel box model to other strange attractors and thereby show that the behaviors are quite generic and robust.

It is known that the widely studied Bray–Liebhafsky reaction typically exhibits complex chemical behaviour. Numerous mathematical systems have been proposed to describe the iodine oscillations that occur during this process. Recently, a four-variable model of the Bray–Liebhafsky reaction has been proposed and analytical and numerical investigations suggested that chaotic solutions may exist. We revisit this four-variable model here and perform what appears to be the first detailed work on this system. We suggest that this model is perhaps not chaotic after all. Informed by these fresh insights, we propose a reduced two-variable model based upon the four-variable system. This model is created with the twin goals of enabling simpler mathematical analysis while retaining the underlying chemical mechanisms. We are able to demonstrate that our reduced problem performs very well when compared with the full model for realistic parameter values. In particular, key regions of parameter space are identified within which temporal oscillations can occur. Moreover, these persistent oscillations are consistent with the available qualitative experimental observations.

Singular vectors (SVs) have long been employed in the initialization of ensemble numerical weather prediction (NWP) in order to capture the structural organization and growth rates of those perturbations or “errors” associated with initial condition errors and instability processes of the large scale flow. Due to their (super) exponential growth rates and spatial scales, initial SVs are typically combined empirically with evolved SVs in order to generate forecast perturbations whose structures and growth rates are tuned for specified lead-times. Here we present a systematic approach to generating finite time or "mixed" SVs (MSVs) based on a method for the calculation of covariant Lyapunov vectors (CLVs) and appropriate choices of the matrix cocycle. We first derive a data-driven reduced order model to characterize persistent geopotential height anomalies over Europe and Western Asia (Eurasia) over the period 1979-present from the NCEPv1 reanalysis. We then characterize and compare the MSVs and SVs of each persistent state over Eurasia for particular lead-times from a day to over a week. Finally, we compare the spatio-temporal properties of SVs and MSVs in an examination of the dynamics of the 2010 Russian heatwave. We show that MSVs provide a systematic approach to generate initial forecast perturbations projected onto relevant expanding directions in phase space for typical NWP forecast lead-times.

The CSIRO Climate retrospective Analysis and Forecast Ensemble system: version 1 (CAFE60v1) provides a large (96 member) ensemble retrospective analysis of the global climate system from 1960 to present with sufficiently many realizations and at spatio-temporal resolutions suitable to enable probabilistic climate studies. Using a variant of the ensemble Kalman filter, 96 climate state estimates are generated over the most recent six decades. These state estimates are constrained by monthly mean ocean, atmosphere and sea ice observations such that their trajectories track the observed state while enabling estimation of the uncertainties in the approximations to the retrospective mean climate over recent decades. For the atmosphere, we evaluate CAFE60v1 in comparison to empirical indices of the major climate teleconnections and blocking with various reanalysis products. Estimates of the large scale ocean structure, transports and biogeochemistry are compared to those derived from gridded observational products and climate model projections (CMIP). Sea ice (extent, concentration and variability) and land surface (precipitation and surface air temperatures) are also compared to a variety of model and observational products. Our results show that CAFE60v1 is a useful, comprehensive and unique data resource for studying internal climate variability and predictability, including the recent climate response to anthropogenic forcing on multi-year to decadal time scales.

A new technique to derive delay models from systems of partial differential equations, based on the Mori–Zwanzig (MZ) formalism, is used to derive a delay-difference equation model for the Atlantic Multidecadal Oscillation (AMO). The MZ formalism gives a rewriting of the original system of equations, which contains a memory term. This memory term can be related to a delay term in a resulting delay equation. Here, the technique is applied to an idealized, but spatially extended, model of the AMO. The resulting delay-difference model is of a different type than the delay differential model which has been used to describe the El Niño Southern Oscillation. In addition to this model, which can also be obtained by integration along characteristics, error terms for a smoothing approximation of the model have been derived from the MZ formalism. Our new method of deriving delay models from spatially extended models has a large potential to use delay models to study a range of climate variability phenomena.

The codes that have been used to produce the Figures in the paper.

The dynamics of the North Atlantic Oscillation (NAO) are analyzed through a data-driven model obtained from atmospheric reanalysis data. We apply a regularized vector autoregressive clustering technique to identify recurrent and persistent states of atmospheric circulation patterns in the North Atlantic sector (110W-0E, 20N-90N). In order to analyze the dynamics associated with the resulting cluster-based models, we define a time-dependent linear delayed map with a switching sequence set a priori by the cluster affiliations at each time step. Using a method for computing the covariant Lyapunov vectors (CLVs) over various time windows, we produce sets of mixed singular vectors (for short windows) and approximate the asymptotic CLVs (for longer windows). The growth rates and alignment of the resulting time-dependent vectors are then analyzed. We find that the window chosen to compute the vectors acts as a filter on the dynamics. For short windows, the alignment and changes in growth rates are indicative of individual transitions between persistent states. For long windows, we observe an emergent annual signal manifest in the alignment of the CLVs characteristic of the observed seasonality in the NAO index. Analysis of the average finite-time dimension reveals the NAO- as the most unstable state relative to the NAO+, with persistent AR states largely stable. Our results agree with other recent theoretical and empirical studies that have shown blocking events to have less predictability than periods of enhanced zonal flow.

The basis and challenge of strongly coupled data assimilation (CDA) is the accurate representation of cross-domain covariances between various coupled subsystems with disparate spatio-temporal scales, where often one or more subsystems are unobserved. In this study, we explore strong CDA using ensemble Kalman filtering methods applied to a conceptual multiscale chaotic model consisting of three coupled Lorenz attractors. We introduce the use of the local attractor dimension (i.e. the Kaplan–Yorke dimension, dimKY) to prescribe the rank of the background covariance matrix which we construct using a variable number of weighted covariant Lyapunov vectors (CLVs). Specifically, we consider the ability to track the nonlinear trajectory of each of the subsystems with different variants of sparse observations, relying only on the cross-domain covariance to determine an accurate analysis for tracking the trajectory of the unobserved subdomain. We find that spanning the global unstable and neutral subspaces is not sufficient at times where the nonlinear dynamics and intermittent linear error growth along a stable direction combine. At such times a subset of the local stable subspace is also needed to be represented in the ensemble. In this regard the local dimKY provides an accurate estimate of the required rank. Additionally, we show that spanning the full space does not improve performance significantly relative to spanning only the subspace determined by the local dimension. Where weak coupling between subsystems leads to covariance collapse in one or more of the unobserved subsystems, we apply a novel modified Kalman gain where the background covariances are scaled by their Frobenius norm. This modified gain increases the magnitude of the innovations and the effective dimension of the unobserved domains relative to the strength of the coupling and timescale separation. We conclude with a discussion on the implications for higher-dimensional systems.

Animation of Figure 7; Animation of Figure 12 red trajectory; Animation of Figure 12 black trajectory; Animation of Figure 14 red trajectory; Animation of Figure 14 black trajectory

The Atlantic meridional overturning circulation (AMOC) transports substantial amounts of heat into the North Atlantic sector, and hence is of very high importance in regional climate projections. The AMOC has been observed to show multi-stability across a range of models of different complexity. The simplest models find a bifurcation associated with the AMOC ‘on’ state losing stability that is a saddle node. Here, we study a physically derived global oceanic model of Wood et al. with five boxes, that is calibrated to runs of the FAMOUS coupled atmosphere-ocean general circulation model. We find the loss of stability of the ‘on’ state is due to a subcritical Hopf for parameters from both pre-industrial and doubled CO 2 atmospheres. This loss of stability via subcritical Hopf bifurcation has important consequences for the behaviour of the basin of attraction close to bifurcation. We consider various time-dependent profiles of freshwater forcing to the system, and find that rate-induced thresholds for tipping can appear, even for perturbations that do not cross the bifurcation. Understanding how such state transitions occur is important in determining allowable safe climate change mitigation pathways to avoid collapse of the AMOC.

The Atlantic Meridional Overturning Circulation (AMOC) transports substantial amounts of heat into the North Atlantic sector, and hence is of very high importance in regional climate projections. The AMOC has been observed to show multi-stability across a range of models of different complexity. The simplest models find a bifurcation associated with the AMOC `on' state losing stability that is a saddle node. Here we study a physically derived global oceanic model of Wood et al with five boxes, that is calibrated to runs of the FAMOUS coupled atmosphere-ocean general circulation model. We find the loss of stability of the `on' state is due to a subcritical Hopf for parameters from both pre-industrial and doubled CO_2 atmospheres. This loss of stability via subcritical Hopf bifurcation has important consequences for the behaviour of the basin of attraction close to bifurcation. We consider various time-dependent profiles of freshwater forcing to the system, and find that rate-induced thresholds for tipping can appear, even for perturbations that do not cross the bifurcation. Understanding how such state transitions occur is important in determining allowable safe climate change mitigation pathways to avoid collapse of the AMOC.

Models incorporating delay have been frequently used to understand climate variability phenomena, but often the delay is introduced through an ad hoc physical reasoning, such as the propagation time of waves. In this paper, the Mori-Zwanzig formalism is introduced as a way to systematically derive delay models from systems of partial differential equations and hence provides a better justification for using these delay-type models. The Mori-Zwanzig technique gives a formal rewriting of the system using a projection onto a set of resolved variables, where the rewritten system contains a memory term. The computation of this memory term requires solving the orthogonal dynamics equation, which represents the unresolved dynamics. For nonlinear systems, it is often not possible to obtain an analytical solution to the orthogonal dynamics and an approximate solution needs to be found. Here, we demonstrate the Mori-Zwanzig technique for a two-strip model of the El Niño Southern Oscillation (ENSO) and explore methods to solve the orthogonal dynamics. The resulting nonlinear delay model contains an additional term compared to previously proposed ad hoc conceptual models. This new term leads to a larger ENSO period, which is closer to that seen in observations.

We present a study of a delay differential equation (DDE) model for the Mid-Pleistocene Transition (MPT). We investigate the behavior of the model when subjected to periodic forcing. The unforced model has a bistable region consisting of a stable equilibrium along with a large amplitude stable periodic orbit. We study how forcing affects solutions in this region. Forcing based on astronomical data causes a sudden transition in time and under increase of the forcing amplitude, moving the model response from a non-MPT regime to an MPT regime. Similar transition behavior is found for periodic forcing. A bifurcation analysis shows that the transition is not due to a bifurcation but instead to a shifting basin of attraction. While determining the basin boundary we demonstrate how one can accurately compute the intersection of a stable manifold of a saddle with a slow manifold in a DDE by embedding the algorithm for planar maps proposed by England et al. (SIADS 2004(3)) into the equation-free framework by Kevrekidis et al. (Rev. Phys. Chem. 2009 (60)).

The codes needed to reproduce the results of the paper. For every figure there is a folder with the corresponding code.

The Mid-Pleistocene Transition, the shift from 41 kyr to 100 kyr glacial-interglacial cycles that occurred roughly 1 Myr ago, is often considered as a change in internal climate dynamics. Here we revisit the model of Quaternary climate dynamics that was proposed by Saltzman and Maasch (1988). We show that it is quantitatively similar to a scalar equation for the ice dynamics only when combining the remaining components into a single delayed feedback term. The delay is the sum of the internal times scales of ocean transport and ice sheet dynamics, which is on the order of 10 kyr. We find that, in the absence of astronomical forcing, the delayed feedback leads to bistable behaviour, where stable large-amplitude oscillations of ice volume and an equilibrium coexist over a large range of values for the delay. We then apply astronomical forcing. We perform a systematic study to show how the system response depends on the forcing amplitude. We find that over a wide range of forcing amplitudes the forcing leads to a switch from small-scale oscillations of 41 kyr to large-amplitude oscillations of roughly 100 kyr without any change of other parameters. The transition in the forced model consistently occurs near the time of the Mid-Pleistocene Transition as observed in data records. This provides evidence that the MPT could have been primarily a forcing-induced switch between attractors of the internal dynamics. Small additional random disturbances make the forcing-induced transition near 800 kyr BP even more robust. We also find that the forced system forgets its initial history during the small-scale oscillations, in particular, nearby initial conditions converge prior to transitioning. In contrast to this, in the regime of large-amplitude oscillations, the oscillation phase is very sensitive to random perturbations, which has a strong effect on the timing of the deglaciation events.