Rate-induced tipping in natural and human systems

Over the last 2 decades, tipping points in open systems subject to changing external conditions have become a topic of a heated scientific debate due to the devastating consequences that they may have on natural and human systems. Tipping points are generally believed to be associated with a system bifurcation at some critical level of external conditions. When changing external conditions across a critical level, the system undergoes an abrupt transition to an alternative, and often less desirable, state. The main message of this paper is that the rate of change in external conditions is arguably of even greater relevance in the human-dominated Anthropocene but is rarely examined as a potential sole mechanism for tipping points. Thus, we address the related phenomenon of rate-induced tipping: an instability that occurs when external conditions vary faster, or sometimes slower, than some critical rate, usually without crossing any critical levels (bifurcations). First, we explain when to expect rate-induced tipping. Then, we use three illustrative and distinctive examples of differing complexity to highlight the universal and generic properties of rate-induced tipping in a range of natural and human systems.

1 Introduction

In this paper, we consider tipping instabilities in nonlinear open systems (Ashwin et al., 2012). By open, we mean systems that are influenced by changing external conditions which we refer to as external forcings. In a mathematical dynamic model of an open system, such external forcings are represented by time-varying input parameters.

Figure 1 Schematic illustration of threshold instability. (a) The stability landscape of the system at an initial level of forcing. The well represents the base state, the hill top defines the threshold (indicated by the vertical dashed red line), and the (blue) ball indicates the current state of the system. To the left of the hill top (threshold), the ball rolls into the well, meaning the system converges to the base state. To the right of the hill top, the ball runs away, indicating tipping to an alternative state. (b) The stability landscape at a new forcing level. Note that the initial base state is on the other side of the new hill top (threshold). Therefore, if the system is at the initial base state and the forcing switches sufficiently fast to the new level, the ball will run away, meaning the system will tip to an alternative state.

Large and abrupt changes in the state of an open system may occur when the external forcing exceeds some critical level (Scheffer, 2010; Lenton, 2011; Kuehn, 2011). The points in time, or in the level of forcing, at which such changes occur are commonly referred to as bifurcation-induced tipping points (Ashwin et al., 2012). They have been identified in many domains, including ecosystems (Scheffer et al., 1993, 2001, 2009; Siteur et al., 2014; Dakos et al., 2019; Pierini and Ghil, 2021) and the human brain (Rinzel and Ermentrout, 1998; Moehlis, 2008; Screen and Simmonds, 2010; Mitry et al., 2013; Maturana et al., 2020), and are of particular concern under anthropogenic climate change (Lenton et al., 2008; Ashwin and von der Heydt, 2020; Arias et al., 2021; Ritchie et al., 2021; Boers and Rypdal, 2021; Boulton et al., 2022). Furthermore, it has recently been recognised that critical levels can be exceeded temporarily without causing tipping (van der Bolt et al., 2018; Ritchie et al., 2019; Alkhayuon et al., 2019; O'Keeffe and Wieczorek, 2020). This occurs when the time of exceedance is short compared to the inherent timescale of the system (O'Keeffe and Wieczorek, 2020; Ritchie et al., 2021; Alkhayuon et al., 2023).

