# Stabilising millennial oscillations in large-scale ocean circulation with a delayed feedback due to a circumpolar current

Andrew Keane, Alexandre Pohl, Henk A Dijkstra, Andy Ridgwell

► To cite this version:

Andrew Keane, Alexandre Pohl, Henk A Dijkstra, Andy Ridgwell. Stabilising millennial oscillations in large-scale ocean circulation with a delayed feedback due to a circumpolar current. Physica D: Nonlinear Phenomena, 2025, 476, pp.134680. 10.1016/j.physd.2025.134680. hal-05048600

HAL Id: hal-05048600

https://hal.science/hal-05048600v1

Submitted on 28 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Stabilising millennial oscillations in large-scale ocean circulation with a delayed feedback due to a circumpolar current

Andrew Keane<sup>a,b,cID,*</sup>, Alexandre Pohl<sup>cID</sup>, Henk A. Dijkstra<sup>d,eID</sup>, Andy Ridgwell<sup>fID</sup>

## ARTICLE INFO

Communicated by Stefan Wiggins

Dataset link: 10.5281/zenodo.14606515, 10.52 81/zenodo.1407657, https://github.com/andr ewkeane/kpdr_oscillations

**Keywords:**

Delay dynamics

Bifurcation analysis

Ocean circulation

Oscillation

## ABSTRACT

The global ocean circulation plays a pivotal role in the regulation of the Earth's climate. The specific pattern and strength of circulation also determines how carbon and nutrients are cycled and via the resulting distribution of dissolved oxygen, where habitats suitable for marine animals occur. However, evidence from both geological data and models suggests that state transitions in circulation patterns have occurred in the past. Understanding the controls on marine environmental conditions and biodiversity requires a full appreciation of the nature and drivers of such transitions. Here we present stable millennial oscillations of meridional overturning circulation in an Earth system model of intermediate complexity, cGENIE, that appear to only occur in the presence of a circumpolar current. To demonstrate that a circumpolar current can act as a driver of stable oscillations, we adapt a simple ocean box model to include a delayed feedback to represent the effect of a circumpolar current on meridional overturning circulation. We investigate the millennial oscillatory solutions that arise in the box model by bifurcation analysis and show that the model can reproduce the same bifurcation structure observed in the Earth system model. Our results provide new insights into the nature of oscillations that could have occurred under certain continental configurations in the geological past, and also highlight the potential influence of the changing Antarctic circumpolar current speed on the stability of the Atlantic meridional overturning circulation.

## 1. Introduction

The ocean circulation plays a critical role in the Earth's global climate system. Wind-driven currents acting over the first hundreds of metres impose the first-order pattern and magnitude of ocean heat transport, and therefore largely determine superficial climate and notably the latitudinal gradient of temperature [1]. Upper-ocean currents also control equatorial and coastal upwelling position and strength, thus the spatial patterns of primary productivity and extension of oxygen minimum zones [2]. The deep ocean circulation, on the other hand, is affected by density differences and is important for the redistribution of nutrients in the global ocean (and critically, their return back to the ocean surface) [3] as well as the sequestration of atmospheric carbon in the ocean interior [4]. Overturning circulation also controls the transport of oxygen equilibrated in surface seawater into the ocean interior ('ventilation') and, ultimately, the capacity of the ocean to host animal life at depth [5,6].

Abundant sedimentological, geochemical, and paleontological data demonstrate that the global ocean circulation has significantly varied through geological time. At the (multi-)million-year scale, it has varied in response to continental rearrangement. For instance, Trabucho Alexandre et al. [7] propose that intense deoxygenation and the deposition of organic-rich sediments developed in the mid-Cretaceous North Atlantic, around 100 million years ago, as a consequence of a specific ocean circulation pattern making this ocean basin act as a nutrient trap. A similar mechanism was suggested by Meyer et al. [8] to explain photic-zone euxinia expansion in the Paleo-Tethys ocean at the Permian-Triassic boundary. At a higher frequency, the ocean circulation, both shallow and deep, is also known to vary in concert with global climate. Today, vigorous deep-water formation takes place in the North Atlantic (and Southern) Ocean while limited ventilation is observed at similar northern latitudes in the Pacific Ocean [9]. During the Last Glacial Maximum, sediment core proxy data for ocean

<footnote><sup>a</sup> School of Mathematical Sciences, University College Cork, Cork, Ireland</footnote>

<footnote><sup>b</sup> Environmental Research Institute, University College Cork, Cork, Ireland</footnote>

<footnote><sup>c</sup> Biogéosciences UMR 6282, Université Bourgogne Europe, CNRS, F-21000, Dijon, France</footnote>

<footnote><sup>d</sup> Institute for Marine and Atmospheric Research Utrecht, Department of Physics, Utrecht University, Utrecht, The Netherlands</footnote>

<footnote><sup>e</sup> Center for Complex Systems Studies, Department of Physics, Utrecht University, Utrecht, The Netherlands</footnote>

<footnote><sup>f</sup> Department of Earth and Planetary Sciences, University of California, Riverside, CA, USA</footnote>

<footnote>* Corresponding author.
E-mail address: andrew.keane@ucc.ie (A. Keane).</footnote>
ventilation suggest that the North Pacific Ocean was better ventilated at intermediate depths [9], while the Atlantic circulation was probably more sluggish (shallower and/or less vigorous) [10]. Similarly, Freeman et al. [9] used radiocarbon measurements to reconstruct the benthic–planktonic $^{14}$C age offset (i.e. ‘ventilation age’) during the last deglaciation and demonstrated a see-saw in the ventilation of the intermediate Atlantic and Pacific Oceans.

In order to assimilate the proxy database in a quantitative framework and gain mechanistic understanding of the drivers of observed changes in ocean properties, (paleo)climatologists use ocean-atmosphere general circulation models (GCMs) [11–13]. These models provide complex representations of ocean physics and dynamics and the coupling with the atmospheric component, sea ice, and possibly ocean biogeochemistry and land surface [14]. Models of the Coupled Model Intercomparison Project phase 6 (CMIP6), for instance, aim to offer a satisfactory representation of the climate system [15]. However, these models are computationally expensive and, consequently, are generally run for only a few hundred to a thousand years. Because this duration approaches the ocean mixing timescale, it is assumed that it allows the model to approximately reach deep-ocean equilibrium, i.e., the physical characteristics of the whole ocean (temperature and salinity), and thus ocean dynamics, are supposed stable from that point in time [11,16]. By doing so, it is also implicitly assumed that the ocean circulation does in fact eventually reach a stable steady state. Yet, various attempts to run global climate models for a longer time period demonstrated more complex behaviours of the ocean circulation, which sometimes exhibit regimes of stable self-sustained oscillations. Oscillations have notably been simulated in ocean-only GCMs [17,18] and ocean GCMs coupled to simple atmospheric components (Earth system models of intermediate complexity; EMICs) [19,20]. To our knowledge, the only studies reporting such self-sustained oscillations based on ocean-atmosphere GCM simulations are Refs. [21,22] using the fully coupled HadCM3 GCM, Ref. [23] using the MIROC4m GCM, Ref. [24] using the Max Planck Institute Earth system model and the last glacial CESM simulations of Vettoretti et al. [25]. The absence of further documentation of stable oscillatory regimes in other fully coupled GCMs is not surprising considering that a simulation duration of at least 6 thousand years seems to be required to identify stable oscillations, while a typical run-time for GCMs used to explore past (geological) climate states is typically only about 3 thousand years (e.g. [13,26,27]).

The existence of regimes of stable self-sustained oscillations in ocean circulation has recently been demonstrated by Pohl et al. [3] for specific combinations of past continental configuration and climatic state using the EMIC cGENIE. The oscillations have been invoked to explain both the poorly-oxygenated oceans that pervaded several hundreds of millions of years ago [28], when marine life started diversifying in the marine realm [29], as well as enigmatic thousand-year scale cycles identified in marine sediments of the deep geological past [30]. Based on recent studies that pinpoint the tight coupling between marine redox conditions and marine biodiversity [6,31–33], the authors conclude that the simulated high-frequency variability in ocean redox conditions would have constituted harsh conditions for the development of marine life at that time. The study by Pohl et al. [3] exemplifies the implications that such oscillations have for our understanding of the evolution of both global climate and marine biodiversity, but also highlights the need to better characterise the physical mechanisms causing these oscillations. It focused on the consequences of the resulting state transitions in ocean circulation for ocean oxygenation and the interpretation of geological redox indicators, rather than on the underlying mechanisms of the oscillations. In addition, the experimental setup used by the authors, involving complex land mass shapes and configurations, as well as the presence in the model of a variety of climatic and biogeochemical feedbacks, creates challenges for a detailed analysis. As a result, the physical mechanisms causing the stable self-sustained

oscillations and state transitions in ocean dynamics in cGENIE remain to be fully elucidated.