However, there is another, less obvious potential consequence of changes in external forcing. When an external forcing changes faster than some critical rate rather than necessarily by a large amount, this can lead to rate-induced tipping points (Stocker and Schmittner, 1997; Luke and Cox, 2011; Wieczorek et al., 2011; Ashwin et al., 2012; Ritchie and Sieber, 2016; Siteur et al., 2016; Suchithra et al., 2020; Arumugam et al., 2020; Pierini and Ghil, 2021; Wieczorek et al., 2023; Longo et al., 2021; Kuehn and Longo, 2022; Kaur and Sharathi Dutta, 2022; Hill et al., 2022; Arnscheidt and Rothman, 2022). In contrast to bifurcation-induced tipping, rate-induced tipping occurs due to fast-enough changes in external forcing and usually does not exceed any critical levels as a result of external forcing. Such tipping points are much less widely known and yet are arguably even more relevant to contemporary issues such as climate change (Lohmann and Ditlevsen, 2021; Clarke et al., 2021; O'Sullivan et al., 2022), ecosystem collapse (Scheffer et al., 2008; Vanselow et al., 2019; van der Bolt and van Nes, 2021; Neijnens et al., 2021; Vanselow et al., 2022), and the resilience of human systems (Witthaut et al., 2021).

This paper combines a review writing style with new results to make the concept of rate-induced tipping points accessible to a wide scientific audience. Even though the phenomenon is rarely discussed by scientists and policy makers, we argue that it is ubiquitous and likely to be prevalent in many open systems. The rationale is that the current human-dominated era of Earth history, which has been called the Anthropocene (Crutzen, 2002; Crutzen and Stoermer, 2021), is characterised by systems (e.g. climate, ecosystems, infrastructures, and economy) that are subject to fast-changing external conditions, and are thus kept far from the changing equilibrium as a result of human activity. These circumstances of rapidly changing external forcing are precisely the conditions that can lead to rate-induced tipping. Rate-induced tipping is therefore especially relevant to the contemporary period, even though the phenomenon is not widely known or understood. By contrast, bifurcation-induced tipping in a system that is forced slowly towards a threshold and that thus stays close to or tracks the changing equilibrium (Ashwin et al., 2012) is much more widely understood but is less relevant in a rapidly forced system (Ritchie et al., 2021). We demonstrate this through an analysis of the following three distinct dynamic models of natural and human systems of differing complexity: a predator–prey ecosystem, the large-scale ocean circulation, and an electrical power grid network.

Rate-induced tipping occurs when the system deviates too much from the changing stable equilibrium and crosses some threshold. Here, we focus on examples of what Wieczorek et al. (2023) call a regular threshold. Easily verifiable criteria for the occurrence of rate-induced tipping, such as threshold or basin instability, are identified in all of these examples. Furthermore, we uncover universal features of rate-induced tipping. These include multiple critical rates of change due to the interaction of different timescales of the external forcing with the inherent timescales of the system. Finally, we highlight important phenomena, such as return tipping, that are non-obvious and can be easily overlooked. For further reading on different threshold types, see O'Sullivan et al. (2022) for an example of an elusive quasithreshold and Lohmann and Ditlevsen (2021) for an example of what appears to be a fractal-like irregular threshold.

2 When to expect rate-induced tipping

A system is known to be susceptible to rate-induced tipping if the state the system currently resides in is threshold unstable (O'Keeffe and Wieczorek, 2020; Wieczorek et al., 2023). Suppose that, for a given initial level of external forcing, the system resides in a stable equilibrium (although, in general, it can be a stable limit cycle or an even more complicated attractor). This equilibrium will be referred to as the base state. One way of depicting threshold instability is with a moving-stability landscape, as illustrated by Fig. 1, where the base state is represented by the blue ball in the well in Fig. 1a¹. The hill top defines the position of the threshold (indicated by the vertical dashed red line). If the ball is to the left of the hill top (threshold), it will roll into the well, and the system will converge to the base state, whereas if the ball is to the right of the hill top, it will roll in the opposite direction, and the system will tip to some alternative state. The alternative state may be a different stable state for a multi-stable system (Scheffer et al., 2008; O'Keeffe and Wieczorek, 2020; Halekotte and Feudel, 2020; Lohmann et al., 2021; Slyman and Jones, 2023) or a transient state for an excitable (possibly mono-stable) system (Wieczorek et al., 2011; Vanselow et al., 2019; O'Sullivan et al., 2022).

In contrast to bifurcation-induced tipping, the change in the forcing usually does not cause any qualitative change in the stability landscape but instead shifts its position. If the threshold moves past the initial position of the base state for a new forcing level, as shown in Fig. 1b, the base state is said to be threshold unstable in terms of varying the forcing (Wieczorek et al., 2023). In the case when the threshold is a basin boundary of two attractors in a multi-stable system, the system is said to be basin unstable (O'Keeffe and Wieczorek, 2020). The threshold or basin instability condition gives the forcing shift magnitude that enables rate-induced tipping. In general, one can prove that threshold or basin instability is sufficient for the occurrence of rate-induced tipping: there is an external forcing that gives rate-induced tipping if the system is threshold or basin unstable (Kiers and Jones, 2020; Wieczorek et al., 2023). In many examples, including those considered here, we find that threshold instability appears to be both necessary and sufficient for the occurrence of rate-induced tipping: there is an external forcing that gives rate-induced tipping if and only if the system is threshold or basin unstable.

The rigorous result can be understood intuitively as follows. Consider a change in the level of the forcing that gives threshold or basin instability as depicted in Fig. 1. If the forcing changes from the initial level to the new level at a sufficiently slow rate, the ball remains in the well, and the system is said to track the moving base state. If the forcing switches at a sufficiently fast rate, the initial ball finds itself on the other side of the hill top (threshold) and tips to an alternative state. Thus, there will be at least one intermediate critical rate of change at which there is a transition between tracking and tipping. Once it is known that a system is threshold or basin unstable and thus susceptible to rate-induced tipping, the goal is to find the critical rate, or even multiple critical rates, for a given profile (shape) of external forcing. In the next section, we give a more precise description of critical rates.

3 Defining critical rates

Let us denote the time-varying external forcing with λ. The level of the forcing at a time t is simply the value of λ at this time t. However, defining critical rates of change in external forcing is more subtle. On the one hand, different external forcings will have different physical units and be different, often nonlinear, functions of time. On the other hand, we would like to quantify critical rates of change in a uniform way that is independent of the physical units and the temporal profile of the forcing. Therefore, we introduce a rate parameter r in units per inverse second (or day, year, etc.); we write the external forcing as λ(rt), where u=rt is dimensionless; and we work with r as the main input parameter. Most importantly, we define a critical rate as a special value of r at which rate-induced tipping occurs, while the shift magnitude of λ(rt) remains fixed.

To avoid confusion between the rate parameter r and the rate of change of external forcing $\mathrm{d}\mathit{\lambda }/\mathrm{d}t$, we note that

$${\displaystyle \frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}u}}{\displaystyle \frac{\mathrm{d}u}{\mathrm{d}t}}=r{\displaystyle \frac{\mathrm{d}\mathit{\lambda }\left(u\right)}{\mathrm{d}u}}$$

has units of λ per second, depends on r and on the profile of λ(u), and may itself be a function of time. In other words, the rate parameter r quantifies the rate of change of external forcing with a given profile. Furthermore, if the forcing λ itself is a physical rate of some sort (e.g. freshwater flux into the North Atlantic, measured in Sverdrups (millions of m³ s⁻¹) or population growth rate, measured in individuals per unit area per year, as per the examples in Sect. 5), λ(rt) will be the level of this rate at time t, referred to as the level of the forcing, and r will quantify the rate of change of this rate, referred to as the rate of change of the forcing.

Figure 2 Illustration of bifurcation-induced, rate-induced, and return tipping. (a–c) Time profiles of ramp external forcing and (d–f) corresponding system response to these forcing profiles. (a, d) Bifurcation-induced tipping – a slow change in external forcing past a critical level (fold bifurcation) causes tipping. The tipping occurs for any rate of change of the forcing past the fold. (b, e) Rate-induced tipping – the system fails to adapt to a too-fast change in external forcing (red) even though the forcing never crosses the critical level. For a slow-enough change of the forcing, the system tracks the moving base state and avoids tipping (green). (c, f) Return tipping – avoiding bifurcation-induced tipping by reversing the trend in the external forcing too quickly can lead to rate-induced tipping on a decrease of the forcing. Branches of stable equilibria are denoted by solid black curves, and branches of unstable equilibria are denoted by dashed black curves. Stable and unstable branches meet at a fold bifurcation (black dot). The system starting from the grey dot is basin unstable for forcing shift magnitudes that end in the grey region of basin instability.

4 Rate-induced tipping in a simple model

In natural and human systems, tipping points are often associated with crossing a critical level of the forcing, defined by a dangerous (e.g. fold) bifurcation for the frozen system with fixed-in-time forcing, causing a catastrophic, abrupt, and irreversible change to the state of the system (Thompson et al., 1994; Thompson and Sieber, 2011). This type of tipping is commonly referred to as bifurcation-induced tipping (or B-tipping) and is illustrated by Fig. 2a and d (Ashwin et al., 2012). Suppose the system starts from the base state on the upper branch of stable equilibria for the frozen system, as indicated by the grey dot in Fig. 2d. Initially, as external forcing changes slowly (Fig. 2a), the state of the system (the red trajectory in Fig. 2d) tracks the moving base state (the branch of stable equilibria). However, once external forcing reaches the fold that defines the critical level of the forcing, the base state disappears (the branch of stable equilibria terminates), and the system subsequently undergoes a catastrophic transition to the alternative stable state. Crucially, the tipping occurs for any rate of change of the forcing past the fold. The alternative state is often a less desired state, such as an extinction state in an ecosystem (O'Keeffe and Wieczorek, 2020), collapse of an ocean circulation (Alkhayuon et al., 2019), or a blackout on a power grid network (Budd and Wilson, 2002). However, it could also be a more desired state, such as a well-being state for developing countries (Mirza et al., 2019).

Figure 2 introduces a subtle but crucial difference to previous examples that have considered B-tipping (Lenton et al., 2008; Scheffer et al., 2009; Ritchie et al., 2021), and that is to apply a tilt² to the bifurcation structure (O'Keeffe and Wieczorek, 2020, Sect. 7); see the methods section in Appendix A for further details. This important distinction introduces the possibility of a different form of tipping known as rate-induced tipping (or R-tipping). Unlike B-tipping in Fig. 2a and d, where crossing a critical level of the external forcing causes a catastrophic transition, R-tipping occurs when the system fails to adapt to a too-rapidly-changing external forcing, usually without crossing any critical levels. Figure 2b and e consider two scenarios where the change in the level of external forcing is the same but occurs at different rates. Most importantly, the forcing stops in the (grey) region of basin instability and never crosses the critical level defined by the fold bifurcation point. For a slow rate of change in external forcing (green trajectory), the system is able to continually adapt to the moving base state and tracks the stable branch of equilibria without tipping. However, for a slightly faster rate of change in external forcing (red trajectory), the system is unable to adapt to the moving base state and undergoes R-tipping to the alternative stable state.

If a system is thought to be approaching a B-tipping event, then a natural option would be to reverse the external forcing to avoid crossing a largely unknown critical level. However, fast reversals in the forcing could introduce a new problem that has been largely overlooked, namely return tipping (O'Keeffe and Wieczorek, 2020). Figure 2c and f illustrate such a scenario for a fold bifurcation structure tilted down. Suppose that external forcing has caused a system to approach close to the fold bifurcation. Reversing the forcing slowly will allow the system to closely track the moving base state (the branch of stable equilibria) as shown by the green trajectory. However, a too-fast reversal may give rise to R-tipping on return if the system is basin unstable on reversal of the forcing. Then, the end result is opposite to what was intended. Although B-tipping is avoided, the system, rather surprisingly, R-tips to the alternative stable state (red trajectory). Therefore, in general, reversing external forcing as quickly as possible does not guarantee avoiding tipping.

Figure 3 Tipping for a nonlinear shift ramp and return profiles in the tilted-fold example. (a) Time profiles of ramp and return (impulse) external forcings and (b) system response to these forcing profiles. Forcing profiles vary between the same minimum and maximum levels but at different rates: slow (blue), medium (orange), and fast (purple). Ramp forcing profiles in (a) are given by a concatenation of the left half of a colour curve and the dashed black curve. The corresponding system responses in (b) are given by a concatenation of the left half of the solid curve and the dashed curve of the same colour. Return forcing profiles in (a) and the corresponding system responses in (b) are given by the solid colour curves. (c) Tipping diagram for ramp and return profiles; note the logarithmic scale for the rate parameter. Critical boundaries separate regions of tracking from tipping for the ramp forcing profile (black dashed curve) and for the return forcing profile (solid black curve). White region – tracking for ramp and return profiles; green region – tipping for ramp profile, tracking for return profile; red region – tipping for ramp and return profiles.

Figure 3 provides a more in-depth analysis of the tilted saddle-node model considered in Fig. 2d and e. Let us assume initially that the external forcing has the profile of a nonlinear shift ramp; see the methods section in Appendix A for details. Three sample time series of shift forcing profiles between the same levels but at different rates are given by a concatenation of the left half of a colour curve and the dashed black curve in Fig. 3a; these are similar to external forcings used in Fig. 2a and b. The corresponding response of the system is depicted in Fig. 3b. For the slowest change in external forcing (the dashed blue trajectory), the system is able to adapt to and track the changing base state. However, if the rate of change in external forcing becomes too fast (the dashed orange and purple trajectories), then the system fails to adapt to the changing base state and R-tips to the alternative state.

The critical rate, which determines the onset of R-tipping, will depend on how much the external forcing is changed by. The dashed black curve in Fig. 3c shows the critical boundary for the ramp external forcing, separating regions of tipping (coloured) from no tipping (white), in the plane of the rate parameter against the change in the forcing level, referred to as the peak change. B-tipping occurs if the external forcing crosses the fold bifurcation without returning. Indeed, for very small rate parameters (slow rates), the critical boundary asymptotes to the distance required to reach the fold (indicated by the thin black line). However, for larger rate parameters (faster rates), tipping can occur before the fold is transgressed because of R-tipping. The three coloured dots correspond to the external forcing parameters used in Fig. 3a and b. Notice that the blue dot is in the white region, signifying tracking, whereas the other two dots are within the coloured regions, denoting tipping for the ramp forcing. For very large rate parameters, the critical boundary asymptotes to the basin instability boundary – the smallest change in the level of external forcing that gives basin instability, as described in Fig. 1. In summary, the ramp forcing with a peak change past the basin instability boundary and below the fold level gives rise to R-tipping with a single critical rate. This critical rate decreases with the peak change.

The dynamics become more interesting when the external forcing is reversed back to its initial level after reaching its peak level. To illustrate how, we will now consider external forcing that has a profile of a symmetric impulse, referred to as a return forcing profile³; see the methods section in Appendix A for further details. Three sample time series of return forcing profiles with the same peak change but different rates are given by the colour curves in Fig. 3a. The critical boundary for the return external forcing is given by the solid black curve in Fig. 3c, which separates regions of tracking (green and white) from tipping (red) and is very different from the dashed curve for the ramp external forcing. To be more specific, the green (points of return) region corresponds to scenarios where tipping is prevented by reversing the external forcing. The red (points of no return) region corresponds to scenarios where tipping still occurs despite reversing the external forcing. Previously, it has been shown for B-tipping that safely overshooting a critical level by a given distance can be achieved, provided the reversal in external forcing is faster than some critical rate. However, the added possibility of R-tipping owing to the tilted bifurcation structure, combined with the symmetric return forcing (see the methods in Appendix A section for further details), means that multiple critical rates can arise for a fixed peak change in the return forcing profiles; this is illustrated by the S-shaped solid black curve in Fig. 3c. In a symmetric return forcing, multiple critical rates emerge because there is competition between the sufficiently slow approach towards the fold required to avoid R-tipping and the sufficiently fast reversal required for safe overshoots of the fold.

In the green region to the left of the vertical fold line, reversing the forcing prevents the system from an impending R-tipping that would occur if the forcing were not to be reversed. An example is given by the solid orange curve in Fig. 3b. Interestingly, there is also a small red region to the left of the vertical fold line. This region gives rise to two critical rates for a fixed peak change, which bound a (red) sub-interval of the rate parameter where R-tipping is not prevented by return forcing. An example is given by the solid purple curve in Fig. 3b.

For small overshoots of the fold, even greater complexity is possible with the potential of three critical rates and two (red) tipping sub-intervals for a fixed peak change of return forcing. For very small rate parameters (the red region to the right of the lower part of the vertical fold line), B-tipping occurs and cannot be prevented by reversing the forcing. However, for the same peak change and larger rate parameters (the slightly wider green region to the right of the lower part of the vertical fold line), it becomes possible to prevent B-tipping and avoid R-tipping upon return. Keeping the peak change fixed and increasing the rate parameter even more (the red region to the right of the middle part of the vertical fold line) prevents B-tipping but triggers R-tipping, meaning that the system tips again despite reversing the forcing. Then, for the same peak change and very large rate parameters (the green region to the right of the upper part of the vertical fold line), both B-tipping and R-tipping upon return can be prevented again but for a different reason – the system processes are too slow to react to a fast forcing impulse.

5 Rate-induced tipping in ecology and climate

We now consider an example from ecology, namely that of a predator–prey system, which models the time evolution of plant and herbivore biomass densities (Scheffer et al., 2008). The model has been proposed to conceptually study tipping points in bistable ecosystems with a non-monotone functional response. Examples of such systems can be the dominance shift between submerged macrophytes and phytoplankton (Scheffer et al., 1993) and between coral reefs and macro-algae (Hughes, 1994) or the transition of kelp forests into sea urchin barrens that are dominated by crustose coralline algae (Steneck et al., 2002).

In this example, changes in environmental conditions affect the plant growth rate and herbivore mortality rate simultaneously and play the role of external forcing. Such a scenario can be considered to be possible under climate change, where plants benefit from the fertilisation effect due to increasing levels of CO₂ (Reich et al., 2014) but where the resulting increased temperatures are detrimental to herbivores (Lacetera, 2019). The plant growth rate and herbivore mortality rate vary within a range where the ecosystem has two stable equilibria. The stable coexistence equilibrium is the base state. The stable plant-only equilibrium with no herbivores is the alternative stable state.

Figure 4 Rate-induced tipping in the plant–herbivore (a–c) and AMOC (d–f) models. Time series of the external ramp forcing profiles – (a) dimensionless environmental conditions and (d) freshwater hosing; see the methods section in Appendix A for further details. Time series of the system responses to the changing (b) environmental conditions and (e) freshwater hosing for three different rates of change (different colours). The system responses to external ramp forcing profiles in the phase plane of (c) the plant and herbivore biomass and (f) the North and Tropical Atlantic salinities. See supplementary videos “PH_Rtipping_full.mp4” and “AMOC_Rtipping_full.mp4” (Ritchie et al., 2023) for animated versions of (a–c) and (d–f), respectively.

We consider three sample time profiles of a nonlinear shift ramp external forcing, shown in Fig. 4a, with the same change in the level of environmental conditions but at different rates. The resulting impacts on the herbivore biomass are shown in Fig. 4b. For the slowest (green) change in environmental conditions, the ecosystem tracks the moving base state. Herbivores increase slightly, which represents the continued presence of coexisting plants and herbivores. Contrast this to the fastest (red) change in environmental conditions, which causes R-tipping to the alternative stable state. In this scenario, the herbivore population declines to zero, and the ecosystem becomes entirely dominated by plants.

Figure 4c illustrates the underlying dynamics in the phase plane of the plant and herbivore biomass. For the initial level of environmental conditions, the base state is indicated by the grey dot. The shift in the level of environmental conditions changes the position of the base state, and this change is indicated by the dotted black line. The black dot at the other end of the dotted line indicates the base state for the final level of environmental conditions. Further ramifications of the shift in environmental conditions include changes in the basin of attraction of the base state. The basin of attraction shifts (from the dotted curve to the solid grey curve) such that the initial base state (the grey dot) is not contained in the basin of attraction of the final base state (the region above the solid grey curve). A consequence of basin instability is that the behaviour of solutions starting near the initial base state depends on the rate of change of environmental conditions. For slow rates, the system is able to continually adapt and remain within the changing basin of attraction of the base state so that solutions converge to the final base state (the green trajectory). However, for a sufficiently fast rate, the ecosystem fails to track the fast moving base state and R-tips to the final alternative stable state indicated by the second black dot (the red trajectory). This rate-sensitive behaviour is expected due to basin instability. What may be surprising is that, at a critical rate of change in environmental conditions, the corresponding solution converges to an unstable edge state (black circle) (Wieczorek et al., 2023) on the basin boundary of the final base state (the blue trajectory).

Next, we consider an example from climate, specifically the possible collapse of the Atlantic meridional overturning circulation (AMOC) under global warming. The AMOC forms part of the global thermohaline circulation, which is a large-scale ocean circulation current driven by temperature and salinity gradients. The AMOC contributes to the relatively mild climate in western Europe by transporting heat from the tropics to the North Atlantic. Once the warm salty waters reach the North Atlantic, they cool down and become denser. The higher density of these waters causes sinking (or overturning), followed by a return to the tropics along the bottom of the ocean. However, the AMOC can be easily disturbed by contemporary climate change. The amount of freshwater added to the North Atlantic (referred to as freshwater hosing) may increase under climate change, for example due to the melting of the Greenland Ice Sheet, changes in precipitation patterns, or both. Specifically, Stocker and Schmittner (1997) and Lohmann and Ditlevsen (2021) show in coupled climate models that the AMOC can collapse under sufficiently fast rates of change in either CO₂ emissions or freshwater hosing. Additionally, R-tipping of the AMOC has been observed in a global oceanic box model (Alkhayuon et al., 2019).

Here, we work with the global oceanic box model for the AMOC (Wood et al., 2019), in which changing freshwater hosing plays the role of external forcing. The forcing varies within a range where the AMOC model has two stable equilibria. The stable AMOC-On equilibrium is the base state. The stable AMOC-Off equilibrium is the alternative stable state.

Figure 5 Tipping diagrams for (a) the plant–herbivore and (b) AMOC models. Critical boundaries separate regions of tracking from tipping for the ramp forcing profile (dashed black curve) and for the return forcing profile (solid black curve). White region – tracking for ramp and return forcing profiles; green region – tipping for ramp forcing profile, tracking for return forcing profile; red region – tipping for ramp and return forcing profiles.

Figure 4d shows three sample time profiles of ramp forcing that increase the freshwater hosing from zero to the same non-zero level but each at a different rate. The response of the AMOC to these freshwater-forcing scenarios is shown in Fig. 4e. For the slowest (green) change in freshwater hosing, there is tracking of the moving base state. The AMOC strength suffers a slight drop but ultimately remains in the AMOC-On state. However, for the fast (red) change in freshwater hosing, there is R-tipping to the alternative stable state. The AMOC strength declines to a point of complete collapse.

The phase portrait in the plane of the salinities in the North and Tropical Atlantic boxes in Fig. 4f illustrates the underlying dynamics. The grey dot shows the base state for the initial level of freshwater hosing. Increasing the freshwater hosing shifts the base state as indicated by the dotted black line. The black dot at the other end of the dotted line indicates the base state at the final level of freshwater hosing. Notice that the basin of attraction of the final base state (the region enclosed by the blue periodic orbit) does not contain the initial base state (the grey dot). Therefore, the initial base state is basin unstable, and R-tipping from the base state to the alternative stable state (lower black dot) will occur for sufficiently fast shifts in the level of freshwater hosing. For slow rates of increase in freshwater hosing, the system continually adapts and remains within the changing basin of attraction of the base state so that solutions converge to the final base state (the green trajectory). On the other hand, for sufficiently fast rates of increase in freshwater hosing, the system is unable to adapt and falls outside the changing basin of attraction of the base state, causing solutions to converge to the final alternative state (the red trajectory). At the critical rate of increase in freshwater hosing, surprising behaviour is again observed as the corresponding solution converges to a repelling periodic orbit that defines the basin boundary of the final base state (the blue trajectory).

To validate our choice of the simple tilted saddle-node model in Sect. 4, we plot in Fig. 5 the tipping diagrams for the plant–herbivore and AMOC models, showing regions of tracking (white), points of return (green), and points of no return (red). The tipping diagrams for the plant–herbivore model in Fig. 5a and the AMOC model in Fig. 5b are very similar to the tipping diagram for the simple model in Fig. 3b, although there are some differences. In both examples, there is a region of basin instability that gives rise to R-tipping for shift forcing profiles with a peak change that does not cross any critical levels. This means that the simple model with a tilted bifurcation structure indeed captures non-obvious tipping phenomena found in higher-dimensional systems. In Fig. 5a, the (red) region of points of no return extends to the left of the critical level (the vertical black fold line), but there are only up to two critical rates for the return forcing profile with a fixed peak change, and the (green) region of points of return vanishes for small rate parameters. In Fig. 5b, there are up to three critical rates for the return forcing profile with a fixed peak change, but the (red) region of points of no return does not extend to the left of the critical level (the vertical black Hopf line). Additionally the (black) boundary has small smooth wiggles that appear as non-smooth corners. A similar wiggling effect near a Hopf bifurcation has been observed in O'Keeffe and Wieczorek (2020).

Figure 6 Tipping for ramp and return forcing profiles in the power grid model. (a) Time profiles of ramp and return (impulse) power demand and (b) the resulting voltage response for these power demand (forcing) profiles. Power demand profiles vary between the same minimum and maximum levels but at different rates: slow (blue), medium (orange), and fast (purple). Ramp forcing profiles are given by a concatenation of the left half of a colour curve and the dashed black curve. (b) The corresponding system responses given by a concatenation of the left half of the solid curve and the dashed curve of the same colour. Return forcing profiles in (a) and the corresponding system responses in (b) are given by the solid colour curves. (c) Tipping diagram of ramp and return forcing profiles. Critical boundaries separate regions of tracking from tipping (blackout) for both the ramp profile (dashed black curve) and the return profile (solid black curve). For return forcing profiles, reversible tipping is also possible through a phase slip. White region – tracking for ramp and return profiles; green region – irreversible tipping to blackout for ramp profiles, tracking for return profiles; light-red region – irreversible tipping to blackout for ramp profiles, reversible tipping to phase slips for return profiles; dark-red region – irreversible tipping to blackout for ramp and return profiles.

6 Rate-induced tipping in power grid networks

R-tipping instabilities are not confined to natural systems but can occur in any system, including human systems. One example is the energy sector and power grid networks (Suchithra et al., 2020). Crucially, electricity needs to be used as soon as it is produced, since it cannot be stored easily (Mokrian and Stephen, 2006). Therefore, providing a near-constant voltage of electricity to millions of homes while power demand varies seasonally and daily and may spike during major events is a technological challenge (Mahmud and Zahedi, 2016). Some noticeable examples of power blackouts or near misses are the northeastern USA blackout in 2003, caused by a series of faults in local control systems (Hu et al., 2016), and the near-miss blackout in England following the conclusion of the Euro 1990 semi-final (Swarup, 2007).

The latter example in particular was arguably the result of R-tipping effects. The power demand on the network, following the conclusion of the football match, was expected to be high. Hence, the national grid took measures to ensure that the network would be able to cope with the high power demand. However, the national grid failed to envisage the match going to extra time and penalties. Thus, the rapid increase in power demand, following the eventual conclusion of the match, gave controllers insufficient time to react. Therefore, in this case, the limiting factor was not the peak in demand but instead the rate at which the demand on the network rose.

Here, we use a conceptual model of a power grid network (Dobson et al., 1988; Dobson and Chiang, 1989; Budd and Wilson, 2002), where changing power demand plays the role of external forcing. The situation with stable states in this model is more complicated than in the previous examples. The power demand changes within a range where there are infinitely many stable equilibria with the same fixed voltage magnitude and a 2π difference in the phase angle. Since each of these stable equilibria has the same voltage magnitude, they can be thought of as a single base state of the system. Furthermore, there are two alternative states. An alternative transient state is a temporary drop in the voltage magnitude accompanied by a 2π shift in the phase angle, caused by the coupling within the system. This state corresponds to a transition between two neighbouring stable equilibria within the base state. An alternative stable state is at zero voltage magnitude and corresponds to an electrical blackout⁴.

We now demonstrate that a rapid increase in the level of power demand can lead to R-tipping in the form of different disruptions in power supply. First, we note that, for slow-enough increases in the level of power demand, the power grid network always tracks the moving base state (not shown here). Then, in Fig. 6, we show the response of the power grid network to ramp and return forcing profiles with higher rates of change.

We start with three ramp shifts that vary between the same levels of power demand without crossing any critical levels but at different rates. These are given by a concatenation of the left half of a colour curve and the dashed black curve in Fig. 6a. Owing to the higher rates of change, all three ramp shifts cause R-tipping to the alternative stable state, resulting in blackout; see the dashed curves in Fig. 6b. However, reversing the power demand (see solid colour pulses in Fig. 6a) can avoid blackout and restore the power grid network to the base state if the reversal is fast enough (see the corresponding solid colour responses in Fig. 6b). The slowest (blue) reversal is too slow to avoid irreversible R-tipping to blackout. The medium-rate (orange) reversal avoids blackout but involves reversible R-tipping to the alternative transient state – a temporary voltage drop followed by a long recovery towards the base state. The fastest (purple) reversal avoids both types of R-tipping and does not cause any noticeable disruptions to power supply.

The tipping diagram for ramp and return forcing profiles in Fig. 6c has, in addition to multiple critical rates already observed in the previous examples, an important new feature. Owing to the existence of two alternative states, there are two different red regions corresponding to different types of R-tipping for a return forcing profile: irreversible R-tipping to blackout (dark-red region), where the voltage magnitude drops to zero permanently, and reversible R-tipping (light-red region) involving a temporary drop in the voltage magnitude accompanied by a 2π phase slip. The three scenarios illustrated in Fig. 6a and b correspond to the three coloured dots in the tipping diagram in Fig. 6c. Note that all three dots are located past the basin instability boundary and below the critical level of power demand (the vertical black fold line), meaning that these transitions are purely rate induced.

7 Conclusions

We have shown that many natural and human systems can experience rate-induced tipping (R-tipping). Such instabilities usually occur for sufficiently fast increases in external forcing despite never crossing any critical levels of external forcing. In other words, systems are able to continually adapt to a moving base state and therefore avoid tipping when external forcing changes sufficiently slowly but fail to adapt to or track a moving base state when external forcing changes faster than some critical rate. Reversing the external forcing can prevent a system from suffering R-tipping, but the rates required for this add an additional layer of complexity in the presence of bifurcation-induced tipping (B-tipping). Previously, it has been shown that safe overshoots of critical levels for B-tipping require fast rates of change (Ritchie et al., 2021). However, faster rates of change make a system more susceptible to R-tipping. To make the concept of R-tipping accessible to a wide scientific audience, we

• described an easily verifiable criterion of threshold or basin instability for R-tipping to occur

• demonstrated basin instability and ensuing R-tipping in conceptual models of natural and human systems, including irreversible and reversible R-tipping

• highlighted interesting phenomena, such as multiple critical rates and return tipping, that can arise from an interplay between R-tipping and B-tipping for non-monotone forcing profiles.

The dynamic models used in this study for representing a predator–prey ecosystem, the Atlantic overturning circulation, and an electrical power grid network are relatively simple. Since our forcing profiles are very idealised by design, we focus on the qualitative behaviour that can arise for different rates of forcing rather than on quantitative predictions. For quantitative predictions, further research on R-tipping is required in more-realistic higher-complexity models, such as state-of-the-art global circulation models, and with more-realistic forcing profiles. The base state in such models may not necessarily be a steady state (an equilibrium) but, for example, could take the form of a periodic orbit or even a chaotic state. This could lead to more complex tipping behaviour, such as phase tipping (P-tipping) (Alkhayuon and Ashwin, 2018; Kaszás et al., 2019; Alkhayuon et al., 2021; Ashwin and Newman, 2021; Alkhayuon et al., 2023).

R-tipping is likely to be prevalent in many systems given contemporary rates of change such as unprecedented anthropogenic climate change. This paper highlights the importance of considering how fast external forcing is changing as opposed to solely focusing on levels of change. Consequently, the actions taken to control the rate of change in forcing are equally as important as the actions taken to control the level at which forcing is halted.

Appendix A: Methods

A1 Forcing profiles

In our analysis, we consider two types of dimensionless nonlinear external forcing profiles, denoted a(rt). The first is a ramp forcing profile that starts close to 0 and subsequently increases continuously until reaching a peak level of 1 at time T. This forcing profile subsequently remains at 1.

$$\begin{array}{}\text{(A1)}& a\left(rt\right)=\left\{\begin{array}{ll}\text{sech}\left(r\right(t-T\left)\right),& \mathrm{0}\le t\le T,\\ \mathrm{1},& t>T.\end{array}\right.\end{array}$$

For R-tipping scenarios with a fixed peak change, the rate of change in the forcing, quantified by the rate parameter r, determines if tipping occurs.

The second forcing profile we consider is a return (impulse) forcing, primarily used to examine the possibility of avoiding tipping. Equation (A1) is modified by removing the piecewise element of the forcing for t>T, such that the forcing returns to its initial level (at a mirrored rate of the approach) after reaching the peak level of 1 at time T:

$$\begin{array}{}\text{(A2)}& a\left(rt\right)=\text{sech}\left(r\right(t-T\left)\right).\end{array}$$

One example of a return forcing profile, given by Eq. (A2), is an idealised scenario to reverse the impact of anthropogenic climate change back to initial levels, i.e. via the development of technologies to remove CO₂ from the atmosphere (Huntingford et al., 2017).

A2 Conceptual model

We use a conceptual model to illustrate some of the universal features associated with R-tipping alone and with an interplay between R-tipping and B-tipping. Using a modified (tilted) version of the normal form for a saddle-node bifurcation, a state variable x is modelled by the following single ordinary differential equation:

$$\begin{array}{}\text{(A3)}& {\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}=-x\left[(x-A-s\mathit{\lambda }(rt){)}^{\mathrm{2}}+\mathit{\lambda }(rt)\right],\end{array}$$

where s is a tilt parameter, A provides the distance between the fold point and the alternative state, and the external forcing is given by

$$\mathit{\lambda }\left(rt\right)={\mathit{\lambda }}_{-}+{\mathrm{\Delta }}_{\mathit{\lambda }}a\left(rt\right),$$

where a(rt) takes the form of a ramp or a return profile defined above. Figure 2a and d are obtained using the parameter values T=500, s=4, A=3.2, ${\mathit{\lambda }}_{-}=-\mathrm{0.5}$, Δ_λ=0.505, and r=0.05 and the initial condition at the location of the base state for $\mathit{\lambda }={\mathit{\lambda }}_{-}$. Figure 2b and e are obtained using the parameter values T=500, s=4, A=3.2, ${\mathit{\lambda }}_{-}=-\mathrm{0.5}$, Δ_λ=0.4, r=1 (green trajectory), and r=3 (red trajectory) and the initial condition at the location of the base state for $\mathit{\lambda }={\mathit{\lambda }}_{-}$. Figure 2c and f are obtained using the parameter values T=500, $s=-\mathrm{4}$, A=1.5 ${\mathit{\lambda }}_{-}=-\mathrm{0.1}$, ${\mathrm{\Delta }}_{\mathit{\lambda }}=-\mathrm{0.4}$, r=1 (green trajectory), and r=3 (red trajectory) and the initial condition at the location of the base state for $\mathit{\lambda }={\mathit{\lambda }}_{-}$. Figure 3 was obtained using a value of T where a(0) is no larger than 10⁻⁴; other parameter values, namely s=4, A=3.2, ${\mathit{\lambda }}_{-}=-\mathrm{0.5}$, r=0.5 (blue trajectory), r=0.9 (orange trajectory), and r=1.7 (purple trajectory); and the initial condition at the location of the base state for $\mathit{\lambda }={\mathit{\lambda }}_{-}$.

A3 Plant–herbivore model

The time evolution of plant P and herbivore H biomass densities [g m⁻²] can be modelled as the following two coupled ordinary differential equations (Scheffer et al., 2008; O'Keeffe and Wieczorek, 2020):

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}=\mathit{\rho }\left(rt\right)P-C{P}^{\mathrm{2}}-Hg\left(P\right),\text{(A5)}& {\displaystyle \frac{\mathrm{d}H}{\mathrm{d}t}}=H\left(E{e}^{-bP}g\right(P)-m(rt\left)\right),\end{array}$$

with a non-monotone functional response,

$$\begin{array}{}\text{(A6)}& {\displaystyle }g\left(P\right)=c_{max}{\displaystyle \frac{P^{\mathrm{2}}}{P^{\mathrm{2}}+a^{\mathrm{2}}}}{e}^{-b_{c}P}.\end{array}$$

Following the approach of O'Keeffe and Wieczorek (2020), we fix six of the eight parameters (see Table 2.1, O'Keeffe and Wieczorek, 2020) and allow the plant growth rate ρ(rt) [1 d⁻¹] and herbivore mortality rate m(rt) [1 d⁻¹] to vary over time, subject to environmental conditions. In this system, the external forcing is given by

$$\begin{array}{}\text{(A7)}& \mathit{\lambda }\left(rt\right)=\left(\begin{array}{c}\mathit{\rho }\left(rt\right)\\ m\left(rt\right)\end{array}\right)=\left(\begin{array}{c}{\mathit{\rho }}_{-}+{\mathrm{\Delta }}_{\mathit{\rho }}a\left(rt\right)\\ m_{-}+{\mathrm{\Delta }}_{m}a\left(rt\right)\end{array}\right),\end{array}$$

where a(rt) takes the form of a ramp or a return profile, and ρ(rt) and m(rt) are varied proportionally to each other according to

$$\begin{array}{}\text{(A8)}& m\left(rt\right)=g\left(\mathit{\rho }\right(rt)-{\mathit{\rho }}_{-})+m_{-},\end{array}$$

where $g=\mathrm{1}/\mathrm{30}$. The forcing profile of the environmental conditions e(rt), corresponding to a normalised shift in ρ and m, is given by

$$\begin{array}{}\text{(A9)}& e\left(rt\right)={\mathrm{\Delta }}_{e}a\left(rt\right).\end{array}$$

The peak change in environmental conditions, used in Fig. 5a, is given by

$$\begin{array}{}\text{(A10)}& {\mathrm{\Delta }}_e={\displaystyle \frac{{\mathrm{\Delta }}_{\mathit{\rho }}}{{\mathrm{\Delta }}_{\mathit{\rho },\text{Fold}}}}={\displaystyle \frac{{\mathrm{\Delta }}_m}{{\mathrm{\Delta }}_{m,\text{Fold}}}},\end{array}$$

where ${\mathrm{\Delta }}_{\mathit{\rho },\text{Fold}}=\mathrm{0.1966}$, and ${\mathrm{\Delta }}_{m,\text{Fold}}=\mathrm{0.00655}$.

Figure 4a–c are obtained using the parameter values T=200, ${\mathit{\rho }}_{-}=\mathrm{0.5}$, $m_{-}=\mathrm{0.125}$, Δ_ρ=0.1, Δ_m=0.0033 (Δ_e=0.51), r=0.03 (green trajectory), r≈0.0453 (blue trajectory), and r=0.1 (red trajectory) and the initial condition at the location of the base state for $\mathit{\rho }={\mathit{\rho }}_{-}$ and $m=m_{-}$.

Figure 5a is obtained using a value of T where a(0) is no larger than 10⁻⁵ and the initial condition at the location of the base state for $\mathit{\rho }={\mathit{\rho }}_{-}$ and $m=m_{-}$ where it is equal to 1 at the fold bifurcation.

A4 AMOC model

Wood et al. (2019) proposed a five-box model to model the global oceanic current. The model is driven by salinity fluxes [psu] S_i in the main ocean waters: the North Atlantic (N), Tropical Atlantic (T), Indo-Pacific (IP), Southern Ocean (S), and Bottom waters (B). The flow strength q [Sv] of the Atlantic meridional overturning circulation (AMOC) is subsequently determined by a (fixed) temperature gradient ΔT and variable salinity gradient $\mathrm{\Delta }S=S_{\mathrm{N}}-S_{\mathrm{S}}$ between the North Atlantic and Southern Ocean boxes as follows:

$$\begin{array}{}\text{(A11)}& q={\displaystyle \frac{\mathit{\lambda }(\mathit{\alpha }\mathrm{\Delta }T+\mathit{\beta }\mathrm{\Delta }S)}{\mathrm{1}+\mathit{\lambda }\mathit{\alpha }\mathit{\mu }}}.\end{array}$$

Alkhayuon et al. (2019) empirically highlighted that the salinities of the Southern Ocean and Bottom waters vary much slower than the other boxes. Therefore, assuming these salinities are fixed and the salinity in the Indo-Pacific S_IP can be determined from a conservation of salinity, the original model reduces to the following two-dimensional model:

$$\begin{array}{}\text{(A12)}& {\displaystyle }{V}_{\mathrm{N}}{\displaystyle \frac{\mathrm{d}{S}_{\mathrm{N}}}{\mathrm{d}t}}=q(S_{\mathrm{T}}-S_{\mathrm{N}})+K_{\mathrm{N}}(S_{\mathrm{T}}-S_{\mathrm{N}})-F_{\mathrm{N}}\left(rt\right)S_{\mathrm{0}},\text{(A13)}& {\displaystyle }\begin{array}{rl}{V}_{\mathrm{T}}{\displaystyle \frac{\mathrm{d}{S}_{\mathrm{T}}}{\mathrm{d}t}}& =q[\mathit{\gamma }{S}_{\mathrm{S}}+(\mathrm{1}-\mathit{\gamma })S_{\mathrm{IP}}-S_{\mathrm{T}}]+K_{\mathrm{S}}(S_{\mathrm{S}}-S_{\mathrm{T}})\\ & +K_{\mathrm{N}}(S_{\mathrm{N}}-S_{\mathrm{T}})-F_{\mathrm{T}}\left(rt\right)S_{\mathrm{0}},\end{array}\end{array}$$

for q≥0 and

$$\begin{array}{}\text{(A14)}& {\displaystyle }{V}_{\mathrm{N}}{\displaystyle \frac{\mathrm{d}{S}_{\mathrm{N}}}{\mathrm{d}t}}=\left|q\right|(S_{\mathrm{B}}-S_{\mathrm{N}})+K_{\mathrm{N}}(S_{\mathrm{T}}-S_{\mathrm{N}})-F_{\mathrm{N}}\left(rt\right)S_{\mathrm{0}},\text{(A15)}& {\displaystyle }\begin{array}{rl}{V}_{\mathrm{T}}{\displaystyle \frac{\mathrm{d}{S}_{\mathrm{T}}}{\mathrm{d}t}}& =\left|q\right|(S_{\mathrm{N}}-S_{\mathrm{T}})+K_{\mathrm{S}}(S_{\mathrm{S}}-S_{\mathrm{T}})\\ & +K_{\mathrm{N}}(S_{\mathrm{N}}-S_{\mathrm{T}})-F_{\mathrm{T}}\left(rt\right)S_{\mathrm{0}},\end{array}\end{array}$$

for q<0. For more details about the parameters and their values, we refer the reader to Tables 3 and 4 in Alkhayuon et al. (2019).

In this study, we are interested in the surface freshwater fluxes F_i [Sv] and i in $\mathit{\{}N,T,\mathrm{IP},S\mathit{\}}$ as input parameters for the system. Following the approach of Alkhayuon et al. (2019) and Wood et al. (2019), these fluxes are defined as linear functions of a hosing parameter H [Sv], such that the total flux is 0 for all H and H=0 corresponds to the baseline values of F_i (Table 3, Alkhayuon et al., 2019). In this system, the external forcing is given by

$$\begin{array}{}\text{(A16)}& \mathit{\lambda }\left(rt\right)=\left(\begin{array}{c}{F}_{\mathrm{N}}\left(rt\right)\\ F_{\mathrm{T}}\left(rt\right)\\ F_{\mathrm{IP}}\left(rt\right)\\ F_{\mathrm{S}}\left(rt\right)\end{array}\right)=\left(\begin{array}{c}\mathrm{0.486}+\mathrm{0.1311}H\left(rt\right)\\ -\mathrm{0.997}+\mathrm{0.6961}H\left(rt\right)\\ -\mathrm{0.754}-\mathrm{0.5646}H\left(rt\right)\\ \mathrm{1.265}-\mathrm{0.2626}H\left(rt\right)\end{array}\right),\end{array}$$

where

$$H\left(rt\right)={\mathrm{\Delta }}_{H}a\left(rt\right),$$

and a(rt) takes the form of a ramp or a return profile. Figure 4d–f are obtained using the parameter values T=1000, Δ_H=0.365, r=0.005 (green trajectory), r≈0.01645 (blue trajectory), and r=0.017 (red trajectory) and the initial condition at the location of the base state (the AMOC-On state) for H=0. Figure 5b was obtained using a value of T where a(0) is no larger than 10⁻⁵ and the initial condition at the location of the base state for H=0.

A5 The power grid model

We use a three-bus power system model from Dobson et al. (1988) and Dobson and Chiang (1989) to represent the dynamics that can be observed on a power grid. Two generators supply power to a P–Q load in parallel with a capacitor and induction motor (Revel et al., 2006). The model consists of the following four differential equations for the generator phase angle δ_m, the angular velocity ω_m, the phase angle δ, and the magnitude V of the load voltage (Ajjarapu and Lee, 1992):

$$\begin{array}{}\text{(A17)}& \begin{array}{rl}{\dot{\mathit{\delta }}}_{\mathrm{m}}& ={\mathit{\omega }}_{\mathrm{m}},\\ {\dot{\mathit{\omega }}}_{\mathrm{m}}& ={\displaystyle \frac{\mathrm{1}}{M}}(-d_{\mathrm{m}}{\mathit{\omega }}_{\mathrm{m}}+P_{\mathrm{m}}-P_{\mathrm{e}}),\\ \dot{\mathit{\delta }}& ={\displaystyle \frac{\mathrm{1}}{K_{q\mathit{\omega }}}}\left(-K_{qv}V-K_{q{v}^{\mathrm{2}}}{V}^{\mathrm{2}}+Q_{\mathrm{l}}-Q_{\mathrm{0}}-Q_{\mathrm{1}}\left(rt\right)\right),\\ \dot{V}& ={\displaystyle \frac{\mathrm{1}}{T{K}_{q\mathit{\omega }}{K}_{pv}}}(K_{p\mathit{\omega }}{K}_{q{v}^{\mathrm{2}}}{V}^{\mathrm{2}}+(K_{p\mathit{\omega }}{K}_{qv}-K_{q\mathit{\omega }}{K}_{pv})V\\ & +K_{p\mathit{\omega }}(Q_{\mathrm{0}}+Q_{\mathrm{1}}(rt)-Q_{\mathrm{l}})-K_{q\mathit{\omega }}(P_{\mathrm{0}}+P_{\mathrm{1}}-P_{\mathrm{l}})).\end{array}\end{array}$$

where P_e,P_l, and Q_l are given by

$$\begin{array}{}\text{(A18)}& \begin{array}{rl}{P}_{\mathrm{l}}& =-V_{\mathrm{0}}V{T}_{\mathrm{0}}\sin (\mathit{\delta }+{\mathit{\theta }}_{\mathrm{0}})-V_{\mathrm{m}}V{Y}_{\mathrm{m}}\sin (\mathit{\delta }-{\mathit{\delta }}_{\mathrm{m}}+{\mathit{\theta }}_{\mathrm{m}})\\ & +\left(Y_{\mathrm{0}}\sin \right({\mathit{\theta }}_{\mathrm{0}})+Y_{\mathrm{m}}sin({\mathit{\theta }}_{\mathrm{m}}\left)\right)V^{\mathrm{2}},\\ Q_{\mathrm{l}}& =V_{\mathrm{0}}V{Y}_{\mathrm{0}}\cos (\mathit{\delta }+{\mathit{\theta }}_{\mathrm{0}})+V_{\mathrm{m}}V{Y}_{\mathrm{m}}\cos (\mathit{\delta }-{\mathit{\delta }}_{\mathrm{m}}+{\mathit{\theta }}_{\mathrm{m}})\\ & -\left(Y_{\mathrm{0}}\cos \right({\mathit{\theta }}_{\mathrm{0}})+Y_{\mathrm{m}}\cos ({\mathit{\theta }}_{\mathrm{m}}\left)\right)V^{\mathrm{2}}),\\ P_{\mathrm{e}}& =-V_{\mathrm{m}}V{Y}_{\mathrm{m}}\sin (\mathit{\delta }-{\mathit{\delta }}_{\mathrm{m}}{\mathit{\theta }}_{\mathrm{m}})-V_{\mathrm{m}}^{\mathrm{2}}{Y}_{\mathrm{m}}\sin \left({\mathit{\theta }}_{\mathrm{m}}\right).\end{array}\end{array}$$

In this system, the external forcing is in the form of the reactive power demand (referred to as power demand in the main text) of the load Q₁(rt) and is given by

$$\begin{array}{}\text{(A19)}& {\displaystyle }\mathit{\lambda }\left(rt\right)=Q_{\mathrm{1}}\left(rt\right)=Q_{{\mathrm{1}}_{-}}+{\mathrm{\Delta }}_{Q_{\mathrm{1}}}a\left(rt\right),\end{array}$$

where a(rt) takes the form of a ramp or a return profile. Figure 6a and b are obtained using the parameter values T=0.15, $Q_{{\mathrm{1}}_{-}}=\mathrm{5}$, ΔQ=5, r=50 (blue trajectory), r=130 (orange trajectory), and r=400 (purple trajectory) and the initial condition at the location of the base state at $Q_{\mathrm{1}}=Q_{{\mathrm{1}}_{-}}$. Figure 6c is obtained using a value of T where a(0) is no larger than 10⁻², the parameter value $Q_{{\mathrm{1}}_{-}}=\mathrm{5}$, and the initial condition at the location of the base state at $Q_{\mathrm{1}}=Q_{{\mathrm{1}}_{-}}$

For a full description of the other parameters and the values used in this study, we refer the reader to Budd and Wilson (2002). However, we choose to set P₁=5, the real power demand of the load, such that when Q₁(rt) is varied according to Eq. (A19), the system can only cross a fold bifurcation and does not encounter a Hopf bifurcation, as observed for P₁=0 (Wang et al., 1994).