Here, we analyse the mechanisms at play by combining a series of new Earth system model experiments conducted using an idealised ocean configuration, with a box model designed specifically to reproduce and identify key characteristics of the cGENIE simulations. Experiments using cGENIE with the idealised continental configuration serve as a stepping-stone between results of oscillations in simulations using “real-world” continental configurations [3] and understanding gained through the analysis of simple ocean box models.

Box models have been often used to identify physical mechanisms of oscillatory behaviour in the ocean-climate system. Regarding variability of the meridional overturning circulation, two basic types of oscillatory behaviour have been identified. Overturning (or loop) oscillations can be found in a four-box model [34] and involve the advective propagation of a buoyancy anomaly along the overturning loop. Such variability is of centennial time scale when only the Atlantic is considered and millennial scale when the global overturning is considered [18]. Oscillatory variability may also be induced by convective processes. The basics of the behaviour can be understood from the so-called ‘flip-flop’ oscillation [35], where oscillations between convective and non-convective states occur. The origin of the millennial time-scale oscillations, also referred to as ‘deep-decoupling’ oscillations or ‘flushes’ is, however, in more detailed models best explained by a box model which combines both advective and convective processes [36].

Here, we present a model that reproduces the oscillations obtained in the cGENIE Earth system model [3] in a very simple way, by embedding key climatic features in a simple 3-box model. The novel element is that a time delay is introduced, which serves as an efficient representation of the effects of a circumpolar current on the oceanic salt transport.

## 2. EMIC experiments

To gain insight into the existence of self-sustained oscillations in the ocean circulation, numerical simulations are conducted using an idealised setup of the EMIC cGENIE [37]. cGENIE consists of a reduced physics (frictional geostrophic) 3D ocean circulation model coupled to a 2D energy-moisture-balance atmospheric component plus sea-ice module. Although we do not show the results of it here, cGENIE also includes a state-of-the-art representation of marine biogeochemistry. Compared to CMIP6 Earth system models, the simplified climatic component of cGENIE offers a rapid model integration time, permitting the model to be run for tens of thousands of years.

The adopted idealised cGENIE setup is “Drake World”. This differs slightly from the classic MIT GCM implementation [38], but is primarily characterised by a thin strip of land extending southwards to 45°S and a 5-km deep flat-bottomed ocean everywhere else (Fig. 1(a)). The model uses a 36 × 36 equal-area horizontal grid, with 16 unevenly-spaced z-coordinate vertical levels in the ocean. We apply idealised boundary conditions of zonally averaged wind stress and speed, plus a zonally averaged planetary albedo, following Vervoort et al. [39] and van de Velde et al. [40] (Fig. 1(b)). The solar constant is set to the modern value (1368 W m⁻²). The physical parameters controlling the model climatology follow Cao et al. [41]. Additional experiments are conducted using an alternative, idealised continental configuration with a strip of land extending from the North Pole to the South Pole (“Ridge World” [3,38]) and everything else unchanged.

Fig. 1(c) shows the barotropic stream-function (units of Sverdrup, 1 Sv = $10^6$ m³ s⁻¹), as a measure of the vertically averaged flow, in Drake World with atmospheric CO₂ at 16 times the pre-industrial atmospheric level (PAL, 1 PAL = 280 ppm). As with all the following results, this shows the state of the model following spin-up (i.e. after all transient behaviour has died away). The most notable feature of the barotropic stream-function is the strong zonal flow in the southern channel, similar to the Antarctic Circumpolar Current observed in the
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** cGENIE Drake World setup. (a) Model bathymetry. Landmasses (above sea level; a.s.l.) are shaded white. (b) Zonal wind stress field used as boundary condition. (c) Barotropic streamfunction (Sverdrup, $1 \text{ Sv} = 10^6 \text{ m}^{-3} \text{ s}^{-1}$) simulated at 16 times the pre-industrial atmospheric CO$_2$ concentration of 280 ppm. The boundary is periodic in the horizontal direction.</figcaption></figure>

Earth's modern continental configuration. It is a mainly wind-driven current that constantly transports water around the pole and persists also for the other levels of atmospheric CO$_2$ considered here.

We find that when atmospheric CO$_2$ forcing is varied in the model (from 4 to 16 PAL), ocean dynamics exhibits a complex response. Overall, we observe a gradual shift in the locus of deep-water formation from the (oceanic) North Pole (Fig. 2(a)) to the South Pole (Fig. 2(b)). Interestingly, during this shift we encounter two families of stable oscillations.

The blue points in Fig. 2(c) show mean salinity for an ensemble of simulations in Drake World. Throughout this paper, we take the mean salinity across the southern channel, comprising all depths, longitudes and 45–90°S. Because the simulations take a long time to run (around 6 days for a 60-kyr long cGENIE simulation as shown in Fig. 2(c)), they are run in parallel for different atmospheric CO$_2$ levels using the same initial conditions. As CO$_2$ increases, the ocean transitions from a stable steady state with northern deep-water formation to a regime of stable oscillations between northern deep-water formation and a weakened circulation state (represented by the envelopes in Fig. 2(c)). This weakened state becomes stable near 8.5 PAL. Further increasing of CO$_2$ gives rise to another regime of stable oscillations between a largely collapsed circulation state and strong circulation with southern deep-water formation. Finally, the southern deep-water formation steady state becomes stable for higher levels of atmospheric CO$_2$. Example time series of the oscillations are shown in Fig. 2(d). Although results are only shown here for salinity across the southern channel, these changes in ocean circulation have been shown to induce significant and global periodic changes in many other observables, such as deep-ocean ventilation and oxygenation (see Ref. [3], extended data Fig. 10).

An additional series of experiments was conducted on Ridge World, i.e., under identical boundary conditions but with the strip of land extending to the southernmost grid points. The results reveal a transition from a stable steady state with northern deep-water formation

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** cGENIE results. Meridional overturning streamfunction for Drake World at (a) 4 PAL and (b) 16 PAL. (c) Envelope of mean salinity (across the southern channel: comprising all depths, longitudes and 45–90°S) simulated at various atmospheric CO$_2$ levels, illustrating regimes of stable equilibria (where no envelope is visible) and stable oscillations (where the envelope can be seen), for Drake World (blue) and Ridge World (black). (d) Example time series of stable oscillations (legend indicates PAL) simulated using Drake World with transients (initial spin-up of 30 kyr) removed.</figcaption></figure>

to one with southern deep-water formation. This transition is sudden, represented by the jump of black points in Fig. 2(c) near 4 PAL, and is akin to the classic “tipping point” behaviour of the Atlantic meridional overturning circulation observed in many studies [42–44]. However, importantly, this transition does not exhibit regimes of stable oscillations, at least not for the levels of CO$_2$ tested. A primary question addressed in this paper is then: Why do these oscillations appear once the ocean has an open-ocean passage around the South Pole?

### 3. Simple box model

We begin with the simplified ocean box model from section 3 of Ref. [45], which consists of two polar boxes and one upper equatorial box, as illustrated below in Fig. 3. Assuming salt conservation, the model is written in terms of the salt content of the southern polar box (box 1) and the northern polar box (box 2). The model, as it appears in Ref. [45], represents only meridional transport and its behaviour corresponds to the meridional overturning circulation of Ridge World. In Drake World, on the other hand, the presence of a circumpolar current leads to the zonal transport of salt in the southern channel.

Here, we adapt the ocean box model to capture the effect of this zonal transport by including a delayed feedback term on the salinity of
<figure><img src="image_5.png" /><figcaption>Fig. 3. Boxes representing regions of the ocean with salinity variables $S_{i \in \{1,2,3\}}$ with arrows representing meridional water flux of strength $m$ and zonal flux of the southern channel of strength $\sigma$. Parameters $F_1$ and $F_2$ represent freshwater fluxes.</figcaption></figure>

box 1. In particular, a parcel of water in the southern polar box will leave its position at time $t$, travel around the southern pole, and return to the same position at time $t + \tau$. This creates a delayed effect on the average salinity in box 1. The model becomes a set of delay differential equations (DDEs) for the salinity variables $S_{i \in \{1,2\}}$:

$$V \dot{S}_1 = S_0 F_1 + m (S_2 - S_1) + \sigma (S_1^\tau - S_1) \quad (1)$$

$$V \dot{S}_2 = -S_0 F_2 + m (S_3 - S_2).$$

Parameter $V$ is the (equal) volume of the boxes and $S_0$ is a reference salinity. Parameters $F_1$ and $F_2$ represent surface freshwater fluxes between the boxes due to atmospheric water vapour transport and wind-driven surface currents. The delayed variable $S_1^\tau = S_1(t - \tau)$ has delay time $\tau$ and $\sigma$ is the feedback strength. In other words, we have water with salt density $S_1$ leaving and salt density $S_1^\tau$ entering box 1 with $\sigma$ representing zonal flux. The variable $m$ represents the volume transport between boxes due to the density gradient across the ocean:

$$m = k [\beta (S_2 - S_1) - \alpha T^*]. \quad (2)$$

Parameter $k$ is a hydraulic constant, while $\beta$ and $\alpha$ are salinity and temperature expansion coefficients, respectively. The model assumes a fixed temperature gradient $T^*$ between the northern and southern boxes. In Fig. 3, a positive $m$ corresponds to clockwise circulation. Note that here we only discuss the positive $m$ solutions, and therefore consign the negative $m$ equations to Appendix A for completeness. Assuming that salt is conserved, we have $S_3 = 3S_0 - S_1 - S_2$; see Appendix B for a discussion of the validity of this salt conservation rule.

In this paper we take parameter values from Ref. [45] and let $k = 23 \times 10^{17} \text{ m}^3 \text{ yr}^{-1}$, $\alpha = 1.7 \times 10^{-4} \text{ K}^{-1}$, $\beta = 0.8 \times 10^{-3} \text{ psu}^{-1}$, and $S_0 = 35 \text{ psu}$. We let $V = 3.5 \times 10^{17} \text{ m}^3$, which is approximately the volume of the southern channel in Drake World. Furthermore, we set $F_2 = 1 \text{ Sv} (10^6 \text{ m}^3 \text{ s}^{-1})$ and $T^* = 0 \text{ K}$, and note that qualitatively similar results to those presented below are found when varying these parameters. Finally, we approximate delay times for each latitude within the southern channel using averaged zonal ocean speeds from cGENIE simulations across longitude and depth. The average of these delay times gives us a rough estimate for $\tau$ of about 1000 years. Therefore, we use similar values for $\tau$ in the results presented below.

The presence of a circumpolar current creates a delayed effect on the average salinity in the southern polar region. We acknowledge that a delayed feedback is not the only way to model the effect of the circumpolar current. For example, one could instead use multiple boxes to represent zonal flow across the southern channel. We suggest that a delayed feedback is a relatively elegant modelling approach, since it is mathematically concise, requiring only one additional term to the model, and allows for some simple physical interpretation. In the case of zero wind speed across the southern channel, there is no significant zonal flux, such that $\sigma$ approaches zero and the delayed feedback effect becomes negligible. In the case of $\sigma$ equal to zero, the model reverts to

representing only meridional flow and corresponds to the behaviour of Ridge World rather than Drake World. For non-zero wind speed across the southern channel, the wind speed drives the zonal flow and is, therefore, inversely proportional to the delay time $\tau$. The delayed feedback term does not change the steady-state solutions of the system, since $S_1^\tau = S_1$ for all steady states, such that the term vanishes. It may, nonetheless, affect the stability of a steady state. When $\tau$ is sufficiently small the delayed feedback may have a stabilising effect. For example, if $S_1(t)$ experiences a positive perturbation away from its steady-state value, then $S_1 > S_1^\tau$ such that the delayed feedback works to decrease $S_1$, since $S_1 \propto (S_1^\tau - S_1)$. However, if $\tau$ and $\sigma$ are sufficiently large, this restoring behaviour of the delayed feedback may overshoot the steady state and instead promote oscillatory dynamics. This type of phenomenon is typical of delay-induced dynamics; see, for example, chapter 2 of Ref. [46]. Of course, we could use additional delayed feedback terms to represent the influence of further periodic flows, such as the surface gyres seen in Fig. 1(c). However, these surface effects would translate to delayed feedback terms with very small $\sigma$ and $\tau$ values, resulting in negligible effects for the meridional overturning, and are, therefore, neglected. It is only the periodic flow in the southern channel that, due to the lack of land mass, possesses enough momentum to generate a deep flow with sufficiently long delay times, resulting in large $\tau$ values.

Clearly, the box model representation is very much highly simplified compared to a dynamic 3D ocean circulation model. As such, we do not attempt to find quantitative agreement between the two models. Rather, we demonstrate that this simple model possesses surprisingly rich dynamics and is sufficient for a qualitative reproduction of the interesting ocean circulation behaviour observed in cGENIE experiments.

In Ref. [45] the model without the delayed feedback was shown to exhibit, not only the typical bistability between upper ($m > 0$) and lower ($m < 0$) branches of solutions, but also a Bogdanov-Takens bifurcation. This means that instead of the upper branch losing stability at a fold bifurcation, it may for certain values of freshwater flux already lose stability at a Hopf bifurcation, which is always subcritical [47]. The resulting unstable periodic orbits terminate at homoclinic bifurcations. We now investigate the effect that the delayed feedback has on this bifurcation structure. As shown in Appendix C, the analysis of model (1) is complicated by the delay term, such that we primarily rely on numerical methods. In the following section, solutions to the model and their stability properties are found using the continuation package DDE-Biftool [48–51].

## 4. Stable oscillations

The solutions to the box model shown in Fig. 4 demonstrate how the delayed feedback effect allows the existence of stable periodic orbits. Black curves show the steady-state solutions of $S_1$ as parameter $F_1$ is varied, while the blue curves indicate the maximum of $S_1$ for periodic solutions. Solid and dashed curves indicate stable and unstable solutions, respectively. In panel (a) the parameter $\sigma$ is set to zero; that is, there is no delayed feedback effect. In this case, as $F_1$ increases, the stable steady-state solution loses stability at a subcritical Hopf bifurcation, resulting in a branch of unstable periodic orbits. When the delayed feedback is active, the Hopf bifurcation may become supercritical, resulting in a branch of stable periodic orbits, as shown in panel (b). In this particular example, the stable periodic solutions lose stability at a fold of periodic orbits. The periods of the resulting stable oscillations have a millennial timescale as demonstrated by the example time series shown in panel (c).

Some readers familiar with control theory may not be surprised by the ability of the delayed feedback to stabilise the pre-existing periodic orbits. Curiously, the additional term in our model that represents the zonal flow takes the same mathematical form as Pyragas control [52]. This control method is used to stabilise unstable periodic orbits, as has
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Solutions of model (1) depending on freshwater forcing parameter $F_1$. Solid/dashed curves represent stable/unstable solutions. Black are equilibria and blue are maxima of periodic solutions. White and black filled circles represent Hopf and fold bifurcations, respectively. (a) is without the delayed feedback (i.e. $\sigma = 0$), while (b) is with $\sigma = 11$ Sv. The delay time $\tau$ is set to 900 years. (c) shows an example time series with $\sigma = 11$ Sv and $F_1 = -0.208$ Sv.</figcaption></figure>

been demonstrated in many applications, for example, in semiconductor laser experiments [53]. To be precise, the delayed feedback effect in our model is not equivalent to Pyragas control. In the latter, the control is applied to all variables of the system (we only apply it to the salinity in one box) and the delay time is set equal to the period of the unstable periodic orbit to be stabilised. Nonetheless, the method may still be successful for other delay times, as investigated in Ref. [54]. In model (1) we find that the periods of the stable oscillations are generally at least twice the delay time.

Whether the delayed feedback induces a change in criticality of the Hopf bifurcation depends on the feedback parameters $\sigma$ and $\tau$. Fig. 5 illustrates this dependence, where curves of Hopf bifurcations are calculated in the $(F_1, \sigma)$-plane for different values of $\tau$. The thick parts of the curves indicate that the Hopf bifurcation is supercritical. Notice that, depending on the feedback strength and delay time, it is possible that the system experiences multiple Hopf bifurcations as $F_1$ increases.

Fig. 6(a) is a bifurcation diagram, or more precisely an “attractor diagram”, showing the amplitude of oscillations in terms of mean salinity (across the southern channel of Drake World) constructed from simulation results of cGENIE. We observe that as $CO_2$ is increased the steady-state solution loses stability near 5 PAL. The amplitude

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Curves of Hopf bifurcations in the $(F_1, \sigma)$-plane. At a Hopf bifurcation the steady-state solution loses (or gains) stability and produces a branch of periodic solutions. Where the curves are thin/thick the Hopf bifurcation produces a branch of unstable/stable periodic solutions. The colours refer to different delay times $\tau$ shown in the legend.</figcaption></figure>

growth near this bifurcation point resembles a square-root law, typical of a Hopf bifurcation [55]. The family of periodic solutions disappears at, what appears to be, another supercritical Hopf bifurcation near 8 PAL. The bifurcation point must be very close to 8 PAL, since it was necessary to run the 8 PAL simulation for an extremely long time (up to 120 thousand years) while the amplitude of the oscillations continued to decrease very slowly. Finally, the steady-state solution again loses stability after 8.5 PAL.

As illustrated in Fig. 5 we may choose parameter values so that the box model will pass through two supercritical Hopf bifurcations as $F_1$ increases. Realistically, $F_1$ is not the only parameter value that would change due to changes in atmospheric $CO_2$. Parameters $F_2, T^*, \sigma$ and $\tau$ would all be susceptible to influence from atmospheric changes. Nonetheless, we only consider varying $F_1$ in an effort to keep the analysis simple and demonstrate qualitative changes in behaviour. Fig. 6(b) shows that the effect of the circumpolar current is sufficient for reproducing a family of periodic solutions between two supercritical Hopf bifurcations. In the box model the steady-state solution finally loses stability at a subcritical Hopf bifurcation.

Fig. 7(a) is the same as Fig. 6(a), only for higher atmospheric $CO_2$. Near 8.5 PAL we observe a jump to large-amplitude oscillations, as well as bistability between steady-state and oscillatory solutions. Both the jump to large-amplitude oscillations and bistability are reproduced by the box model in panel (b). Here, the steady-state solution loses stability in a subcritical Hopf bifurcation. The unstable periodic orbits become stable via a fold bifurcation of periodic orbits, before later losing stability at a period-doubling bifurcation.

The termination of the large-amplitude oscillations in Fig. 7(a) involves co-existing oscillatory solutions at 10 PAL and a gradual transition to a state of strong convection at the southern pole. These solutions are shown as filled blue circles in panel (a). Following these oscillatory solutions, for larger values of $CO_2$, we only find steady-state solutions that correspond to $m < 0$ solutions in the box model, which we do not consider here for simplicity. The re-creation of transitions between strong and weak convecting states in the box model would require upper and lower boxes, as well as the introduction of an additional convection timescale into the dynamics as in Ref. [36], and is therefore beyond the scope of model (1).
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Amplitudes of oscillations with steady-state and periodic solutions shown in black and blue, respectively. (a) cGENIE experiments with varying atmospheric CO₂. (b) Model (1) solutions while varying *F*₁ with *σ* = 9.5 Sv and *τ* = 1100 years. Solid/dashed curves represent stable/unstable solutions. White filled circles represent Hopf bifurcations.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Amplitudes of oscillations with steady-state and periodic solutions shown in black and blue, respectively. (a) cGENIE experiments with varying atmospheric CO₂. Filled blue circles indicate solutions not captured by the box model. (b) Model (1) solutions while varying *F*₁ with *σ* = 9 Sv and *τ* = 850 years. Solid/dashed curves represent stable/unstable solutions. White, black and grey filled circles represent Hopf, fold and period-doubling bifurcations, respectively.</figcaption></figure>

## 5. Discussion

The delayed feedback provides the simple box model with the capacity to mimic important qualitative behaviour observed in cGENIE experiments using an idealised configuration. It hence provides a simple mechanism for the occurrence of stable oscillations in the ocean circulation and its dependence on a zonal flow across the southern channel. The analysis, conducted using the numerical continuation software DDE-Biftool, highlights the need for each system parameter to be within a certain range of values in order for oscillations to be observed. For the cGENIE model, this implies that oscillations could easily be missed if, for example, the wind speed over the southern channel is too large (see Appendix D for a discussion on the role of wind speed).

The results presented here provide a possible explanation for observations in cGENIE experiments using realistic continental configurations of the geological past, where stable oscillations coincide with configurations having a large open ocean pole [3]. These cGENIE oscillations, which occur on the millennial timescale, could shed light on periods of high variability observed in high-resolution proxy data [30].

Despite some evidence of stable oscillatory behaviour in GCMs [22], the extended run-time of such models makes them ill-suited for the explorations through parameter space required when investigating possible causes and attributing factors of oscillatory behaviour. cGENIE, on the other hand, represents a very practical compromise for establishing connections between interesting dynamics and fundamental mechanisms that drive them. As such, this study is conducted in a similar vein to some previous studies that have conducted parameter sweeps in order to identify possible bifurcations and regions of qualitatively different climate states (for example, see Refs. [56–58]).

The way that the stable oscillations presented here are generated by delayed feedback (in turn, by continental configuration) is different to previous studies of ocean models. In some past studies convection and associated deep-water flushes have been the driving mechanism for oscillations [35,36]. In other studies sea-ice dynamics are required as part of a feedback mechanism that generates oscillations [59–61]. One logical next step is to investigate how the delayed feedback interacts with these other sources of oscillations. There is evidence from cGENIE simulations (not shown) of quasi-periodic behaviour, which could be the result of an interaction between multiple mechanisms for oscillations.

The simple box model with delayed feedback does not reproduce all dynamical features observed in the cGENIE experiments, such as the oscillations shown as blue filled circles in Fig. 7(a). Furthermore, although the simple model captures the millennial timescale of the stable oscillations, the period of the oscillations observed in cGENIE can easily exceed 10 thousand years (for example, the large-amplitude oscillations around 9 PAL in Fig. 2(c)), which is not found in the simple model. A further modelling possibility that would effect the timescales in the box model includes adding more boxes to the box model. Finally, deep-water flushes could be captured by a salinity delay due to a deep-ocean adjustment to convective events, which would further increase the delay time of the model.

Results from our box model may also have implications for our understanding of future changes in the global ocean circulation. It has been reported that part of the Antarctic Circumpolar Current (ACC) is accelerating due to anthropogenic ocean warming [62]. This would correspond to a shortened delay time in the feedback of model (1). Fig. 8 shows a curve of Hopf bifurcations that demonstrates the sensitivity of an example *m* > 0 solution (*F*₁ = -0.22 Sv) to changes in the feedback. Interestingly, it is stable only for intermediate values of feedback strength and delay time. In both scenarios of increasing and decreasing delay time, the steady-state solution loses stability at a Hopf bifurcation. For the parameters used in Fig. 8, there are curve segments of supercritical Hopf bifurcations (bold), which give rise to stable oscillations, as well as, subcritical Hopf bifurcations (non-bold).
<figure><img src="image_10.png" /><figcaption>Fig. 8. Stability diagram of the steady-state solution at $F_1 = -0.22$ Sv in the $(\tau, \sigma)$-plane. Non bold/bold curves show subcritical/supercritical Hopf bifurcations.</figcaption></figure>

The latter correspond to so-called “dangerous” bifurcations [63], where the system jumps suddenly and irreversibly to a remote state with an alternative circulation pattern. Therefore, the dynamics discussed in this paper could also be relevant in a modern context. Of course, the Drake Passage of modern continental configuration is smaller than the passage of the Drake World setup used here (Fig. 1(a)). Nonetheless, it is worthwhile to consider the possibility that changes in zonal flow, such as that of the ACC, could play a part in destabilising the meridional overturning circulation of the Atlantic Ocean.

### CRediT authorship contribution statement

**Andrew Keane:** Writing – review & editing, Writing – original draft, Visualization, Methodology, Investigation, Formal analysis, Data cura- tion, Conceptualization. **Alexandre Pohl:** Writing – review & editing, Writing – original draft, Visualization, Software, Resources, Methodol- ogy, Investigation, Data curation, Conceptualization. **Henk A. Dijkstra:** Writing – review & editing, Writing – original draft, Methodology, Investigation, Conceptualization. **Andy Ridgwell:** Writing – review & editing, Writing – original draft, Visualization, Software, Resources, Methodology, Investigation, Data curation, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing finan- cial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

This is a contribution to UMR 6282 Biogeosciences Team SEDS and UNESCO project IGCP 735 “Rocks and the Rise of Ordovician Life” (Rocks n’ ROL). AP acknowledges the support of the French Agence Nationale de la Recherche (ANR) under references ANR-22-CE01-0003 (project ECO-BOOST) and ANR-23-CE01-0003 (project CYCLO-SED) and the Leverhulme Trust for supporting the project ‘Earth System dynamics at the dawn of the animal-rich biosphere’ (RPG-2022-233). Calculations were performed using HPC resources from DNUM CCUB (Centre de Calcul de l’Université de Bourgogne). The project has re- ceived funding from the ERC, project 101055096 (HD). AR acknowl- edges support from National Science Foundation grant EAR-2121165, as well as from the Heising-Simons Foundation, United States. This work benefited from a research visit of AK to University of California- Riverside funded by a Charlemont grant from the Royal Irish Academy. We are grateful to two anonymous reviewers for their feedback that helped to improve the quality and presentation of the paper.

## Appendix A. Model equations for negative m

In the original model, without the delayed feedback term and as described in Ref. [45], one set of equations is applied when the ocean flow is clockwise, while another set of equations is applied when the flow is counter-clockwise. The direction of the flow is determined by the model variable m, as defined in Eq. (2) above. Since the oscillatory dynamics that we are interested in occur for positive m, we only provide the model equations for the positive m case in the main text above. When m is negative, the meridional flux between the boxes (shown in Fig. 3) is reversed and Eq. (1) becomes:

$$V \dot{S}_1 = S_0 F_1 - m (S_3 - S_1) + \sigma (S_1^\tau - S_1) \quad (A.1)$$

$$V \dot{S}_2 = -S_0 F_2 - m (S_1 - S_2).$$

## Appendix B. Salt conservation approximation

The salt conservation rule used in the above analysis

$$S_3(t) = 3S_0 - S_1(t) - S_2(t) \quad (B.1)$$

takes the same form as in previous studies of ocean box models [45,47] and has the significant advantage of reducing the model to a relatively simple set of two differential equations (i.e. Eqs. (1)). This rule can be interpreted as saying that any salt from the total $3V S_0$ that is not in box 1 or box 2 must be in box 3. It can be derived by summing the equations of the full three box model, without the delay term,

$$V \dot{S}_1 = S_0 F_1 + m (S_2 - S_1)$$

$$V \dot{S}_2 = -S_0 F_2 + m (S_3 - S_2)$$

$$V \dot{S}_3 = S_0(F_2 - F_1) + m(S_1 - S_3)$$

to find

$$\dot{S}_1 + \dot{S}_2 + \dot{S}_3 = 0$$

then integrating from $t_0$ to t

$$S_1(t) + S_2(t) + S_3(t) - S_1(t_0) - S_2(t_0) - S_3(t_0) = 0$$

$$\Rightarrow S_1(t) + S_2(t) + S_3(t) - 3S_0 = 0$$

$$\Rightarrow S_3(t) = 3S_0 - S_1(t) - S_2(t).$$

However, once the delay term is included in the model, this calcu- lation only holds when $S_1^\tau = S_1$. Otherwise, Eq. (B.1) serves only as an approximation. To derive a more precise salt conservation rule, we begin with the three box model with the delay term

$$V \dot{S}_1 = S_0 F_1 + m (S_2 - S_1) + \sigma (S_1^\tau - S_1) \quad (B.2)$$

$$V \dot{S}_2 = -S_0 F_2 + m (S_3 - S_2)$$

$$V \dot{S}_3 = S_0(F_2 - F_1) + m(S_1 - S_3).$$

One might be inclined to compare our analysis in Section 4 with an analysis of model (B.2). However, numerical continuation of model (B.2) proves awkward, since for $S_1^\tau = S_1$ (i.e. for all fixed points) the Jacobian of the system becomes singular and numerical continuation cannot be applied. Instead, we derive an updated salt conservation rule by summing Eqs. (B.2) to get

$$\dot{S}_1 + \dot{S}_2 + \dot{S}_3 = -\frac{\sigma}{V} (S_1 - S_1^\tau). \quad (B.3)$$

In order to find a salt conservation rule that is relatively easy to work with (e.g. without any derivatives or integrals) we apply the fundamental theorem of calculus to get

$$S_1(t) - S_1(t - \tau) = \frac{d}{dt} \int_{t-\tau}^{t} S_1(t') dt' \quad (B.4)$$

and approximate the integral using the composite trapezoidal rule

$$\int_{t-\tau}^{t} S_1(t') dt' = \frac{0.1\tau}{2} (S_1(t) + 2S_1(t - 0.1\tau) + 2S_1(t - 0.2\tau))$$
<figure><img src="image_11.png" /><figcaption>Fig. B.9. Same as Fig. 5 with more realistic salt conservation rule (Eq. (B.6)).</figcaption></figure>

$$ 
\begin{aligned}
& + 2S_1(t - 0.3\tau) + 2S_1(t - 0.4\tau) + 2S_1(t - 0.5\tau) \\
& + 2S_1(t - 0.6\tau) + 2S_1(t - 0.7\tau) + 2S_1(t - 0.8\tau) \\
& + 2S_1(t - 0.9\tau) + S_1(t - \tau)).
\end{aligned} 
\quad (B.5) $$

Combining Eqs. (B.3)–(B.5) gives

$$ 
\begin{aligned}
\dot{S}_1(t) + \dot{S}_2(t) + \dot{S}_3(t) &= -\frac{\sigma}{V} (S_1(t) - S_1(t - \tau)) \\
&= -\frac{\sigma}{V} \frac{0.1\tau}{2} (\dot{S}_1(t) + 2\dot{S}_1(t - 0.1\tau) + 2\dot{S}_1(t - 0.2\tau) \\
&\quad + 2\dot{S}_1(t - 0.3\tau) + 2\dot{S}_1(t - 0.4\tau) + 2\dot{S}_1(t - 0.5\tau) \\
&\quad + 2\dot{S}_1(t - 0.6\tau) + 2\dot{S}_1(t - 0.7\tau) + 2\dot{S}_1(t - 0.8\tau) \\
&\quad + 2\dot{S}_1(t - 0.9\tau) + \dot{S}_1(t - \tau)).
\end{aligned} $$

This leads to the updated salt conservation rule

$$ 
\begin{aligned}
S_3(t) = 3S_0 - S_1(t) - S_2(t) &- \frac{\sigma}{V} \frac{0.1\tau}{2} (S_1(t) + 2S_1(t - 0.1\tau) \\
& + 2S_1(t - 0.2\tau) + 2S_1(t - 0.3\tau) + 2S_1(t - 0.4\tau) \\
& + 2S_1(t - 0.5\tau) + 2S_1(t - 0.6\tau) + 2S_1(t - 0.7\tau) \\
& + 2S_1(t - 0.8\tau) + 2S_1(t - 0.9\tau) + S_1(t - \tau)).
\end{aligned} 
\quad (B.6) $$

Eq. (B.6) may appear cumbersome, but replacing Eq. (B.1) with Eq. (B.6) in DDE-Bif tool is very simple. It only takes several lines of code to define the additional delay times required for the new salt conservation rule.

We analyse the effect that this new rule has on the results presented above. To keep the comparison brief we show the equivalent of Fig. 5 with the new salt rule, since this figure shows the main effect of the delayed feedback stabilising periodic orbits in a concise way.

Fig. B.9 displays curves of Hopf bifurcations in the $(F_1, \sigma)$-plane for different values of $\tau$ using model (1) with Eq. (B.6). Again, thin and thick curves indicate subcritical and supercritical bifurcations. We observe that, although there are some quantitative differences between Figs. 5 and B.9, the updated salt conservation rule does not impede the role of the delayed feedback in stabilising periodic orbits or changing the criticality of Hopf bifurcations from subcritical to supercritical.

## Appendix C. Steady states of model (1)

Here we include calculations on the linear stability analysis of steady-state solutions to model (1). Due to the nature of DDEs, the characteristic equation is transcendental, which, in this case, limits the analysis to finding a pair of expressions describing the Hopf bifurcations.

First, we nondimensionalise the model by rescaling the dependent variables

$$ \hat{S}_1 = \frac{S_1}{S_0}, $$

$$ \hat{S}_2 = \frac{S_2}{S_0}, $$

$$ \hat{S}_3 = 3 - \hat{S}_1 - \hat{S}_2. $$

We rescale time, $t$, such that

$$ \hat{t} = \frac{kt}{V}, $$

$$ \hat{S}_1^\dagger = \hat{S}_1(\hat{t} - \hat{t}) \text{ with } \hat{t} = \frac{k\tau}{V}. $$

The parameters are rescaled as:

$$ \hat{F}_1 = \frac{F_1}{k}, $$

$$ \hat{F}_2 = \frac{F_2}{k}, $$

$$ \hat{\sigma} = \frac{\sigma}{k}, $$

$$ \hat{\beta} = \beta S_0, $$

$$ \hat{T}^* = \alpha T^*. $$

The resulting nondimensionalised model is

$$ \frac{d\hat{S}_1}{d\hat{t}} = \hat{F}_1 + [\hat{\beta}(\hat{S}_2 - \hat{S}_1) - \hat{T}^*](\hat{S}_2 - \hat{S}_1) + \hat{\sigma}(\hat{S}_1^\dagger - \hat{S}_1) $$

$$ \frac{d\hat{S}_2}{d\hat{t}} = -\hat{F}_2 + [\hat{\beta}(\hat{S}_2 - \hat{S}_1) - \hat{T}^*](\hat{S}_3 - \hat{S}_2). $$

To make further calculations easier, we rewrite the variables as

$$ \hat{S}_1 = S_1, $$

$$ \hat{S}_{21} = \hat{S}_2 - \hat{S}_1, $$

$$ \hat{S}_3 = 3 - 2\hat{S}_1 - \hat{S}_{21}. $$

By letting $\hat{F}_{21} = \hat{F}_2 + \hat{F}_1$, the model can be written as

$$ \frac{d\hat{S}_1}{d\hat{t}} = \hat{F}_1 + (\hat{\beta}\hat{S}_{21} - \hat{T}^*)\hat{S}_{21} + \hat{\sigma}(\hat{S}_1^\dagger - \hat{S}_1) $$

$$ \frac{d\hat{S}_{21}}{d\hat{t}} = -\hat{F}_{21} + 3(\hat{\beta}\hat{S}_{21} - \hat{T}^*)(1 - \hat{S}_1 - \hat{S}_{21}) - \hat{\sigma}(\hat{S}_1^\dagger - \hat{S}_1). $$

For simplicity, we let $\hat{T}^* = 0$, as is done for the results in the main text, and solve for steady-state solutions, $(\hat{S}_1^*, \hat{S}_{21}^*)$, with $\hat{S}_1^* = \hat{S}_1$:

$$ \hat{S}_{21}^* = \pm\sqrt{-\frac{\hat{F}_1}{\hat{\beta}}}, $$

$$ \hat{S}_1^* = 1 - \hat{S}_{21}^* - \frac{\hat{F}_{21}}{3\hat{\beta}\hat{S}_{21}^*}. $$

To consider the stability of the steady states, we linearise the system around $(\hat{S}_1^*, \hat{S}_{21}^*)$:

$$ \begin{pmatrix} \frac{d\hat{S}_1}{d\hat{t}} \\ \frac{d\hat{S}_{21}}{d\hat{t}} \end{pmatrix} = \begin{pmatrix} -\hat{\sigma} & 2\hat{\beta}\hat{S}_{21}^* \\ -3\hat{\beta}\hat{S}_{21}^* + \hat{\sigma} & 3\hat{\beta} - 3\hat{\beta}\hat{S}_1^* - 6\hat{\beta}\hat{S}_{21}^* \end{pmatrix} \begin{pmatrix} \hat{S}_1 \\ \hat{S}_{21} \end{pmatrix} + \begin{pmatrix} \hat{\sigma} & 0 \\ -\hat{\sigma} & 0 \end{pmatrix} \begin{pmatrix} \hat{S}_1^\dagger \\ \hat{S}_{21}^\dagger \end{pmatrix}. $$

As derived in Ref. [46], the characteristic equation of a linear system of n DDEs

$$ \frac{d}{dt}x(t) = Ax(t) + Bx(t - \tau), \quad x \in \mathbb{R}^n, $$

is given by

$$ \det \{ \lambda \mathbb{I} - A - e^{-\lambda \tau} B \}. $$
By evaluating at the steady state, noting that the model is designed for $\dot{S}_{21} > 0$, and letting

$$a = \sqrt{-\beta \hat{F}_1}, \\ b = \beta \hat{F}_2,$$

the characteristic equation becomes

$$\det \begin{pmatrix} \lambda + \hat{\sigma} - \hat{\sigma}e^{-\lambda\hat{\tau}} & -2a \\ 3a - \hat{\sigma} + \hat{\sigma}e^{-\lambda\hat{\tau}} & \lambda + 4a - \frac{b}{a} \end{pmatrix} = 0.$$

We are specifically interested in Hopf bifurcations. Therefore, we set $\lambda = i\hat{\omega}$, where $\hat{\omega}$ corresponds to the nondimensionalised frequency of oscillations around the steady state. By separating the resulting expression into real and imaginary parts, we find the following pair of expressions:

$$0 = 6a^3 + 2a^2\hat{\sigma}(1 - \cos(\hat{\omega}\hat{\tau})) - a\hat{\omega}(\hat{\omega} + \hat{\sigma}\sin(\hat{\omega}\hat{\tau})) - b\hat{\sigma}(1 - \cos(\hat{\omega}\hat{\tau})), \\ 0 = 2a^2(2\hat{\omega} + \hat{\sigma}\sin(\hat{\omega}\hat{\tau})) + a\hat{\omega}\hat{\sigma}(1 - \cos(\hat{\omega}\hat{\tau})) - b(\hat{\omega} + \hat{\sigma}\sin(\hat{\omega}\hat{\tau})).$$

These expressions can be verified with the Hopf bifurcations found numerically in the main text. However, using them for further analysis is impractical and at this point it makes sense to rely on numerical techniques.

## Appendix D. Oscillation dependence on wind speed

The wind speed in cGENIE is physically related to the delay time in model (1). In particular, the zonal flow across the southern channel (i.e. the circumpolar current) is primarily wind-driven. Therefore, the wind speed across the latitudes of the southern channel in cGENIE strongly effects the average zonal flow velocity of the circumpolar current. Concurrently, in model (1) the delay time represents the transportation time of water around the southern pole, averaged across all depths and latitudes of the southern channel. For each latitude, the transportation time is inversely proportional to the average zonal flow velocity. Hence, there is an inverse relationship between the wind speeds across the southern channel in cGENIE and the average transportation time, represented by the delay time in model (1).

To highlight the sensitivity of the stable oscillations to model parameters of cGENIE, we show that the oscillations begin to disappear as the wind speed over the southern channel is slightly increased. Fig. D.10(a) shows the amplitude of the stable oscillations in the cGENIE model, in terms of mean salinity (across the southern channel), of the solution at 6 PAL. The horizontal axis represents the wind speed, across the southern channel, as a factor of the wind speed used in the rest of the paper, as illustrated in Fig. 1(b). In other words, the solution with wind speed equal to one corresponds to the same 6 PAL solution shown in Figs. 2 and 6. As the wind speed is increased we observe that the amplitude of the oscillations decreases and appears to be approaching zero.

Since the circulation around the pole is wind-driven, the wind speed determines the velocity of the flow, and therefore determines the delay time associated with the delayed feedback in the box model. To demonstrate this relationship between wind speed of the cGENIE model and the delay time of the box model, we consider the box model solution shown in Fig. 6(b) at $F_1 = -0.21$ that approximately corresponds to the cGENIE oscillation at 6 PAL. In Fig. D.10(b) we observe that as the delay time is decreased, the amplitude of the oscillations decreases until reaching zero at a Hopf bifurcation.

## Data availability

The code for the version of the muffin release of the cGENIE Earth system model used in this paper, is tagged as v0.9.58, and is assigned a DOI: 10.5281/zenodo.14606515. Configuration files for the specific experiments presented in the paper can be found in the directory: genie-userconfigs/PUBS/submitted/Keane_et_al.PhysicaD.2025. Details of the

<figure><img src="image_12.png" /><figcaption>Fig. D.10. Amplitudes of oscillations with steady-state and periodic solutions shown in black and blue, respectively. (a) cGENIE experiments for 6 PAL with varying wind force strength across the southern channel as a factor of the strength used throughout the paper. (b) Model (1) solutions while varying $\tau$ with $\sigma = 9.5$ Sv and $F_1 = 0.21$ Sv. White filled circle represents Hopf bifurcations.</figcaption></figure>

experiments, plus the command line needed to run each one, are given in the readme.txt file in that directory. All other configuration files and boundary conditions are provided as part of the code release. A manual detailing code installation, basic model configuration, tutorials covering various aspects of model configuration, experimental design, and output, plus the processing of results, is assigned a DOI: 10.5281/zenodo.1407657.

We also provide code for the numerical continuation of model (1): https://github.com/andrewkeane/kpdr_oscillations. It demonstrates the calculation of curves of steady-state solutions and Hopf bifurcations using the Matlab-based software DDE-Biftool.

## References

[1] B.E. Rose, D. Ferreira, Ocean heat transport and water vapor greenhouse in a warm equable climate: A new look at the low gradient paradox, J. Clim. 26 (6) (2013) 2117-2136, http://dx.doi.org/10.1175/JCLI-D-11-00547.1.

[2] A. Pohl, D.A.T. Harper, Y. Donnadieu, G. Le Hir, E. Nardin, T. Servais, Possible patterns of marine primary productivity during the Great Ordovician Biodiversification Event, Lethaia 5 (2) (2017) 187-197, http://dx.doi.org/10.1111/let.12247.

[3] A. Pohl, A. Ridgwell, R.G. Stockey, C. Thomazo, A. Keane, E. Vennin, C.R. Scotese, Continental configuration controls ocean oxygenation during the Phanerozoic, Nature 608 (7923) (2022) 523-527, http://dx.doi.org/10.1038/s41586-022-05018-z.

[4] N. Bouttes, D. Paillard, D.M. Roche, Impact of brine-induced stratification on the glacial carbon cycle, Clim. Past 6 (5) (2010) 575-589, http://dx.doi.org/10.5194/cp-6-575-2010.

[5] A. Pohl, Z. Lu, W. Lu, R.G. Stockey, M. Elrick, M. Li, A. Desrochers, Y. Shen, R. He, S. Finnegan, A. Ridgwell, Vertical decoupling in Late Ordovician anoxia due to reorganization of ocean circulation, Nat. Geosci. 14 (11) (2021) http://dx.doi.org/10.1038/s41561-021-00843-9.

[6] R. Stockey, A. Pohl, A. Ridgwell, S. Finnegan, E.A. Sperling, Decreasing Phanerozoic extinction intensity as a consequence of Earth surface oxygenation and metazoan ecophysiology, Proc. Natl. Acad. Sci. USA 118 (41) (2021) e2101900118, http://dx.doi.org/10.1073/pnas.2101900118.
[7] J. Trabucho Alexandre, E. Tuenter, G.A. Henstra, K.J. Van Der Zwan, R.S. Van De Wal, H.A. Dijkstra, P.L. De Boer, The mid-Cretaceous North Atlantic nutrient trap: Black shales and AOEs, *Paleoceanography* 25 (4) (2010) PA4201, http://dx.doi.org/10.1029/2010PA001925.

[8] K.M. Meyer, L.R. Kump, A. Ridgwell, Biogeochemical controls on photic-zone euxinia during the end-Permian mass extinction, *Geology* 36 (9) (2008) 747–750, http://dx.doi.org/10.1130/G24618A.1.

[9] J.W. Rae, W.R. Gray, R.C. Wills, I. Eisenman, B. Fitzhugh, M. Fotheringham, E.F. Littley, P.A. Rafter, R. Rees-Owen, A. Ridgwell, B. Taylor, A. Burke, Overturning circulation, nutrient limitation, and warming in the Glacial North Pacific, *Sci. Adv.* 6 (50) (2020) eabd1654, http://dx.doi.org/10.1126/sciadv.abd1654.

[10] J.N. Howe, A.M. Piotrowski, T.L. Noble, S. Multiza, C.M. Chiessi, G. Bayon, North Atlantic Deep Water Production during the Last Glacial Maximum, *Nat. Commun.* 7 (2016) 1–8, http://dx.doi.org/10.1038/ncomms11765.

[11] Y. Donnadieu, E. Pucéat, M. Moiroud, F. Guillocheau, J.F. Deconinck, A better-ventilated ocean triggered by Late Cretaceous changes in continental configuration, *Nat. Commun.* 7 (2016) 10316, http://dx.doi.org/10.1038/ncomms10316.

[12] D. Ferreira, J. Marshall, B. Rose, Climate determinism revisited: Multiple equilibria in a complex climate model, *J. Clim.* 24 (4) (2011) 992–1012, http://dx.doi.org/10.1175/2010JCLI3580.1.

[13] M. Laugié, Y. Donnadieu, J.-B. Ladant, L. Bopp, C. Ethé, F. Raisson, Exploring the impact of Cenomanian paleogeography and marine gateways on oceanic oxygen, *Paleoceanogr. Paleoclimatology* 36 (7) (2021) e2020PA004202, http://dx.doi.org/10.1029/2020PA004202.

[14] P. Sepulchre, A. Caubel, J.-B. Ladant, L. Bopp, O. Boucher, P. Braconnot, P. Brockmann, A. Cozic, Y. Donnadieu, V. Estella-Perez, C. Ethé, F. Fluteau, M.-A. Fojouts, G. Gastineau, J. Ghattas, D. Hauglustaine, F. Hourdin, M. Kageyama, M. Khodri, O. Marti, Y. Meurdesoif, J. Mignot, A.-C. Sarr, J. Servonnat, D. Swingoudeuw, S. Szopa, D. Tardif, IPSL-CM5A2. An Earth system model designed for multi-millennial climate simulations, *Geosci. Model. Dev.* 13 (2020) 3011–3053, http://dx.doi.org/10.5194/gmd-13-3011-2020.

[15] M.D. Zelinka, T.A. Myers, D.T. McCoy, S. Po-Chedley, P.M. Caldwell, P. Ceppi, S.A. Klein, K.E. Taylor, Causes of higher climate sensitivity in CMIP6 models, *Geophys. Res. Lett.* 47 (1) (2020) 1–12, http://dx.doi.org/10.1029/2019GL085782.

[16] D.J. Lunt, A. Farnsworth, C. Loptson, G.L. Foster, P. Markwick, C.L. O'Brien, R.D. Pancost, S.A. Robinson, N. Wrobel, Palaeogeographic controls on climate and proxy interpretation, *Clim. Past* 12 (5) (2016) 1181–1198, http://dx.doi.org/10.5194/cp-12-1181-2016.

[17] Z. Sirkes, E. Tziperman, Identifying a damped oscillatory thermohaline mode in a general circulation model using an adjoint model, *J. Phys. Oceanogr.* 31 (8) (2001) 2297–2306, http://dx.doi.org/10.1175/1520-0485(2001)031(<>):2297:iadotm(<>):2.0.co;2.

[18] W. Weijer, H.A. Dijkstra, Multiple oscillatory modes of the global ocean circulation, *J. Phys. Oceanogr.* 33 (11) (2003) 2197–2213, http://dx.doi.org/10.1175/1520-0485(2003)033(<>):2197:MOMOTG(<>):2.0.CO;2.

[19] K.J. Meissner, M. Eby, A.J. Weaver, O.A. Saenko, CO₂ threshold for millennial-scale oscillations in the climate system: Implications for global warming scenarios, *Clim. Dyn.* 30 (2–3) (2008) 161–174, http://dx.doi.org/10.1007/s00382-007-0279-0.

[20] R.J. Haarsma, J.D. Opsteegh, F.M. Selten, X. Wang, Rapid transitions and ultra-low frequency behaviour in a 40 kyr integration with a coupled climate model of intermediate complexity, *Clim. Dyn.* 17 (7) (2001) 559–570, http://dx.doi.org/10.1007/s003820000129.

[21] E. Armstrong, K. Izumi, P. Valdes, Identifying the mechanisms of DO-scale oscillations in a GCM: a salt oscillator triggered by the laurendite ice sheet, *Clim. Dyn.* (2022) 1–19, http://dx.doi.org/10.1007/s00382-022-06564-y.

[22] P.J. Valdes, C.R. Scotese, D.J. Lunt, Deep ocean temperatures through time, *Clim. Past* 17 (4) (2021) 1483–1506, http://dx.doi.org/10.5194/cp-17-1483-2021.

[23] Y. Kuniyoshi, A. Abe-Ouchi, S. Sherriff-Tadano, W.-L. Chan, F. Saito, Effect of climatic precession in Dansgaard-Oeschger-like oscillations, *Geophys. Res. Lett.* 49 (6) (2022) e2021GL095695, http://dx.doi.org/10.1029/2021GL095695.

[24] M. Klockmann, U. Mikolajewicz, H. Kleppin, J. Marotzke, Coupling of the subpolar gyre and the overturning circulation during abrupt glacial climate transitions, *Geophys. Res. Lett.* 47 (21) (2020) e2020GL090361, http://dx.doi.org/10.1029/2020GL090361.

[25] G. Vettoretti, P. Ditlevsen, M. Jochum, S.O. Rasmussen, Atmospheric CO₂ control of spontaneous millennial-scale ice age climate oscillations, *Nat. Geosci.* 15 (4) (2022) 300–306, http://dx.doi.org/10.1038/s41561-022-00920-7.

[26] J. Chen, I.P. Montañez, S. Zhang, T.T. Isson, S.I. Macarewich, N.J. Planavsky, F. Zhang, S. Rauzi, K. Daviau, L. Yao, et al., Marine anoxia linked to abrupt global warming during Earth's penultimate icehouse, *Proc. Natl. Acad. Sci.* 119 (19) (2022) e2115231119, http://dx.doi.org/10.1073/pnas.2115231119.

[27] T.W. Wong Hearing, A. Pohl, M. Williams, Y. Donnadieu, T.H. Harvey, C.R. Scotese, P. Sepulchre, A. Franc, T.R. Vandenbroucke, Quantitative comparison of geological data and model simulations constrains early Cambrian geography and climate, *Nat. Commun.* 12 (1) (2021) 1–11, http://dx.doi.org/10.1038/s41467-021-24141-5.

[28] R. Tostevin, B.J. Mills, Reconciling proxy records and models of Earth's oxygenation during the Neoproterozoic and Palaeozoic, *Interface Focus.* 10 (4) (2020) 20190137, http://dx.doi.org/10.1098/rsfs.2019.0137.

[29] J. Alroy, M. Aberhan, D. Bottjer, M. Foote, F.T. Fürsich, P.J. Harries, A.J.W. Hendy, S.M. Holland, L.C. Ivany, W. Kiessling, M.A. Kosnik, C.R. Marshall, A.J. McGowan, A.I. Miller, T.D. Olszewski, M.E. Patzkowsky, S.E. Peters, L. Villier, P.J. Wagner, N. Bonusco, P.S. Borkow, B. Brenneis, M.E. Clapham, L.M. Fall, C.A. Ferguson, V.L. Hanson, A.Z. Krug, K.M. Layou, E.H. Leckey, S. Nürnberg, C.M. Powers, J.A. Sessa, C. Simpson, A. Tomasovych, C.C. Visaggi, Phanerozoic trends in the global diversity of marine invertebrates, *Science* 321 (5885) (2008) 97–100, http://dx.doi.org/10.1126/science.1156963.

[30] T.W. Dahl, M.L. Siggaard-Andersen, N.H. Schovsbo, D.O. Persson, S. Husted, I.W. Hougård, A.J. Dickson, K. Kjæ r, A.T. Nielsen, Brief oxygenation events in locally anoxic oceans during the Cambrian solves the animal breathing paradox, *Sci. Rep.* 9 (1) (2019) 1–9, http://dx.doi.org/10.1038/s41598-019-48123-2.

[31] J.L. Penn, C. Deutsch, J.L. Payne, E.A. Sperling, Temperature-dependent hypoxia explains biogeography and severity of end-Permian marine mass extinction, *Science* 362 (6419) (2018) eaat1327, http://dx.doi.org/10.1126/science.aat1327.

[32] C. Deutsch, A. Ferrel, B. Seibel, H.-O. Pörtner, R.B. Huey, Climate change tightens a metabolic constraint on marine habitats, *Science* 348 (6239) (2015) 1132–1136, http://dx.doi.org/10.1126/science.aaa1605.

[33] J.G. Rubalcaba, W.C. Verberk, A. Jan Hendriks, B. Saris, H. Arthur Woods, Oxygen limitation may affect the temperature and size dependence of metabolism in aquatic thermosters, *Proc. Natl. Acad. Sci. USA* 117 (50) (2020) 31963–31968, http://dx.doi.org/10.1073/pnas.2003292117.

[34] E. Tziperman, J.R. Toggweiler, Y. Feliks, K. Bryan, Instability of the thermohaline circulation with respect to mixed boundary conditions: Is it really a problem for realistic models? *J. Phys. Oceanogr.* 24 (1994) 217–232, http://dx.doi.org/10.1175/1520-0485(1994)024(<>):2017:10TTCW(<>):2.0.CO;2.

[35] P. Welander, A simple heat-salt oscillator, *Dyn. Atmos. Ocean.* 6 (1982) 233–242, http://dx.doi.org/10.1016/0377-0265(82)90030-6.

[36] A. Colin de Verdière, A simple model of millennial oscillations of the thermohaline circulation, *J. Phys. Oceanogr.* 37 (5) (2007) 1142–1155, http://dx.doi.org/10.1175/JPO3056.1.

[37] A. Ridgwell, J.C. Hargreaves, N.R. Edwards, J.D. Annan, T.M. Lenton, R. Marsh, A. Yool, A. Watson, Marine geochemical data assimilation in an efficient Earth system model of global biogeochemical cycling, *Biogeosciences* 4 (1) (2007) 87–104, http://dx.doi.org/10.5194/bg-4-87-2007.

[38] D. Ferreira, J. Marshall, J.M. Campin, Localization of deep water formation: Role of atmospheric moisture transport and geometrical constraints on ocean circulation, *J. Clim.* 23 (6) (2010) 1456–1476, http://dx.doi.org/10.1175/2009JCLI3197.1.

[39] P. Veervort, S. Kirtland Turner, F. Rochholz, A. Ridgwell, Earth system model analysis of how astronomical forcing is imprinted onto the marine geological record : The role of the inorganic (carbonate) carbon cycle and feedbacks, *Paleoceanogr. Paleoclimatology* 39 (3) (2021) http://dx.doi.org/10.1029/2023PA004826, e2023PA004826.

[40] S.J.V.D. Velde, D. Hülse, C.T. Reinhard, A. Ridgwell, Iron and sulfur cycling in the cGENIE. muffin Earth system model (v0.9.21), *Geosci. Model. Dev.* 14 (2021) 2713–2745, http://dx.doi.org/10.5194/gmd-14-2713-2021.

[41] L. Cao, M. Eby, A. Ridgwell, K. Caldeira, D. Archer, A. Ishida, F. Joos, K. Matsumoto, U. Mikolajewicz, A. Mouchet, J.C. Orr, G.K. Plattner, R. Schlitzer, K. Tokos, I. Totterdell, T. Tschumi, Y. Yamanaka, A. Yool, The role of ocean transport in the uptake of anthropogenic CO₂, *Biogeosciences* 6 (3) (2009) 375–390, http://dx.doi.org/10.5194/bg-6-375-2009.

[42] S. Rahmstorf, Bifurcations of the Atlantic thermohaline circulation in response to changes in the hydrological cycle, *Nature* 378 (6553) (1995) 145–149, http://dx.doi.org/10.1038/378145a0.

[43] H. Alkhayuon, P. Ashwin, L.C. Jackson, C. Quinn, R.A. Wood, Basin bifurcations, oscillatory instability and rate-induced thresholds for Atlantic meridional overturning circulation in a global oceanic box model, *Proc. R. Soc. A* 475 (2225) (2019) 20190051, http://dx.doi.org/10.1098/rspa.2019.0051.

[44] R.M. Van Westen, M. Kliphuis, H.A. Dijkstra, Physics-based early warning signal shows that AMOC is on tipping course, *Sci. Adv.* 10 (6) (2024) eadk1189, http://dx.doi.org/10.1126/sciadv.adk1189.

[45] S. Titz, T. Kuhlbrodt, S. Rahmstorf, U. Feudel, On freshwater-dependent bifurcations in box models of the interhemispheric thermohaline circulation, *Tellus A: Dyn. Meteorol. Ocean.* 54 (1) (2002) 89–98, http://dx.doi.org/10.3402/tellus.a.v54i1.12126.

[46] H.L. Smith, An introduction to delay differential equations with applications to the life sciences, vol. 57, springer New York, 2011, http://dx.doi.org/10.1007/978-1-4419-7646-8.

[47] S. Titz, T. Kuhlbrodt, U. Feudel, Homoclinic bifurcation in an ocean circulation box model, *Int. J. Bifurc. Chaos* 12 (04) (2002) 869–875, http://dx.doi.org/10.1142/S0218127402004759.

[48] K. Engelborghs, T. Luzyanina, G. Samaey, DDE-Biftool: a Matlab package for bifurcation analysis of delay differential equations, 2000, TW Rep. 305.

[49] S.G. Janssens, On a normalization technique for codimension two bifurcations of equilibria of delay differential equations (Master's thesis), Department of Mathematics, Utrecht University, 2010.
[50] J. Sieber, K. Engelborghs, T. Luzyanina, G. Samaey, D. Roosevelt, DDE-BIFTOOL manual — Bifurcation analysis of delay differential equations, 2014, arXiv preprint arXiv:1406.7144.

[51] B. Wage, Normal form computations for delay differential equations in DDE-Bif tool (Master's thesis), Department of Mathematics, Utrecht University, 2014.

[52] K. Pyragas, Continuous control of chaos by self-controlling feedback, Phys. Lett. A 170 (6) (1992) 421-428, http://dx.doi.org/10.1016/0375-9601(92)90745-8.

[53] S. Schikora, H. Wünsche, F. Henneberger, Odd-number theorem: Optical feedback control at a subcritical Hopf bifurcation in a semiconductor laser, Phys. Rev. E 83 (2) (2011) 026203, http://dx.doi.org/10.1103/PhysRevE.83.026203.

[54] A.S. Purewal, C.M. Postlethwaite, B. Krauskopf, Effect of delay mismatch in Pyragas feedback control, Phys. Rev. E 90 (5) (2014) 052905, http://dx.doi.org/10.1103/PhysRevE.90.052905.

[55] S.H. Strogatz, Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering, CRC Press, 2018, http://dx.doi.org/10.1201/9780429492563.

[56] K.A. Crichton, A. Ridgwell, D.J. Lunt, A. Farnsworth, P.N. Pearson, Data-constrained assessment of ocean circulation changes since the middle Miocene in an Earth system model, Clim. Past 17 (5) (2021) 2223-2254, http://dx.doi.org/10.5194/cp-17-2223-2021.

[57] R. Marsh, A. Sobester, E.E. Hart, K. Oliver, N. Edwards, S. Cox, An optimally tuned ensemble of the “eb_gos” configuration of GENIE: Parameter sensitivity and bifurcations in the Atlantic overturning circulation, Geosci. Model. Dev. 6 (5) (2013) 1729-1744, http://dx.doi.org/10.5194/gmd-6-1729-2013.

[58] R. Marsh, A. Yool, T. Lenton, M. Gulamali, N. Edwards, J. Shepherd, M. Krznaric, S. Newhouse, S. Cox, Bistability of the thermohaline circulation identified through comprehensive 2-parameter sweeps of an efficient climate model, Clim. Dyn. 23 (7) (2004) 761-777, http://dx.doi.org/10.1007/s00382-004-0474-1.

[59] C. Li, A. Born, Coupled atmosphere-ice-ocean dynamics in Dansgaard-Oeschger events, Quat. Sci. Rev. 203 (2019) 1-20, http://dx.doi.org/10.1016/j.quascirev.2018.10.031.

[60] A. Roberts, R. Saha, Relaxation oscillations in an idealized ocean circulation model, Clim. Dyn. 48 (7) (2017) 2123-2134, http://dx.doi.org/10.1007/s00382-016-3195-3.

[61] Z. Wang, L.A. Mysak, Glacial abrupt climate changes and Dansgaard-Oeschger oscillations in a coupled climate model, Paleooceanography 21 (2) (2006) http://dx.doi.org/10.1029/2005PA001238.

[62] J.-R. Shi, L.D. Talley, S.-P. Xie, Q. Peng, W. Liu, Ocean warming and accelerating Southern Ocean zonal flow, Nat. Clim. Chang. (2021) 1-8, http://dx.doi.org/10.1038/s41558-021-01212-5.

[63] J.M.T. Thompson, J. Sieber, Predicting climate tipping as a noisy bifurcation: a review, Int. J. Bifurc. Chaos 21 (02) (2011) 399-423, http://dx.doi.org/10.1142/S0218127411028519.