# Exploring the mechanisms of Devonian oceanic anoxia: impact of ocean dynamics, palaeogeography, and orbital forcing

Justin Gérard, Loïc Sablon, Jarno J C Huygh, Anne-Christine da Silva, Alexandre Pohl, Christian Vérard, Michel Crucifix

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

## To cite this version:

Justin Gérard, Loïc Sablon, Jarno J C Huygh, Anne-Christine da Silva, Alexandre Pohl, et al. Exploring the mechanisms of Devonian oceanic anoxia: impact of ocean dynamics, palaeogeography, and orbital forcing. Climate of the Past, 2025, 21, pp.239-260. 10.5194/cp-21-239-2025. hal-04916799

HAL Id: hal-04916799

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

Submitted on 28 Jan 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Exploring the mechanisms of Devonian oceanic anoxia: impact of ocean dynamics, palaeogeography, and orbital forcing

Justin Gérard<fnref n="1" />, Loïc Sablon<fnref n="1" />, Jarno J. C. Huygh<fnref n="2" />, Anne-Christine Da Silva<fnref n="2" />, Alexandre Pohl<fnref n="3" />, Christian Vérard<fnref n="4" />,
and Michel Crucifix<fnref n="1" />

**Correspondence:** Justin Gérard (justin.gerard@uclouvain.be) Received:

28 June 2024 – Discussion started: 1 August 2024

Revised: 2 November 2024 – Accepted: 19 November 2024 – Published: 28 January 2025

**Abstract.** The Devonian is a warmer-than-present geological period spanning from 419 to 359 million years ago (Ma) characterized by multiple identified ocean anoxic/hypoxic events. Despite decades of extensive investigation, no consensus has been reached regarding the drivers of these anoxic events. While growing geological evidence has demonstrated a temporal correlation between astronomical forcing and anoxia during this period, underlying physical mechanisms remain unknown, hence questioning causality. Here, we perform multiple sensitivity experiments, using an Earth system model of intermediate complexity (cGENIE), to isolate the influences of specific Devonian climate and palaeogeography components on ocean oxygen levels, contributing to the better understanding of the intricate interplay of factors preconditioning the ocean to anoxia. We quantify the impact of continental configuration, ocean–atmosphere biogeochemistry (global mean oceanic PO4 concentration and atmospheric pO2), climatic forcing (pCO2), and astronomical forcing on background oceanic circulation and oxygenation during the Devonian. Our results indicate that continental configuration is crucial for Devonian ocean anoxia, significantly influencing ocean circulation and oxygen levels while consistently modulating the effects of other Devonian climate components such as oceanic PO4 concentration, atmospheric pO2 and pCO2, and orbital forcing. The evolution of continental configuration provides a plausible explanation for the increased frequency of ocean anoxic events identified during the Middle and Late Devonian periods, as it contributed to the expansion of oxygen-depleted zones. Our simulations also show that both the decreased atmospheric pO2 and increased oceanic PO4 concentration exacerbate ocean anoxia, consistent with established knowledge. The variation of pCO2 reveals a wide range of ocean dynamics patterns, including stable oscillations, multiple convection cells, multistability, and hysteresis, all leading to significant variations of the ocean oxygen levels and therefore strongly impacting the preconditioning of the ocean to anoxia. Furthermore, multistability and important hysteresis (particularly slow ocean time response) offer different mechanisms to account for the prolonged duration of some ocean anoxic events. Finally, we found that astronomical forcing substantially impacts ocean anoxia by altering ocean circulation and oxygen solubility, with obliquity consistently emerging as the primary orbital parameter driving ocean oxygen variations.

## 1 Introduction

Oceanic anoxic events (OAEs) are defined by the occurrence of poorly oxygenated regional to global oceans. They are pivotal episodes in Earth's history marked by profound disruptions in the global carbon cycle and palaeoceanographic conditions, with critical consequences for marine biodiversity (Vandenbroucke et al., 2015; Becker et al., 2020). These events are generally distinguished by the accumulation of organic-rich black shales in marine sediments, serving as prominent markers of the environmental conditions (Jenkyns, 2010; Becker et al., 2020). The Devonian is a warmer-than-present geological period (Joachimski et al., 2009) spanning from 419 to 359 million years ago (Ma) characterized by significant evolutive steps, such as the progressive colonization by forests of new terrestrial habitats (Stein et al., 2020), the colonization of the oceans by large predatory fish (Dahl et al., 2010), and major mass extinctions (Bond and Wignall, 2008; Kaiser et al., 2016). Across its 60 million years, the Devonian witnessed multiple anoxic/hypoxic events such as the Kellwasser and Hangenberg events (Becker et al., 2020), with estimated durations between 10⁴ and 10⁶ years (Kabanov et al., 2023).

Extensive investigations took place over the last decades but no consensus regarding the mechanisms responsible for triggering these OAEs and explaining their durations has been reached (Haddad et al., 2018). Proposed hypotheses invoke volcanic activity (Paschall et al., 2019; Kabanov et al., 2023), variations in continental weathering through plant evolution (Algeo and Scheckler, 1998), and changes in ocean circulation (Dopieralska, 2009; Chen et al., 2013). Most likely, an interplay of physical and biogeochemical factors must be involved to explain the widespread character, magnitude, and duration of anoxic events. There is also no conclusive evidence suggesting that the triggering factors remained unchanged throughout the evolving conditions of the Devonian period. Recent research provides growing evidence that at least some Devonian OAEs were paced by astronomical cycles (De Vleeschouwer et al., 2014; Wichern et al., 2024). The astronomical forcing refers to

<footnote n="1">¹Earth and Life Institute (ELI), Université catholique de Louvain (UCLouvain), 1348 Louvain-la-Neuve, Belgium</footnote>

<footnote n="2">²Department of Geology, University of Liège, Sart Tilman, 4000 Liège, Belgium</footnote>

<footnote n="3">³Biogéosciences, UMR 6282 CNRS, Université de Bourgogne, 21000 Dijon, France</footnote>

<footnote n="4">⁴Department of Earth Sciences, University of Geneva, 13 Rue des Maraîchers, 1205 Geneva, Switzerland</footnote>
changes in the distribu- tion and amount of incoming solar radiation received by the Earth through the modification of its
position and orienta- tion in space. Some of the Devonian OAEs have been linked to astronomical parameters through recent
findings in differ- ent geological records. Specifically, De Vleeschouwer et al. (2017) and Da Silva et al. (2020) found
maximum relative obliquity power along with eccentricity minimum in multiple δ¹³C records related to the Kellwasser events.
Additionally, De Vleeschouwer et al. (2013) noted that the anoxic black shales in Poland of the *Annulata*, Dasberg, and
Hangen- berg events are separated by approximately 2.4 Myr (million years), corresponding to the expected duration between
two eccentricity nodes. Therefore, astronomical forcing provides a potential explanation for the quasi-periodic character of
Devonian anoxic events. While it is suggested to have played a discernible role in the pacing of these events, astronomi- cal
forcing may not have been enough on its own and would most probably have acted as a perturbation driving a precon- ditioned
system (Carmichael et al., 2019) towards an anoxic regime (e.g. Sarr et al., 2022).

The general objective of the present study is to contribute to the better understanding of the intricate interplay of fac- tors preconditioning the ocean to anoxia, which is crucial for unravelling the mechanisms behind Devonian OAEs. A specific challenge for developing plausible anoxic scenar- ios involving ocean dynamics during the Devonian primarily comes from uncertainties associated with various factors, including atmospheric CO2 concentration, orbital configura- tion, and location of continents (Laskar et al., 2011; Foster et al., 2017; Brugger et al., 2019; Marcilly et al., 2022). To meet our general objective, we conduct various sensitivity analy- ses around those parameters. First, we explore the contribu- tion of Devonian continental configuration and bathymetry to the global mean oceanic oxygen content through modi- fication of ocean circulation patterns (Sect. 3.1). Of partic- ular concern are the uncertainties surrounding palaeogeog- raphy, especially when it comes to bathymetry, as they ex- ert considerable influence on ocean circulation patterns. The bathymetry of the Devonian is largely unknown because the ocean floor has long since disappeared, as the crust (litho- sphere) has been destroyed in subduction zones due to plate tectonics activity (Scotese and Wright, 2018), preventing di- rect access to deep-ocean data. Therefore, sedimentary data almost exclusively concern the shallow ocean, even though growing evidence points to changes in deeper ocean layers, underscoring the potential influence of deep-ocean dynam- ics on Devonian anoxic events (Chen et al., 2013; Boyer et al., 2014; Turner and Slatt, 2016; Song et al., 2017; Zhang et al., 2020a, b). Two distinct approaches exist regarding the representation of bathymetry in Devonian continental re- constructions: one involves a flat bathymetry (Scotese and Wright, 2018), while the other provides a more realistic non-flat bathymetry but with minimal constraints (Vérard, 2019a). The importance of the continental configuration is addressed by comparing results obtained using palaeogeо- graphic reconstructions from Scotese and Wright (2018) with those from Vérard (2019a). Second, we aim to produce an in-depth analysis of ocean dynamics under specific continen- tal configurations (Sect. 3.2). We perform sensitivity analy- ses of three biogeochemical quantities: atmospheric pCO2 (partial pressure of CO2), global mean oceanic PO4 con- centration ([PO4]), and atmospheric pO2 (partial pressure of O2). We also present a hysteresis experiment to explore how the system reacts in transient conditions. Third, we analyse the effect of astronomical forcing on the large-scale ocean circulation and separate the impact of precession, obliquity, and eccentricity to distinguish their respective contribution to anoxia expansion (Sect. 3.3). Together, these three steps al- low us to explore the physical contribution of essential com- ponents of the Devonian climate system to oceanic oxygenation.

The global ocean circulation will be simulated using the Earth system model of intermediate complexity cGENIE (Ridgwell et al., 2007), which has been used in many palaeo and future climate studies. For instance, it was used to investigate the Pacific Meridional Overturning Circulation during the Last Glacial Maximum (Rae et al., 2020), to constrain ocean circulation changes since the middle Miocene (Crich- ton et al., 2021), to explore the contribution of ocean circu- lation reorganization to Late Ordovician anoxia (Pohl et al., 2021), and to understand the causes of the Atlantic Meridional Overturning Circulation slowdown under global warm- ing (Gérard and Crucifix, 2024).

## Methods

### Model description

cGENIE is an Earth system model of intermediate com- plexity following the definition given by Claussen et al. (2002). It has multiple modules that can be classified into two main categories: the physical components and the bio- geochemical components. The physical part of cGENIE con- sists of a frictional geostrophic ocean (3D) coupled with dynamic- thermodynamic sea ice (2D) and a simple energy and moisture balance atmosphere (2D). This forms the C- GOLDSTEIN framework (Marsh et al., 2011). Simulation of the circulation involves advection, convection, and mix- ing. The framework relies on parameterized isoneutral dif- fusion and eddy-induced advection to efficiently transport heat, salinity, and biogeochemical tracers. The oceansurface dynamics are driven by zonal and meridional wind stress based on a prescribed wind field. Surface albedo and wind stress are constant throughout a simulation. The atmospheric energy-moisture balance model (EMBM) uses air temper- ature and specific humidity to compute the heat and mois- ture convergence fields in one effective layer. Precipitation occurs when the relative humidity in this layer exceeds a predefined threshold. The current model version lacks a dy- namic hydrological scheme for continents, setting land evap- oration to zero and allocating precipitation to coastal cells based on a runoff map. Similar to the UVic model of Weaver et al. (2001), sea ice is based on the dynamic-thermodynamic model of Semtner (1976) and Hibler (1979). Prognostic vari- ables such as ice thickness, areal fraction, and concentration are used to track the horizontal transport of sea ice and the exchange of heat and freshwater with both the ocean and the atmosphere (Reinhard et al., 2020). Precipitation goes di- rectly into the ocean, ignoring the ice, and sublimated water from sea ice is directly added to the atmosphere. A full model and coupling procedure description is given in Edwards and Marsh (2005) and Marsh et al. (2011).
Within the ocean, the module BIOGEM, detailed in Ridg- well et al. (2007), manages air–sea gas exchange along with the transformations and spatial redistribution of biogeochem- ical compounds at both surface and interior levels. Biogeni- cally induced chemical fluxes are primarily constrained by nutrient limitation processes (Maier-Reimer, 1993; Ward et al., 2018). Concentration storage of atmospheric composi- tion and useful isotopic properties is achieved through the ATCHEM module that allows for the computation of partial pressure for the different chemical constituents (Ridgwell et al., 2007). cGENIE allows for more modules such as sed- iment geochemistry (SEDGEM), geochronology of erosion (ROKGEM), and plankton species interactions (ECOGEM) (Ward et al., 2018). By integrating these diverse modules, cGENIE enables the exploration and analysis of the interac- tions between physical and biogeochemical processes. To en- sure simplicity and computer efficiency, the configuration of cGENIE used here includes the C-GOLDSTEIN framework together with the BIOGEM and ATCHEM biogeochemical modules.

## Palaeogeographic reconstructions

We consider the continental reconstructions produced by Scotese and Wright (2018) and Vérard (2019a). The palaeo- geographic reconstructions of Scotese and Wright (2018) span the Devonian from 360 to 420 Ma at roughly 5 Myr in- tervals, encompassing a total of 13 distinct configurations. They are inferred from a local determination of the palaeoen- vironment primarily using lithostratigraphic and palaeonto- logical analyses, with a focus on shoreline regions. Local re- constructions are then extrapolated to address temporal un- certainties and geological gaps, ultimately contributing to a global palaeogeographic reconstruction through interpola- tion. This process is repeated across multiple sites world- wide to create a detailed understanding of past environmental conditions. The palaeogeographic reconstructions of Vérard (2019a), on the other hand, are obtained through the con- version of a full-plate tectonic model (global coverage with fully closed plate boundaries; see Vérard, 2019b) into a quan- tified topography (i.e. both orography and bathymetry) fol- lowing Vérard et al. (2015). With the preliminary version of Panalesis (Panalesis.v0, Vérard, 2019a), these reconstruc- tions are available at 370, 383, 393, 408, and 420 Ma. The approach of Scotese and Wright (2018) has spatial gaps, es- pecially in the ocean because there is no constraining infor- mation. Consequently, the deliverable of Scotese and Wright (2018) features a flat ocean bottom, whereas, by contrast, the geophysical approach of Vérard (2019a) generates global re- constructions, with a detailed representation of deep-ocean bathymetry (Vérard, 2019a). Even if the bathymetry from Vérard (2019a) is poorly constrained, it allows for testing the consequences of such non-flat bathymetry for large-scale ocean dynamics. Sea level references of the two types of re- constructions also differ. Scotese and Wright (2018) align the sea level with the modelled period, whereas Vérard (2019a) consistently uses the present-day volume of ocean water, which, in different sizes of oceanic basins, leads to various sea levels. As we are interested in comparing the two se- ries of reconstructions side by side, we adjusted all Vérard (2019a) reconstructions to conform to the sea level specifica- tions made by Scotese and Wright (2018).

In our experiments, the ocean geometry is represented on a 36 36 horizontal equal-area grid with 16 logarithmically spaced vertical levels (from 80 to 765 m). The equal-area fea- ture of the horizontal grid results from the uniform spacing in the longitude and sine of the latitude. Despite its low spatial resolution, cGENIE has proven to be reliable in representing ocean dissolved oxygen spatial patterns and values in modern (Ridgwell et al., 2007) and palaeoclimate (Hülse et al., 2021; Pohl et al., 2021, 2023) simulations. We focus here mainly on the contribution of physical mechanisms and, for this rea- son, consider a scheme of sinking organic matter remineral- ization independent of temperature. It allows us to better iso- late the contribution of physical changes to the biogeochem- ical variables. A simple zonally averaged wind stress profile is applied to the ocean surface. The albedo profile depends on the latitude only (see Fig. 1); there is no explicit depen- dence on the land– sea mask. Although variations in conti- nental configurations are known, in reality, to affect albedo, our experiments are constrained by the absence of a land sur- face scheme (i.e. no snow cover, clouds, or ice sheets), lim- iting the realism, similarly to Cermeño et al. (2022) and Ver- voort et al. (2024). Boundary conditions, including a solar constant, were generated using the “muffingen” open-source software suite v0.9.24 (Ridgwell, 2023) and are shown in Fig. 1 for different continental configurations. In this study, a specific continental configuration is usually referred to as follows: Scot370M (see Fig. 1), where the first half indi- cates the source of the reconstruction (Scotese and Wright, 2018, or Vérard, 2019a) and the second half indicates the age of the configuration (expressed in million years before present). The remaining continental configurations of the De- vonian period are shown at cGENIE resolution in Fig. 2. All the ocean and atmosphere parameters are based on Cao et al. (2009) and are identical to those from Pohl et al. (2022) (excluding the temperature dependence in the organic matter remineralization scheme).

## Description of the experiments

We present three sets of experiments where all simulations were run for 10 000 years until the system reached equi- librium unless stated otherwise. In the first set, we test continental configurations and solar constants (following Gough, 1981) spanning the whole Devonian. This includes all 13 Devonian continental reconstructions from Scotese and Wright (2018) and 5 from Vérard (2019a), which we test with preindustrial values of atmospheric $pCO_2$ (280 ppm) and $pO_2$ (20.95 %) as well as global mean oceanic $[PO_4]$ (2.159 µmol kg$^{-1}$). Throughout this study, these preindus- trial values are also referred to as 1 $pCO_2$, 1 $pO_2$, and 1 $PO_4$, respectively. Additional simulations with the same solar constant across all 13 Scotese and Wright (2018) con- tinental configurations are then used to remove the impact attributed to variations in the solar constant. The second set explores the combinations of different atmospheric and oceanic chemical concentrations applied to six specific conti- nental reconstructions (three from Scotese and Wright, 2018, and three from Vérard, 2019a). The selected reconstructions are given and justified in the Results section corresponding to this experimental set. We used a combination of atmo- spheric $pCO_2$, with
1, 2, 4, 8, and 16 times the preindus- trial value; global mean oceanic [PO₄], with 1, 1.25, 1.5, 1.75, and 2 times the modern concentration; and atmospheric pO₂, with 0.4 times, 0.7 times, and 1 time the preindustrial value. The third set investigates the astronomical configura- tions. Values for pCO₂ and pO₂ were selected to capture a significant portion of the range and uncertainties associ- ated with these parameters throughout the Devonian, with pCO₂ levels from Foster et al. (2017) and Brugger et al. (2019) and pO₂ levels based on Cannell et al. (2022); the range for [PO₄] concentrations was chosen following similar studies on ocean anoxic events using cGENIE (Monteiro et al., 2012; Hülse et al., 2021). We have explored the param- eter space of the astronomical quantities: eccentricity values of 0, 0.03, and 0.06; obliquity values of 21, 22, 23, 24, and 25°; and true longitude of the perihelion ranging from 0 to 315° in increments of 45°. Unless explicitly indicated, all the simulations were performed using an eccentricity of 0 (re- moving any contribution from the precession) and obliquity of 23°. Table 1 provides a schematic overview of the three experiment sets.

## Results

### Set 1: continental configuration

The impact of palaeogeographical evolution on oceanic oxy- genation and ventilation using the continental reconstruc- tions of Scotese and Wright (2018) (red dots) and Vérard (2019a) (blue dots) is shown in Fig. 3. The two sets of re- constructions generate fairly similar mean oxygen concen- trations ([O₂]) but with substantial differences in ventilation ages. Ventilation age is defined as the time since a parcel of water last reached the ocean surface. Simulations using reconstructions of Vérard (2019a) generally present higher ventilation ages because the ocean floor is deeper, increasing the time needed by water masses to reach these depths. In our experimental design, seafloor depth emerges as the primary factor influencing average ventilation age, as no significant differences were observed in the number or extent of deep- water formation regions between the Scotese and Wright (2018) and Vérard (2019a) reconstructions. Additionally, when the ocean floor depth was artificially reduced in the Vérard (2019a) reconstructions, to align with that of Scotese and Wright (2018), the average ventilation age decreased, reaching similar values to those observed with the Scotese and Wright (2018) reconstructions. Figure 3a and b show that the global average dissolved [O₂] and ventilation age evolve in opposite directions (r² 0.91 with *p* value 4.3 10⁻⁷ and r² 0.93 with *p* value 0.008 for the set of reconstruc- tions from Scotese and Wright, 2018, and Vérard, 2019a, re- spectively). This suggests that the variations in oxygen con- centration are primarily explained by changes in ventilation, caused by modifications in ocean circulation dynamics in- duced by the prescribed differences in continental configu- ration. Second-order variations in oxygen observed in our simulations can be attributed to changes in oxygen solubil-

<figure><img src="image_3.png" /><figcaption>Figure 1. Representation of the Devonian bathymetry (m) at approximately 370 Ma according to (a) Scotese and Wright (2018) and (b) Vérard (2019a) at cGENIE resolution. The meridional albedo profile is shown on the right, and the zonally averaged wind stress profile (function of latitude only) is represented by the black arrows in each panel.</figcaption></figure>

<table><caption>Table 1. Experimental design. “Type” indicates whether the continental configuration is based on Scotese and Wright (2018) or Vérard (2019a). pCO₂, [PO₄], and pO₂ are expressed as multiples of preindustrial levels. When variables have multiple values, all possible combi- nations have been tested.</caption><thead><tr><th>Set 1</th><th>Type</th><th>Period (Ma)</th><th>pCO<sub>2</sub></th><th>[PO<sub>4</sub>]</th><th>pO<sub>2</sub></th><th>Eccentricity</th><th>Obliquity</th><th>Precession</th></tr></thead><tbody><tr><td></td><td>Scot</td><td>360; 365;<br>370; 375;</td><td>1</td><td>1</td><td>1</td><td>0</td><td>23°</td><td>0°</td></tr><tr><td></td><td></td><td>380; 385;<br>390; 395;</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>
<table><tbody><tr><td></td><td></td><td>400; 405;<br/>410;<br/>415; 420</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Vera</td><td>370; 383;<br/>393;<br/>408; 420</td><td>1</td><td>1</td><td>1</td><td>0</td><td>23°</td><td>0°</td></tr><tr><td rowspan="3">Set 2</td><td>Type</td><td>Period (Ma)</td><td>pCO<sub>2</sub></td><td>[PO<sub>4</sub>]</td><td>pO<sub>2</sub></td><td>Eccentricity</td><td>Obliquity</td><td>Precession</td></tr><tr><td rowspan="2">Scot</td><td>370; 390;<br/>420</td><td>1; 2; 4;</td><td>1; 1.25; 1.5;</td><td>0.4; 0.7; 1</td><td>0</td><td>23°</td><td>0°</td></tr><tr><td></td><td>8; 16</td><td>1.75; 2</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Vera</td><td>370; 393;<br/>420</td><td>1; 2; 4;</td><td>1; 1.25; 1.5;</td><td>0.4; 0.7; 1</td><td>0</td><td>23°</td><td>0°</td></tr><tr><td></td><td></td><td></td><td>8; 16</td><td>1.75; 2</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">Set 3</td><td>Type</td><td>Period (Ma)</td><td>pCO<sub>2</sub></td><td>[PO<sub>4</sub>]</td><td>pO<sub>2</sub></td><td>Eccentricity</td><td>Obliquity</td><td>Precession</td></tr><tr><td rowspan="2">Scot</td><td>370</td><td>2</td><td>1.25</td><td>0.7</td><td>0; 0.06</td><td>21°; 25°</td><td>0°; 90°;<br/>180°;<br/>270°</td></tr><tr><td>390</td><td>2</td><td>1.25</td><td>0.7</td><td>0; 0.03; 0.06</td><td>21°; 22°;<br/>23°;<br/>24°; 25°</td><td>0°; 45°; 90°;<br/>135°; 180°;<br/>225°;<br/>270°; 315°</td></tr><tr><td></td><td>Vera</td><td>393</td><td>2</td><td>1.25</td><td>0.7</td><td>0; 0.06</td><td>21°; 25°</td><td>0°; 90°;<br/>180°;<br/>270°</td></tr></tbody></table>

ity due to alterations in the solar constant and modifications in exported organic matter (biological pump), which are influenced by temperature adjustments and differences in con- tinental configuration. All the simulations converge towards a stable steady state, except Scot390M. The latter is char- acterized by stable self-sustained oscillations with a period of roughly 3000 years, portrayed by the dashed black line in Fig. 3a and b. Continental changes induce non-monotonic variations in [O₂] through time, with some of them reaching 30 and 50 µmol kg⁻¹ between consecutive equilibrium states using the reconstructions of Scotese and Wright (2018) and Vérard (2019a), respectively. Considering the global scale of

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Representation of the Devonian bathymetry (m) for 16 continental configurations of the Devonian period. Together with Fig. 1, they encompass all 18 reconstructions used in the present study.</figcaption></figure>

these variations, they reflect considerable modifications in the amount of dissolved oxygen in the ocean. Notably, the variations in oxygen intensify with depth and culminate in maximal values within the deepest layers of the ocean. At these depths, the oscillation in Scot390M can reach an amplitude of more than 100 µmol kg⁻¹.

Millennial oscillations in the ocean circulation have been reported in several oceanic models (Sakai and Peltier, 1997; Thornalley et al., 2009; Peltier and Vettoretti, 2014; Sével- lec and Fedorov, 2014; Li and Yang, 2022; Vettoretti et al., 2022). They generally arise from an interplay between advective and convective feedbacks. Given our experimental design, the
emergence of these oscillations could be attributed to either the specific choice of continental configuration or the solar constant. To verify the role of the former, additional simulations with the Scotese and Wright (2018) reconstructions were conducted using the same solar constant (the one from Scot390M). The equilibrium state of these simulations is represented by the light red circles in Fig. 3a. Variations in [O₂] are small, on average lower than 3 µmol kg⁻¹, compared with the corresponding simulations with varying solar constants, suggesting that the solar constant only has a second-order impact on global mean oceanic [O₂]. This also demonstrates that changes in continental configurations can, alone, induce variations in the nature of the solution of the system (apparition of an oscillation). Similar millennial oscillations in the ocean circulation have been reported by Pohl et al. (2022) under early Paleozoic continental configurations. Here we demonstrate that such oscillations can also develop with a Devonian continental configuration.

The extent of the oxygen minimum zone (OMZ) of all the continental reconstructions is shown in Fig. 4a and b. Here,

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

**Figure 3.** Global mean (a) dissolved oxygen concentration (µmol kg⁻¹) and (b) ventilation age (years) of the ocean at equilibrium. Ventilation age is defined as the time since a parcel of water last reached the ocean surface. The red circle and the blue square markers correspond to the simulations conducted using the continental reconstructions of Scotese and Wright (2018) and Vérard (2019a), respectively, with a time-evolving solar constant. The light red colour represents the simulations conducted using the continental reconstructions of Scotese and Wright (2018) with a solar constant of 1322.8 W m⁻² (computed for 390 Ma). The extreme values of the limit cycle observed in the continental configuration of Scot390M are highlighted in black in both panels, representing the amplitude of the oscillations.

the OMZ refers to the horizontal layer of the ocean which has the greatest anoxic extent. The anoxic or dysoxic extents are defined as the percentage of grid cells, within a horizontal layer, with [O₂] under 6.5 and 62.5 µmol kg⁻¹, respectively (after Sarr et al., 2022). Variations in anoxic extent across simulations with different continental reconstruction types (Scotese and Wright, 2018, vs. Vérard, 2019a) are generally substantial: choosing one type of reconstruction over another has a sizeable impact on the simulated OMZ. The re-constructions of Scotese and Wright (2018) tend to produce a shallower and much smaller OMZ, while those of Vérard (2019a) produce a very widespread anoxic extent with the OMZ sometimes being very deep and reaching 20 % of the horizontal layer. The increased anoxic levels from the Vérard (2019a) configurations are attributed to a deeper ocean that increases the average age of water masses and a non-flat bathymetry that introduces more complex circulation patterns, generating additional regions with older water masses. Also, the lack of shallow shelves in the Vérard (2019a) reconstruction reduces upper-ocean deoxygenation and shifts anoxia to the deeper ocean. Nevertheless, a common pattern emerges. For both reconstruction types, the Middle (roughly 394–379 Ma) and Late (roughly 379–360 Ma) Devonian ex-habit higher susceptibility to reaching an anoxic event due to the broader average OMZ extent associated with these periods, and this is consistent with the temporal distribution of observed anoxic events (see Fig. 22.13 in Becker et al., 2020).

Now we focus on the impact of a specific small change in the continental configuration on an anoxia extent. We created a variation of the Vera393M continental configuration, represented in Fig. 5b. The only modification is the removal of the three continental grid points represented by hatching in Fig. 5a. It was designed to alter an important seaway, lowering a common practice in climate modelling (Lawver and Gahagan, 2003; Von Der Heydt and Dijkstra, 2006; Yang et al., 2014). Comparison of Fig. 5a and b shows that the alternative native configuration has a smaller OMZ and is, overall, significantly more oxygenated than the original one. The anoxic extent is reduced by 10 % (from 21 % to 11 %) at the OMZ, here at 1.5 km, but the reduction can increase up to 17 % at 2.5 km depth, withdrawing the entire anoxic region. The increase in oxygen content is associated with a negative water age anomaly. Overall, waters are younger and more oxygenated across the entire interior basin, owing to the advantage of young water masses which were previously blocked by continental grid cells. These simulations are evidence that minor continental changes, in the Devonian, may have a large impact on ocean dynamics and oxygen levels.

Set 2: geochemical state and climate forcing
Contribution of PO₄ and pO₂

Having established the contribution of the continental configuration to variations in oceanic oxygenation during the
Devonian, our next experimental set (Set 2) involves the investigation of the system dynamics under the modification of three prescribed quantities: atmospheric $pCO_2$, global mean oceanic [PO₄], and atmospheric $pO_2$. Here, atmospheric $pCO_2$ is used as a climate forcing, where a doubling of $pCO_2$ induces an increase of 4 W m⁻² in the radiative forcing, which corresponds to a warming of approximately 3 °C in our simulations. For each variable, we test between

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Anoxic ($O_2$ < 6.5 µmol kg⁻¹) and hypoxic ($O_2$ < 62.5 µmol kg⁻¹) extent at the OMZ (expressed as a percentage of the horizontal layer) for the simulations conducted using the continental reconstructions of (a) Scotese and Wright (2018) and (b) Vérard (2019a). Here, the OMZ is defined as the horizontal layer of the ocean which has the greatest anoxic extent, and its depth is given by the black line on the right of each panel. The dashed black line represents the amplitude of the oscillation observed in Scot390M.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Impact of a small variation in the continental configuration on the oceanic dissolved oxygen concentrations (µmol kg⁻¹) at 1.5 km depth. Panels (a) and (b) represent the original and modified Vera393M configuration, respectively. The A and H lines respectively outline the boundaries of the anoxic and hypoxic regions, with thresholds set at 6.5 and 62.5 µmol kg⁻¹ (after Sarr et al., 2022). The continental configuration in (a) is valid at 1.5 km depth. Grid cells representing land at this depth are shaded grey. Hence, it differs slightly from that in Fig. 2, which represents surface bathymetry.</figcaption></figure>

three and five different values, as explained in Sect. 2.3, but implementing these modifications across all 18 continental configurations would result in a large number of runs and data. Therefore, we only selected three reconstructions from Scotese and Wright (2018) and Vérard (2019a) to perform the investigation: one for the Early, Middle, and Late Devonian. We selected periods for which both types of reconstructions are available, as we also aim to compare them side by side. Two natural choices are the continental configurations of 370 and 420 Ma as they are the only two periods matching perfectly across both reconstruction sets. We added the Scot390M configuration due to its oscillatory behaviour, which leads us to choose Vera393M as it corresponds to the closest reconstruction from Vérard (2019a).

The impact of global mean oceanic [PO₄] and atmospheric $pO_2$ on benthic ocean anoxia is shown in Fig. 6 for all the chosen continental configurations (the colours refer to different Devonian periods). The benthic ocean is defined as seafloor grid cells deeper than 2 km below sea level. Throughout our investigations, examining the benthic ocean has proven to be sufficient as every change in dynamics that affects the global ocean was systematically captured by the benthic ocean. In our
experimental setup, the contribution of [PO₄] and atmospheric pO₂ to ocean oxygen content is straightforward and systematic (see Fig. 6). An increase in [PO₄] leads to a greater extent of oceanic anoxia because of an increase in primary productivity and the related remineralization through microbial processes. The atmospheric pO₂ decrease has a significant yet obvious effect, where less oxygen correlates with a greater anoxic extent. Here, the continental configuration modulates the impact of [PO₄] and pO₂ on benthic ocean anoxia, but it does not alter their overall dependency. The percentage of the benthic ocean that has reached the anoxic state is shown in Fig. 7 as a function of atmospheric pCO₂. The role of atmospheric pCO₂ in benthic anoxia is less straightforward than that of pO₂ and [PO₄], as continental configuration greatly conditions the amplitude and even the sign of the effect of a pCO₂ increase. For instance, Scot370M and Vera370M (blue curves in Fig. 7a and b) reveal opposite trends of the anoxic extent with rising pCO₂. Furthermore, for Scot390M, in red in Fig. 7a, we obtain non-monotonous dependencies of anoxia on pCO₂. We further analyse these aspects in the next subsections.

## Contribution of pCO₂

For every continental configuration used in Set 2, the rise in pCO₂ induces two processes with opposite effects on oxygen levels. First, the decrease in oxygen solubility at the surface following a temperature elevation leads to negative oxygen anomalies, which then propagate throughout the benthic ocean by the deep circulation. Second, the northern overturning circulation is systematically strengthened with each doublement of pCO₂, improving the ventilation and reducing the extent of benthic anoxia. We interpret this strengthening as the result of a seesaw effect between the southern and northern water masses, whereby the weakening of the southern convection cell allows for a deeper propagation of the northern cell (see Fig. 8a and b). However, the significance of this seesaw phenomenon varies greatly from one reconstruction to another, leading to differences in results emerging from varying relative contributions of temperature vs. circulation to oceanic [O₂] changes. Specifically, Scot370M (blue curve in Fig. 7a) is the result of a scenario where enhanced ventilation is unable to compensate for the surface oxygen solubility loss, resulting in a net increase in the benthic anoxic extent with pCO₂ (see Fig. 8c). Conversely, Vera370M (blue curve in Fig. 7b) is a case where the strengthening of the overturning exceeds the solubility loss, increasing benthic oxygen levels in response to global climate warming. Throughout our testing, we found that configurations of Scotese and Wright (2018) typically align with the first category, whereas the ones from Vérard (2019a) belong to the second. For the reverse constructions of Scotese and Wright (2018), the alterations of the overturning cell are small, reaching an average strengthening of 1 Sv across the entire ocean, while it often reaches 2 to 3 Sv with those from Vérard (2019a) (see Fig. 7c).

In some cases, effects related to nutrient availability may determine the total effect on oxygen content. For example, in Vera370M, anoxia generally decreases when pCO₂ increases because of the increase in circulation. In some simulations, though, we see the opposite (see Fig. 8e). This occurs at low [PO₄] and when the ocean is already well-oxygenated. As [PO₄] decreases, the benthic anoxic boundary shifts southward (see Fig. 8c) and may reach areas less affected by MOC strengthening. In these simulations, even though circulation strengthens in response to rising temperatures, the active circulation cell does not efficiently ventilate the anoxic areas left in the benthic ocean. Consequently, the main cause of [O₂] change is solubility loss, leading to a subsequent decrease in [O₂]. The same reasoning can be applied to ex-plain the small mismatch between Vera370M and Vera420M (blue and green curves in Fig. 7b) when increasing from 2 to pCO₂.

Another mechanism contributing to ocean oxygenation is associated with the melting of sea ice that occurs when increasing from 1 to 2 pCO₂ (see Fig. 8f). This melting leads to large increases in the number of convection sites and vertical mixing, globally improving the ventilation of the benthic ocean. Although this contribution is almost absent from Fig. 7, it can be seen in Fig. 8e where the major oxygenation occurs when increasing from 1 to 2 pCO₂. Given the strong overturning circulation simulated in Vera370M (see Fig. 8d), the positive oxygen anomaly easily spreads throughout the majority of the benthic ocean, leading to the observed significant reduction in anoxia extent. The effect of sea ice melting is also present for Scot370M but generally does not appear when increasing from 1 to 2 pCO₂. This is because the overturning cell of Scot370M is much weaker than the one from Vera370M (see Fig. 8b and d), which prevents the positive oxygen anomalies from reaching the anoxic water masses, mainly located in the southern part of the ocean (see Fig. 8c). While the influence of sea ice melting on the extent of benthic anoxia is not always striking, it consistently impacts the average benthic oxygen levels, highlighting its significant oxygenating effect. Together, all these mechanisms encompass the relevant processes to understand how large-scale ocean dynamics affect benthic anoxia for Scot370M, Vera370M and Vera420M.

## State transitions in global ocean circulation

We now analyse the remaining curves in Fig. 7 in depth, along with their associated ocean characteristics. As previously mentioned, the anoxic extent in Scot390M (red curve in Fig. 7a) decreases with pCO₂ from 1 to 4 pCO₂ and then increases until 16. We find that the limit cycle identified for this continental configuration (see Sect. 3.1) occurs at 1 and 2 pCO₂. For these concentrations, the system has self-sustained oscillations with a period of 3000 years featuring
<figure><img src="image_8.png" /></figure>

**Figure 6.** Evolution of the benthic anoxic extent (expressed as a percentage of the seafloor) as a function of (a) global mean oceanic [PO₄] and (b) atmospheric pO₂ for different continental reconstructions of the Early, Middle, and Late Devonian. Benthic anoxia is computed here as the percentage of the seafloor deeper than 2000 m characterized by [O₂] < 6.5 µmol kg⁻¹. Simulations were performed using 4 × pCO₂ and 0.4 × pO₂ in (a) and using 4 × pCO₂ and 1.5 × [PO₄] in (b). This pCO₂ was chosen to avoid interference from limit cycles, which only occur at 1 and 2 × pO₂.

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

**Figure 7.** Evolution of the benthic anoxic extent (expressed as a percentage of the seafloor) as a function of pCO₂ for Early, Middle, and Late Devonian simulations conducted using the continental reconstructions of (a) Scotose and Wright (2018) and (b) Vérard (2019a). The green line in panel (a) is dotted to indicate multistability, meaning that this curve is one of many possible trajectories for the evolution of benthic anoxia under varying pCO₂; this specific trajectory fits the experimental design of this set. The dashed red lines in panel (a) represent the range of the stable oscillations. In the presence of a limit cycle, the red square denotes the intermediate value between the extrema of the oscillations. Every simulation was performed using 1.5 × [PO₄] and 0.4 × pO₂. (c) Average northern meridional overturning circulation (MOC) strengthening following a pCO₂ doubling using the continental reconstructions of (a) and (b).

variations in their amplitude (see Fig. 9a and b). These oscillations, with some amplitudes reaching up to 100 µmol kg⁻¹, move the benthic ocean back and forth from oxic to globally anoxic conditions. Nevertheless, the limit cycle shifts the system towards more anoxic conditions on average. Hence, the decrease in oscillation amplitude leads to the reduction of the benthic anoxic extent, explaining the trend observed for Scot390M in Fig. 7a between 1 and 4 pCO₂. Upon reaching 4 times the preindustrial levels, we find the same behaviour as in Scot370M, wherein an increase in pCO₂ induces a reduction of available oxygen at the benthic level through the loss of solubility.

The configuration Scot420M (green curve in Fig. 7a) is most puzzling, with the anoxic extent decreasing, increasing, and decreasing again as pCO₂ levels rise. Yet, upon inspection, the circulation does not undergo large-scale changes that could explain this behaviour. Additional simulations, involving different initial conditions, have revealed that with the Scot420M configuration, the system exhibits pronounced multistability: the system appears to almost randomly converge towards a distribution of potential states. Which state the system reaches depends on both the initial pCO₂ condition and the rate at which we converge toward a specific pCO₂, without clear rules. The range of accessible states also depends on the imposed pCO₂. Specifically, we found
<figure><img src="image_10.png" /></figure>

**Figure 8.** Oceanic circulation and oxygenation for the 370 Ma simulations. (a, b) MOC ($1 \text{ Sv} \equiv 10^6 \text{ m}^3 \text{ s}^{-1}$) for the 370 Ma simulation using the reconstruction of Scotese and Wright (2018) at (a) $1 \times \text{pCO}_2$ and (b) $2 \times \text{pCO}_2$. (c) Benthic [O$_2$] anomaly between 4 and $2 \times \text{pCO}_2$ for Scot370M. The dotted and solid black curves represent the extent of the benthic anoxic region with 2 and $4 \times \text{pCO}_2$, respectively. (d) Same as (b) but for Vera370M. (e) Percentage of the benthic ocean that reached the anoxic state as a function of $\text{pCO}_2$ and [PO$_4$] for Vera370M. The row with $1.5 \times [PO_4]$ is identical to the blue curve in Fig. 7b. (f) Benthic [O$_2$] anomaly between 1 and $2 \times \text{pCO}_2$ for Vera370M. The dotted and solid black curves represent the extent of the benthic anoxic region with 2 and $1 \times \text{pCO}_2$, respectively. Except in (e), all the simulations presented here have been performed using $1.5 \times [PO_4]$ and $0.4 \times \text{pO}_2$.

that the anoxic extent obtained with 8 $\text{pCO}_2$ could vary up to more than 20 % against less than 10 % with 16 $\text{pCO}_2$. Within the spectrum of states that the system can reach, a limit cycle exists, different from the one in Scot390M (see Fig. 9a and b). The identified oscillation has a period of roughly 5000 years and can only be accessed under 2 $\text{pCO}_2$ (see Appendix A for a more detailed explanation). Due to this multistability, Scot420M (green curve in Fig. 7a) only depicts one of the many potential trajectories of benthic anoxia across varying $\text{pCO}_2$ levels (therefore represented as dotted). As a consequence, investigating physical mechanisms with this continental configuration and forcings is particularly challenging.

Vera393M (red curve in Fig. 7b) also shows a non-monotonic behaviour. Between 1 and 8 $\text{pCO}_2$, the warming is linked to an increase in the benthic anoxic extent, whereas going from 8 to 16 $\text{pCO}_2$ results in a reduction of benthic anoxia. Vera393M stands out among the six continental configurations examined in this section: it is the only one which generates a scenario where water masses of southern origin dominate the deep ocean (see Fig. 9d). In other configurations, continental cells in the southernmost latitudinal band prevent the formation of deep water there. The southern circulation relies on a reduced number of convection sites due to the position of the continents. As a result, the ocean exhibits, on average, lower oxygenation levels compared to reconstructions with a dominant northern overturning cell. As $\text{pCO}_2$ increases, the northern circulation cell intensifies, gradually taking over the southern cell because of the see-saw effect (see Fig. 9d-f). Although the northern overturning cell strengthens with each increase in $\text{pCO}_2$, full penetration of related deep waters into the benthic ocean is partially prevented by the southern cell. The main contributions affecting the benthic oxygen levels are the slowdown of the southern cell and the decrease in solubility of surface waters, both contributing to the expansion of the anoxic extent. Only upon going from 8 to 16 $\text{pCO}_2$ do deep waters of northern origin efficiently propagate into the benthic ocean, leading to increased oxygen levels and the reduction of benthic anoxic extent.

## Long time response – hysteresis

The presence of multiple equilibria and state transitions in the global ocean circulation has motivated us to conduct a hysteresis simulation to identify possible transient be-

**Figure 9.** Representation of the limit cycles in (a) benthic [O$_2$] (µmol kg$^{-1}$), benthic temperature (°C), (b) maximum overturning strength (Sv), and northern sea ice extent (%) for different continental configurations and atmospheric $\text{pCO}_2$. The northern sea ice extent is defined here as the
average sea ice cover beyond 60° N. These simulations were run for 20 000 years using 1.25 and 1.75 × [PO₄] for Scot390M and Scot420M, respectively, and 0.7 × pO₂. These pO₂ and [PO₄] values were selected to illustrate the full extent of the influence of the limit cycles on the system. (c) Multistability associated with 8 × pCO₂ (blue) and 16 × pCO₂ (red) considering 10 different initial conditions for Scot420M using 1 PO₄ and 0.4 pO₂. For these simulations, the change from initial to final pCO₂ occurs instantaneously (i.e. in one model iteration), except for the simulations represented by light blue dots, where the change occurs gradually over 5000 years. Measure of the MOC strength in sverdrups (Sv) as a function of latitude and depth under (d) 1 × pCO₂, (e) 4 × pCO₂, and (f) 16 × pCO₂ for Vera393M. These simulations were performed using 1.5 × [PO₄] and 0.4 × pO₂.

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

haviours, potential signs of bistability, and their implications for the ocean time response. Bistability arises when a sys- tem, subject to the same external conditions, can settle into one of two distinct stable equilibria depending on its initial conditions. The hysteresis simulation is performed with the Scot390M continental configuration over 150 000 years and involves a linear increase and decrease in pCO₂ between 1 and 16 times the preindustrial levels as shown in Fig. 10a. Figure 10b shows the transient benthic [O₂] evolution for in- creasing (blue) and decreasing (red) pCO₂.

Given that the simulation restarts from a stable state (pre- vious spin-up of 10 000 years) and considering the slow evo- lution of the forcing (0.056 ppm yr⁻¹), the system can be assumed to be almost in equilibrium. We further show the stable steady state of the system for specific pCO₂ values with dots and stars in Fig. 10b, starting from the blue and red curves, respectively. When the system reaches stable oscilla- tions, the dots and stars are placed at the intermediate value between the extremes of the limit cycle. Comparison of these markers reveals no clear signs of bistability, yet discrepan- cies arise for 2 × pCO₂, where the blue and red markers do not coincide. Under 2 pCO₂, the system converges towards stable oscillations when starting from warm conditions and to a stable steady state when starting from cold conditions.

For the Scot390M configuration, oscillations manifest themselves for specific pCO₂ values within an interval sur-
rounding 2 pCO₂. Additional simulations have located the lower and upper boundary of this interval around 1.8 and 2.2
pCO₂, respectively. Our findings indicate that to trigger these oscillations, atmospheric pCO₂ must fall within this range but
that the system also needs to come from warmer conditions (i.e. higher pCO₂). Hence, we qualify the lower boundary of this
oscillation interval as asymmetrical, indi- cating that it only acts as a boundary when the limit cycle is already present. This
feature of the oscillation interval ac- counts for the emergence of two distinct stationary regimes reached under 2 pCO₂, as
shown in Fig. 10b. In addition, benthic [O₂] displays significant hysteresis, indicating a very slow ocean time response when
the CO₂ varies.

**Figure 10.** Results of a hysteresis experiment. (a) Schematic of atmospheric *pCO*₂ forcing for the hysteresis experiment. The blue and red lines show the transient increase and decrease in the forcing, respectively. (b) Hysteresis curves of average benthic [O*₂*] (µmol kg⁻¹) as a function of *pCO*₂ for Scot390M with 1.25 × [PO*₄*] and 0.7 × pO₂. The x axis is logarithmic. The blue and red curves represent the transient evolution of the average benthic [O*₂*] when the *pCO*₂ forcing increases and decreases, respectively. The blue dots and red stars are the stable equilibrium of the system. The dashed black lines represent the amplitude of the oscillations observed for these stable equilibria. Under 2 × pCO₂ the oscillation manifests only for the red star and not the blue dot.
<figure><img src="image_12.png" /></figure>

### 3.3 Set 3: astronomical forcing

The final set of experiments (Set 3) explores the astronomical forcing contribution to variations in dissolved [O₂] through changes in large-scale ocean dynamics. Given the wide range of considered astronomical configurations (see Sect. 2.3), we limited our analysis to the Scot370M, Scot390M, and Vera393M continental configurations. These reconstructions were chosen as they encompass all the relevant dynamics identified in the previous section, namely sea ice melting, the presence of two overturning cells, and oscillations.

Figure 11 shows how the dysoxic extent of the benthic ocean varies across a range of extreme astronomical configurations for Scot370M and Vera393M. These reconstructions were selected arbitrarily to compare the astronomical signal between significantly different continental configurations. The fixed atmospheric pO₂ is high enough to almost entirely remove the benthic anoxic extent, leaving significant variations only in the dysoxic extent. For Scot370M, the results indicate that obliquity has the largest effect on the benthic dysoxic extent, increasing it up to 13% (see LOLE and HOLE in Fig. 11a). The predominance of the obliquity is not so surprising because obliquity has, among the astronomical parameters, the greatest impact on temperature at high latitudes, where deep-oceanic convection takes place. In comparison, eccentricity and precession play a relatively minor role, resulting in maximal variations of around 2% of the dysoxic extent in Fig. 11a; for intermediate values of obliquity (not shown here), this impact can increase to 5%. For Vera393M, the astronomical forcing has a smaller, yet still significant, impact on the benthic dysoxic extent of poorly oxygenated benthic waters, about half the amplitude identified for Scot370M. Across all astronomical configurations considered, the maximal difference in the benthic dysoxic extent is 15.6% for Scot370M against 8.5% for Vera393M. Still, the response remains dominated by obliquity, with a lesser contribution from eccentricity and precession. These results also highlight that the continental configuration can modulate the influence of the astronomical forcing on the system. This modulating aspect was already introduced in Sarr et al. (2022) for the Late Cretaceous, where some astronomical configurations could induce up to 50% of the water volume to become anoxic in enclosed and already poorly oxygenated basins (central Atlantic) compared to 0% in more open basins (south of the Walvis Ridge).

Figure 12 shows the influence of the astronomical forcing on the nature of the solution for Scot390M expressed as eccentricity–precession maps for each obliquity value considered. Each dot represents a simulation at equilibrium for a certain astronomical configuration, with the colour representing the type of equilibrium reached by the system: red for a stable steady state and blue for stable oscillations. Results in Fig. 12 demonstrate that each astronomical parameter can shift the equilibrium type of the system, transitioning from stable oscillations to a stable steady state or vice versa. Obliquity and eccentricity are the main factors determining the equilibrium type, with precession being secondary. The importance of obliquity and eccentricity is well-explained by their influence on annual mean temperature beyond 60° N. Figure 12 also shows the range of astronomical configurations for which stable oscillations occur. At low obliquity and low eccentricity, the high latitudes become too cold to support stable oscillations, providing a cold boundary for the oscillation. At higher obliquity, stable oscillations also
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Percentage of the benthic ocean that reached anoxia and dysoxia under different astronomical configurations for Scot370M (blue) and Vera393M (red). Benthic dysoxia is computed here as the percentage of the seafloor deeper than 2000 m characterized by $O_2 < 62.5 \mu mol kg^{-1}$. LO (low obliquity) and HO (high obliquity) correspond to simulations where obliquity is set to 21 and 25°, respectively. LE (low eccentricity) and HE (high eccentricity) correspond to simulations where eccentricity is fixed at zero to re- move any contribution from precession and 0.06, respectively. In HE simulations, precession effects are maximized, with the differ- ent months indicating when in the year the perihelion is reached. All simulations were performed using $2 \times pCO_2$, $1.25 \times [PO_4]$, and $0.7 \times pO_2$.</figcaption></figure>

tend to disappear, although the threshold depends on the pre- cise value of eccentricity, with precession becoming critical only when obliquity is at 23°. For these configurations, the high latitudes become too warm, indicating a warm boundary for the oscillation. Additional simulations have shown that the oscillation interval also features the asymmetrical nature identified in the hysteresis section: oscillations emerge only when the system is started from initial conditions beyond the warm boundary.

## Discussion

### Implications for current hypotheses about Devonian anoxic events

We used cGENIE to investigate the contribution of continen- tal configuration, geochemical states, and astronomical forc- ing to Devonian anoxia. First, we found that the continen- tal configuration plays a central role throughout our study, not only through its direct impact on ocean dynamics and oxygen levels, but also through modulating the influences of $pCO_2$, oceanic $[PO_4]$, $pO_2$, and astronomical forcings. The differences between reconstructions from Scotese and Wright (2018) and Vérard (2019a) underscore the impor-

tance of bathymetry for circulation dynamics and anoxic ex- tent. Notably, the presence of self-sustained oscillations in the Scot390M configuration emphasizes the nonlinear na- ture of the system's response to small changes in continental geometry. These results raise questions, similarly to Wong Hearing et al. (2021), regarding the assessment of the reli- ability of different palaeogeographic reconstructions in the deep geological past. Our findings also reveal that continen- tal configurations of the Middle and Late Devonian epochs are associated with increased OMZ extents (as opposed to those from the lower Devonian), thereby preconditioning the system for ocean anoxia. Therefore, in our experimen- tal design, the continental configuration offers an explana- tion for the enhanced frequency of detected anoxic events during these periods (see Fig. 22.13 in Becker et al., 2020). While our research identified continental configuration as an essential mechanism driving ocean $[O_2]$ in the Devonian, previous studies (modelling or geological outcomes) have reached similar conclusions for other periods. Donnadieu et al. (2016) investigated, using the Fast Ocean Atmosphere Model (FOAM), the Late Cretaceous climate and empha- sized the potentially important role of palaeogeography as a preconditioning factor to the development of OAEs. Pohl et al. (2022) showed, with cGENIE, how continental rear- rangements, even small, during the Phanerozoic eon led to profound variations in ocean oxygenation through reorgani- zations of large-scale ocean circulation. Recent studies on ni- trogen proxies related to anoxia extent (Percival et al., 2022), along with more accurate dating (Hedhli et al., 2023) of anoxic events, have emphasized the significance of local fac- tors. They suggest that the Hangenberg black shales do not represent a single global anoxic event but rather multiple, di- achronous basins similar to the Black Sea, occurring around the world. To explain this, they propose that the tectonic con- figuration near the Devonian–Carboniferous boundary could have led to the creation of restricted marine basins, which in turn caused diachronous anoxic events associated with the Hangenberg extinction.

Then, we verified the profound yet straightforward effect of atmospheric $pO_2$ and global mean oceanic $[PO_4]$ on the benthic anoxic extent. As a general rule, both an increase in $[PO_4]$ and a decrease in atmospheric $pO_2$ lead to a greater extent of oceanic anoxia. This aligns with Percival et al. (2020), who propose that the upper Kellwasser, *Annulata*, and Hangenberg
events were triggered by an influx of ter- restrial phosphorus from enhanced weathering. Our analysis also reveals that continental configuration significantly influ- ences the amplitude and even sometimes the sign of the con- tribution of a pCO₂ increase to benthic oxygen levels. This is due to the relative importance of two competing mech- anisms: surface solubility loss and overturning strengthen- ing. Studies generally suggest that OAEs of the Devonian occurred under warm conditions; however, some research in- dicates that these events could have also taken place during periods of a relatively cold global climate (Percival et al.,

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Eccentricity–precession maps representing the nature of the solution of the system across a range of obliquity values: blue for stable oscillations and red for a stable steady state. Each dot represents a simulation conducted with a different astronomical configuration, where eccentricity is given by the radius and the longitude of perihelion by the angle on the circle (the corresponding month is provided as Jan: January, Feb: February, etc.). All simulations were performed for Scot390M using 2 × pCO₂, 1.25 × [PO₄], and 0.7 × pO₂.</figcaption></figure>

2020; Kabanov et al., 2023). Our results show that depending on the continental configuration, both are valid hypotheses. Additionally, we identified a wide range of ocean dynam- ics with substantial effects on ocean oxygen levels includ- ing self-sustained oscillations, variations in overturning cell dominance, modification of sea ice cover, emergence of mul- tistability, and prolonged ocean time response indicated by a hysteresis. These different behaviours are associated with significant variations in global oxygen content. Hence, the richness of the ocean dynamics observed here may be rele- vant for understanding the dynamics of anoxia. For instance, the slow time response of the system under the gradual pCO₂ in the hysteresis simulation provides insight into the pro- longed duration of certain OAEs. Applying constraints on at- mospheric pCO₂ and surface temperature, such as those pre- sented in Joachimski et al. (2009) and Becker et al. (2020), would help to better isolate the relevant ocean dynamics for individual events.

Finally, we assessed the role of astronomical forcing in driving variations in benthic anoxia through changes in large- scale ocean dynamics and oxygen solubility for different continental configurations. Obliquity consistently emerges as the dominant orbital parameter influencing benthic ocean anoxia due to its effect on the annual mean temperature at high latitudes, where deep-ocean convection occurs. In our results, high obliquity systematically led to more ex- tensive oxygen-depleted zones, with eccentricity and preces- sion playing a secondary role. Several studies (Meyers et al., 2012; De Vleeschouwer et al., 2017; Da Silva et al., 2020) indicate that obliquity increases prior to the anoxic phases of Kellwasser events and Oceanic Anoxic Event 2 (in the Late Cretaceous). The results also revealed that the continental configuration significantly modulates the influence of astro- nomical forcing on ocean anoxia. Additionally, astronomical forcing can shift the system’s equilibrium, transitioning be- tween stable oscillations and steady states, with obliquity and eccentricity being the primary determinants of these equilib- rium states.

## Model limitations

cGENIE belongs to the category of Earth system models of intermediate complexity. Its spatial resolution is coarse, the hydrological scheme is simplified, and the dynamics of atmospheric currents are not explicitly modelled. Reassur- ingly, cGENIE has repeatedly demonstrated its ability to cap- ture large-scale ocean circulation patterns (Rahmstorf et al., 2005;
Pohl et al., 2022; Gérard and Crucifix, 2024), providing con- fidence in the robustness of the simulated global circula- tion. Furthermore, even with its coarse spatial resolution, the model effectively captured the variations between continen- tal reconstructions to produce very distinct ocean dynamics regimes. Yet, model limitations have some consequences re- garding how we should interpret the implications of our results about the anoxia observed during the Devonian.

First, the phenomenon of multistability might have been exacerbated. It has long been suggested that multiple ac- tive oceanic circulation solutions may coexist with given boundary conditions, at least in somewhat idealized ocean geometries (Hughes and Weaver, 1994; Rahmstorf, 1995). However, authors have voiced the concern that more de- tailed coastline geometry, proper atmospheric feedbacks, and realistic atmospheric variability may, in the end, settle the system in one dominant regime. Nevertheless, considering the modelling work of Ferreira et al. (2011), Brunetti et al. (2019), Margazoglou et al. (2021), Eberhard et al. (2023), and Lohmann et al. (2024) as well as theoretical considera- tions from Weijer et al. (2019), multistability and limit cy- cles remain plausible scenarios. Furthermore, as discussed in Pohl et al. (2022), Valdes et al. (2020) seem to obtain similar millennial oscillations in HadCM3 as those observed in our study.

Next, in its present configuration, cGENIE does not in- clude ice sheet dynamics, with the implication that the sea level is constant. Ice sheets may also have grown to enough extent (De Vleeschouwer et al., 2015; Becker et al., 2020; Dahl et al., 2022) to release some water into the ocean, po- tentially impacting the circulation. Finally, cGENIE's reso- lution is too coarse to properly capture the connection be- tween the deep ocean and epicontinental margins. Yet, the development of deep-ocean anoxia would help explain the synchronous character of anoxia observed in distant geolog- ical transects. However, the pathways by which deep-ocean anoxia may propagate into epicontinental seas or influence the conditions near arc islands require additional investigation.

Due to computational constraints and the extensive sim- ulations required for each configuration, our study focused only on two distinct types of continental reconstructions. We would naturally encourage an investigation using a broader range of distinct continental configuration types. Additional insights may emerge by considering hybrid continental re- constructions, which would combine continental locations and shallow waters from Scotese and Wright (2018) with the deep ocean from Vérard (2019a). Such hybrid config- urations were explored by Cermeño et al. (2022), merg- ing the continental reconstructions from Scotese and Wright (2018) and Merdith et al. (2021). To our knowledge, these approaches still remain largely unexplored and we advocate for further investigations using different climate models, par- ticularly EMIC and Coupled Model Intercomparison Project-class models to elucidate similarities and potential disparities between cGENIE and finer models.

## Conclusions

This study investigates Devonian OAEs, crucial episodes marked by global disruptions in the carbon cycle and severe impacts on marine biodiversity, with significant uncertainties surrounding their timing, duration, and driving mechanisms. Using the cGENIE Earth system model, we carried out sys- tematic sensitivity experiments to isolate the contributions of specific Devonian climate and palaeogeography components to ocean oxygen levels, providing new insights into these ocean anoxic events. The key findings of our paper can be summarized as follows.

Continental configuration consistently modulates the signal of all other considered components of the De- vonian climate (oceanic $[PO_4]$, atmospheric $pO_2$, and $pCO_2$ as well as orbital forcing) on anoxia, highlighting its crucial contribution to ocean anoxic events. Further- more, the continental configuration greatly influences ocean circulation and global oxygen content, with, at one point, the development of stable oceanic oscilla- tions with periods ranging from 3000 to 5000 years. The oscillations, on average, are found to significantly reduce the oxygen concentration in the deep ocean. The evolution of the continental configuration through the Devonian thus offers a plausible explanation for the enhanced frequency of reported anoxic events dur- ing the Middle and Late Devonian periods by creating more extensive oxygen- depleted zones. Our investiga- tion also revealed the presence of pronounced multi- stability (more than one stable state of the system ex- ists for the same climatic configuration) for a specific continental configuration, providing a potential mech- anism to account for the prolonged duration of these events (“lockdown” of the system in a state with strong oxygen-depleted conditions).

The influence of global mean oceanic $[PO_4]$ and atmo- spheric $[O_2]$ on ocean oxygen content is in line with current knowledge. Simulations indicate that both de- creased atmospheric $[O_2]$ and increased oceanic $[PO_4]$ exacerbate ocean anoxia. The variation of atmospheric $pCO_2$ (here used as a climate forcing to modify tem- peratures) revealed important changes in ocean anoxia through modification of oxygen solubility, ocean circu- lation, and sea ice extent. The results indicate a range of ocean dynamics patterns, including stable oscillations, multiple convection cells, and hysteresis (dependence of the state of a system on its history), all highly sen- sitive to continental configuration. The hysteresis fea- ture also provides a plausible mechanism to explain the prolonged duration of some ocean anoxic events. Furthermore, these processes and behaviours lead to impor- tant modifications of the ocean oxygen levels, therefore strongly impacting the preconditioning of the ocean to anoxia.

Astronomical forcing substantially impacts ocean anoxia by altering ocean circulation and oxygen sol- ubility. Obliquity consistently emerges as the primary orbital parameter driving the vast majority of observed oceanic oxygen variations due to its strong annual aver- age influence on high latitudes. In our experiments, high obliquity persistently led to more extensive oxygen- depleted zones. Simulations demonstrate that astronom- ical forcing can initiate stable oscillations
within the system. The modulation of the orbital forcing impact on anoxia by continental configuration could explain why some anoxic events correlate with astronomical cy- cles, whereas others do not, highlighting the temporal variability in the astronomical pacing of ocean anoxia throughout the Devonian.

This study underscores that shifts in ocean circulation, in- fluenced by continental configuration, CO₂ levels, and orbital forcing, exert significant impacts on ocean anoxia. Ocean cir- culation revealed surprisingly complex dynamics, with re- sponse times ranging from years to tens of thousands of years, which are relevant for understanding anoxia in the De- vonian and, more generally, throughout the Phanerozoïc. Our results also highlight the complexity of modelling Devonian ocean anoxic events and the fact that much work remains to reach complete comprehension of these critical episodes in Earth’s history.

## Appendix A: Limit cycle in Scot420M

Our goal is to explore the unusual shape of the limit cy- cle observed in the Scot420M continental setup. Similar to Scot390M, we find that these oscillations arise from the inter- play between advective and convective feedbacks. However, differences between the two oscillations suggest that another mechanism is at play for Scot420M (see Fig. 9a and b). In this configuration, the observed oscillations emerge as a com- posite of a similar limit cycle as in Scot390M and a deep- decoupling oscillation. This additional deep decoupling os- cillation involves the gradual accumulation of heat below the surface, as shown in Fig. A1. Mack (2021) have also ob- served analogous deep decoupling oscillations, manifesting as two quasi-stable states – an unventilated phase and a ven- tilated phase – using the MITgcm ocean model for Devo- nian palaeoclimate conditions. Furthermore, the close simi- larity between Fig. A1 and Fig. 3.8 from Mack (2021) fur- ther demonstrates the ability of cGENIE to capture com- plex large-scale ocean dynamics patterns. The combination of these two limit cycles generates this “two-step” pattern of the oscillations present for the Scot420M continentalconfig- uration (see Fig. 9a and b).

<figure><img src="image_15.png" /><figcaption>**Figure A1.** Representation of the mean vertical (a) temperature and (b) density profile for Scot420M. The simulation has been performed using 2 × pCO₂, 1.75 × [PO₄] and 0.7 × pO₂.</figcaption></figure>

**Code availability.** The code for the version of the “muf- fin” release of the cGENIE Earth system model used in this paper is tagged as v0.9.59 and is assigned the follow- ing DOI: https://doi.org/10.5281/zenodo.14680971 (Ridgwell et al., 2025). Configuration files for the specific experiments pre- sented in the paper can be found in this directory: genie- userconfigs/PUBS/published/ Gerard_et_al_2024. Details of the ex- periments, 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. Amanual detailing code installation, basic model configuration, tuto- rials covering various aspects of model configuration, experimental design, and output, plus the processing of results, is assigned the fol- lowing DOI: https://doi.org/10.5281/zenodo.13377225 (Ridgwell et al., 2024).
**Data availability.** No datasets were used in this article.

**Author contributions.** JG conducted the cGENIE experiments, developed the code for the analysis, performed all computations, and drafted the manuscript. MC, ACDS, and AP contributed to the study design. CV provided continental reconstructions. All authors

contributed to the discussion of results and writing the manuscript, and all authors approved the final version.

**Competing interests.** The contact author has declared that none of the authors has any competing interests.

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

**Acknowledgements.** Computational resources have been provided by the supercomputing facilities of the UCLouvain (CIS- M/UCL) and the Consortium des Équipements de Calcul Intensif en Fédération Wallonie Bruxelles (CÉCI). We are grateful to the editor and reviewers for constructive comments during the editing process.

**Financial support.** This research is funded by the Belgian National Fund of Scientific Research (FNRS), project WarmAnoxia, PDR grant no. T.0037.22. Anne-Christine Da Silva was supported by FNRS CDR grant no. J.0037.21. Alexandre Pohl was supported by 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 programme Tellus of the Institut National des Sciences de l'Univers, CNRS (project ROSETTA).

**Review statement.** This paper was edited by Yannick Donnadieu and reviewed by two anonymous referees.

## References

Algeo, T. J. and Scheckler, S. E.: Terrestrial-marine teleconnexions in the Devonian: links between the evolution of land plants, weathering processes, and marine anoxic events, Philos. T. Roy. Soc. B, 353, 113-130, https://doi.org/10.1098/rstb.1998.0195, 1998.

Becker, R., Marshall, J., Da Silva, A.-C., Agterberg, F., Gradstein, F., and Ogg, J.: The Devonian Period, in: Geologic Time Scale 2020, Elsevier, 733-810, https://doi.org/10.1016/B978-0-12-824360-2.00022-X, 2020.

Bond, D. P. and Wignall, P. B.: The role of sea-level change and marine anoxia in the Frasnian–Famennian (Late Devonian) mass extinction, Palaeogeogr. Palaeoclim., 263, 107-118, https://doi.org/10.1016/j.palaeo.2008.02.015, 2008.

Boyer, D. L., Haddad, E. E., and Seeger, E. S.: The last gasp: trace fossils track deoxygenation leading into the Frasnian–Famennian extinction event, PALAIOS, 29, 646-651, https://doi.org/10.2110/palo.2014.049, 2014.

Brugger, J., Hofmann, M., Petri, S., and Feulner, G.: On the Sensitivity of the Devonian Climate to Continental Configuration, Vegetation Cover, Orbital Configuration, CO₂ Concentration, and Insolation, Paleoceanography and Paleoclimatology, 34, 1375-1398, https://doi.org/10.1029/2019PA003562, 2019.

Brunetti, M., Kasparian, J., and Vérard, C.: Co-existing climate at- tractors in a coupled aquaplanet, Clim. Dynam., 53, 6293-6308, https://doi.org/10.1007/s00382-019-04926-7, 2019.

Cannell, A., Blamey, N., Brand, U., Escapa, I., and Large, R.: A revised sedimentary pyrite proxy for atmospheric oxygen in the Paleozoic: Evaluation for the Silurian-Devonian Carboniferous period and the relationship of the results to the observed biosphere record, Earth-Sci. Rev., 231, 104062, https://doi.org/10.1016/j.earscirev.2022.104062, 2022.

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

Carmichael, S. K., Waters, J. A., Königshof, P., Suttner, T. J., and Kido, E.: Paleogeography and paleoenvironments of the Late Devonian Kellwasser event: A review of its sedimentological and geochemical expression, Global Planet. Change, 183, 102984, https://doi.org/10.1016/j.gloplacha.2019.102984, 2019.

Cermeño, P., García-Comas, C., Pohl, A., Williams, S., Benton, M. J., Chaudhary, C., Le Gland, G., Müller, R. D., Ridgwell, A., and Vallina, S. M.: Post-extinction recovery of the Phanerozoic oceans and biodiversity hotspots, Nature, 607, 507-511, https://doi.org/10.1038/s41586-022-04932-6, 2022.

Chen, D., Wang, J., Racki, G., Li, H., Wang, C., Ma, X., and Whalen, M. T.: Large sulphur isotopic perturbations and oceanic changes during the Frasnian–Famennian transition of the Late Devonian, J. Geol. Soc. London, 170, 465-476, https://doi.org/10.1144/jgs2012-037, 2013.

Claussen, M., Mysak, L., Weaver, A., Crucifix, M., Fichefet, T., Loutre, M.-F., Weber, S., Alcamo, J., Alexeev, V., Berger, A., Calov, R., Ganopolski, A., Goosse, H., Lohmann, G., Lunkeit, F., Mokhov, I., Petoukhov, V., Stone, P., and Wang, Z.: Earth system models of intermediate complexity: closing the gap in the spectrum of climate system models, Clim. Dynam., 18, 579-586, https://doi.org/10.1007/s00382-001-0200-1, 2002.

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

Dahl, T. W., Hammarlund, E. U., Anbar, A. D., Bond, D. P. G., Gill, B. C., Gordon, G. W., Knoll, A. H., Nielsen, A. T., Schovsbo, N. H., and Canfield, D. E.: Devonian rise in atmospheric oxygen correlated to the radiations of terrestrial plants and large predatory fish, P. Natl. Acad. Sci. USA, 107, 17911-17915, https://doi.org/10.1073/pnas.1011287107, 2010.

Dahl, T. W., Harding, M. A., Brugger, J., Feulner, G., Norrman, K., Lomax, B. H., and Junium, C. K.: Low atmospheric CO₂ levels before the rise of forested ecosystems, Nat. Commun., 13, 7616, https://doi.org/10.1038/s41467-022-35085-9, 2022.

Da Silva, A.-C., Sinnesael, M., Claeys, P., Davies, J. H. F. L., De Winter, N. J., Percival, L. M. E., Schaltegger, U., and De Vleeschouwer, D.: Anchoring the Late Devonian mass extinction in absolute time by integrating climatic controls and radio-isotopic dating, Sci. Rep., 10, 12940, https://doi.org/10.1038/s41598-020-69097-6, 2020.

De Vleeschouwer, D., Rakociński, M., Racki, G., Bond, D. P., Sobień, K., and Claeys, P.: The astronomical rhythm of Late Devonian climate change (Kowala section, Holy Cross Mountains, Poland), Earth Planet. Sc. Lett., 365, 25-37, https://doi.org/10.1016/j.epsl.2013.01.016, 2013.

De Vleeschouwer, D., Crucifix, M., Bounceur, N., and Claeys, P.: The impact of astronomical forcing on the Late Devonian greenhouse climate, Global Planet. Change, 120, 65-80, https://doi.org/10.1016/j.gloplacha.2014.06.002, 2014.

De Vleeschouwer, D., Boulvain, F., Da Silva, A.-C., Pas, D., Labaye, C., and Claeys, P.: The astronomical calibration of the Givetian (Middle Devonian) timescale (Dinant Synclinorium, Belgium), Geological Society, London, Special Publications, 414, 245-256, https://doi.org/10.1144/SP414.3, 2015.

De Vleeschouwer, D., Da Silva, A.-C., Sinnesael, M., Chen, D., Day, J. E., Whalen, M. T., Guo, Z., and Claeys, P.: Timing and pacing of the Late Devonian mass extinction event regulated by eccentricity and obliquity, Nat. Commun., 8, 2268, https://doi.org/10.1038/s41467-017-02407-1, 2017.

Donnadieu, Y., Pucéat, E., Moiroud, M., Guillocheau, F., and De-coninck, J.-F.: A better-ventilated ocean triggered by Late Creaceous changes in continental configuration, Nat. Commun., 7, 10316, https://doi.org/10.1038/ncomms10316, 2016.

Dopieralska, J.: Reconstructing seawater circulation on the Moroccan shelf of Gondwana during the Late Devonian: Evidence from Nd isotope composition of conodonts, Geochem. Geophy. Geosy., 10, 2008GC002247, https://doi.org/10.1029/2008GC002247, 2009.

Eberhard, J., Bevan, O. E., Feulner, G., Petri, S., van Hunen, J., and Baldini, J. U. L.: Sensitivity of Neoproterozoic snowball-Earth inceptions to continental configuration, orbital geometry, and volcanism, Clim. Past, 19, 2203-2235, https://doi.org/10.5194/cp-19-2203-2235, 2023.

Edwards, N. R. and Marsh, R.: Uncertainties due to transport-parameter sensitivity in an efficient 3-D ocean-climate model, Clim. Dynam., 24, 415-433, https://doi.org/10.1007/s00382-004-0508-8, 2005.

Ferreira, D., Marshall, J., and Rose, B.: Climate Determinism Re-visited: Multiple Equilibria in a Complex Climate Model, J. Climate, 24, 992-1012, https://doi.org/10.1175/2010JCLI3580.1, 2011.

Foster, G. L., Royer, D. L., and Lunt, D. J.: Future climate forcing potentially without precedent in the last 420 million years, Nat. Commun., 8, 14845, https://doi.org/10.1038/ncomms14845, 2017.

Gérard, J. and Crucifix, M.: Diagnosing the causes of AMOC slow-down in a coupled model: a cautionary tale, Earth Syst. Dynam., 15, 293-306, https://doi.org/10.5194/esd-15-293-2024, 2024.

Gough, D.: Solar interior structure and luminosity variations, in: Physics of Solar Variations: Proceedings of the 14th ESLAB Symposium held in Scheveningen, the Netherlands, 16-19 September 1980, Springer, 21-34, https://doi.org/10.1007/978-94-010-9633-1_4, 1981.

Haddad, E. E., Boyer, D. L., Droser, M. L., Lee, B. K., Lyons, T. W., and Love, G. D.: Ichnofabrics and chemostratigraphy argue against persistent anoxia during the Upper Kellwasser Event in New York State, Palaeogeogr. Palaeoclim., 490, 178-190, https://doi.org/10.1016/j.palaeo.2017.10.025, 2018.

Hedhli, M., Grasby, S., Henderson, C., and Davis, B.: Multiple diachronous “Black Seas” mimic global ocean anoxia during the latest Devonian, Geology, 51, 973-977, https://doi.org/10.1130/G51394.1, 2023.

Hibler, W. D.: A dynamic thermodynamic sea ice model, J. Phys. Oceanogr., 9, 815-846, https://doi.org/10.1175/1520-0485(1979)009<0815:ADTSIM>2.0.CO;2, 1979.

Hughes, T. M. C. and Weaver, A. J.: Multiple Equilibrium of an Asymmetric Two-Basin Ocean Model, J. Phys. Oceanogr., 24, 619-637, https://doi.org/10.1175/1520-0485(1994)024<0619:MEOAAT>2.0.CO;2, 1994.

Hülse, D., Lau, K. V., Van De Velde, S. J., Arndt, S., Meyer, M., and Ridgwell, A.: End-Permian marine extinction due to temperature-driven nutrient recycling and euxinia, Nat. Geosci., 14, 862-867, https://doi.org/10.1038/s41561-021-00829-7, 2021.

Jenkyns, H. C.: Geochemistry of oceanic anoxic events, Geochem. Geophy. Geosy., 11, 2009GC002788, https://doi.org/10.1029/2009GC002788, 2010.

Joachimski, M., Breisig, S., Buggisch, W., Talent, J., Mawson, R., Gereke, M., Morrow, J., Day, J., and Weddige, K.: Devonian climate and reef evolution: Insights from oxygen isotopes in apatite, Earth Planet. Sc. Lett., 284, 599-609, https://doi.org/10.1016/j.epsl.2009.05.028, 2009.

Kabanov, P., Hauck, T. E., Gouwy, S. A., Grasby, S. E., and Van Der Boon, A.: Oceanic anoxic events, photic-zone euxinia, and controversy of sea-level fluctuations during the Middle-Late Devonian, Earth-Sci. Rev., 241, 10445, https://doi.org/10.1016/j.earscirev.2023.10445, 2023.

Kaiser, S. I., Aretz, M., and Becker, R. T.: The global Hangenberg Crisis (Devonian-Carboniferous transition): review of a first-order mass extinction, Geological Society, London, Special Publications, 423, 387-437, https://doi.org/10.1144/SP423.9, 2016.

Laskar, J., Fienga, A., Gastineau, M., and Manche, H.: La2010: a new orbital solution for the long-term motion of the Earth, Astron.
Astrophys., 532, A89, https://doi.org/10.1051/0004-6361/201116836, 2011.

Lawver, L. A. and Gahagan, L. M.: Evolution of Cenozoic seaways in the circum-Antarctic region, Palaeogeogr. Palaeocl., 198, 11–37, https://doi.org/10.1016/S0031-0182(03)00392-4, 2003.

Li, Y. and Yang, H.: A Theory for Self-Sustained Multicentennial Oscillation of the Atlantic Meridional Overturning Circulation, J. Climate, 35, 5883–5896, https://doi.org/10.1175/JCLI-D-21-0685.1, 2022.

Lohmann, J., Dijkstra, H. A., Jochum, M., Lucarini, V., and Ditlevsen, P. D.: Multistability and intermediate tipping of the Atlantic Ocean circulation, Science Advances, 10, eadi4253, https://doi.org/10.1126/sciadv.adi4253, 2024.

Lunt, D. J., Williamson, M. S., Valdes, P. J., Lenton, T. M., and Marsh, R.: Comparing transient, accelerated, and equilibrium simulations of the last 30 000 years with the GENIE-1 model, Clim. Past, 2, 221–235, https://doi.org/10.5194/cp-2-221-2006, 2006.

Mack, C. L. J.: Modelling the Late Devonian Climate: Implications for evolution and extinction, PhD thesis, University of Oxford, https://ora.ox.ac.uk/objects/uuid:95ac6924-cf4b-43bb-9c02-fb9d37acb596 (last access: 20 January 2025), 2021.

Maier-Reimer, E.: Geochemical cycles in an ocean general circulation model. Preindustrial tracer distributions, Global Bio-geochem. Cy., 7, 645–677, https://doi.org/10.1029/93GB01355, 1993.

Marcilly, C. M.-, Maffre, P., Le Hir, G., Pohl, A., Fluteau, F., Godéris, Y., Donnadieu, Y., H. Heimdal, T., and Torsvik, T. H.: Understanding the early Paleozoic carbon cycle balance and climate change from modelling, Earth Planet. Sc. Lett., 594, 117717, https://doi.org/10.1016/j.epsl.2022.117717, 2022.

Margazoglou, G., Grafke, T., Laio, A., and Lucarini, V.: Dynamical landscape and multistability of a climate model, P. R. Soc. A, 477, 20210019, https://doi.org/10.1098/rspa.2021.0019, 2021.

Marsh, R., Müller, S. A., Yool, A., and Edwards, N. R.: Incorporation of the C-GOLDSTEIN efficient climate model into the GENIE framework: “eb\_go\_gs” configurations of GENIE, Geosci. Model Dev., 4, 957–992, https://doi.org/10.5194/gmd-4-957-2011, 2011.

Marsh, R., Sóbester, A., Hart, E. E., Oliver, K. I. C., Edwards, N. R., and Cox, S. J.: An optimally tuned ensemble of the “eb\_go\_gs” configuration of GENIE: parameter sensitivity and bifurcations in the Atlantic overturning circulation, Geosci. Model Dev., 6, 1729–1744, https://doi.org/10.5194/gmd-6-1729-2013, 2013.

Merdith, A. S., Williams, S. E., Collins, A. S., Tetley, M. G., Mulder, J. A., Blades, M. L., Young, A., Armistead, S. E., Cannon, J., Zahirovic, S., and Müller, R. D.: Extending full-plate tectonic models into deep time: Linking the Neoproterozoic and the Phanerozoic, Earth-Sci. Rev., 214, 103477, https://doi.org/10.1016/j.earscirev.2020.103477, 2021.

Meyers, S. R., Sageman, B. B., and Arthur, M. A.: Obliquity forcing of organic matter accumulation during Oceanic Anoxic Event 2, Paleoceanography, 27, 2012PA002286, https://doi.org/10.1029/2012PA002286, 2012.

Monteiro, F. M., Pancost, R. D., Ridgwell, A., and Donnadieu, Y.: Nutrients as the dominant control on the spread of anoxia and euxinia across the Cenomanian-Turonian oceanic anoxic event (OAE2): Model-data comparison, Paleoceanography, 27, 2012PA002351, https://doi.org/10.1029/2012PA002351, 2012.

Paschall, O., Carmichael, S. K., Königshof, P., Waters, J. A., Ta, P. H., Komatsu, T., and Dombrowski, A.: The Devonian-Carboniferous boundary in Vietnam: Suspended ocean anoxia with a volcanic trigger for the Hangenberg Crisis?, Global Planet. Change, 175, 64–81, https://doi.org/10.1016/j.gloplacha.2019.01.021, 2019.

Peltier, W. R. and Vettoretti, G.: Dansgaard-Oeschger oscillations predicted in a comprehensive model of glacial climate: A “kicked” salt oscillator in the Atlantic, Geophys. Res. Lett., 41, 7306–7313, https://doi.org/10.1002/2014GL061413, 2014.

Percival, L., Bond, D., Rakociński, M., Marynowski, L., Hood, A., Adatte, T., Spangenberg, J., and Föllmi, K.: Phosphorus-cycle disturbances during the Late Devonian anoxic events, Global Planet. Change, 184, 103070, https://doi.org/10.1016/j.gloplacha.2019.103070, 2020.

Percival, L. M., Marynowski, L., Baudin, F., Goderis, S., De Vleeschouwer, D., Rakociński, M., Narkiewicz, K., Corradini, C., Da Silva, A.-C., and Claeys, P.: Combined nitrogen-isotope and cyclostratigraphy evidence for temporal and spatial variability in Frasnian–Famennian Environmen tal Change, Geochem. Geophy. Geosy., 23, e2021GC010308, https://doi.org/10.1029/2021GC010308, 2022.

Pohl, A., Lu, Z., Lu, W., Stockey, R. G., Elrick, M., Li, M., Desrochers, A., Shen, Y., He, R., Finnegan, S., and Ridgwell, A.: Vertical decoupling in Late Ordovician anoxia due to re-organization of ocean circulation, Nat. Geosci., 14, 868–873, https://doi.org/10.1038/s41561-021-00843-9, 2021.

Pohl, A., Ridgwell, A., Stockey, R. G., Thomazo, C., Keane, A., Vennin, E., and Scotese, C. R.: Continental configuration controls ocean oxygenation during the Phanerozoic, Nature, 608, 523–527, https://doi.org/10.1038/s41586-022-05018-z, 2022.

Pohl, A., Stockey, R. G., Dai, X., Yohler, R., Le Hir, G., Hülse, D., Brayard, A., Finnegan, S., and Ridgwell, A.: Why the Early Paleozoic was intrinsically prone to marine extinction, Science Advances, 9, eadg7679, https://doi.org/10.1126/sciadv.adg7679, 2023.

Rae, J. W. B., Gray, W. R., Wills, R. C. J., Eisenman, I., Fitzhugh, B., Fotheringham, M., Littley, E. F. M., Rafter, P. A., Rees-Owen, R., Ridgwell, A., Taylor, B., and Burke, A.: Overturning circulation, nutrient limitation, and warming in the Glacial North Pacific, Science Advances, 6, eabd1654, https://doi.org/10.1126/sciadv.abd1654, 2020.

Rahmstorf, S.: Multiple Convection Patterns and Ther-mohaline Flow in an Idealized OGCM, J. Climate, 8, 3028–3039, https://doi.org/10.1175/1520-0442(1995)008<3028:MCPATF>2.0.CO;2, 1995.

Rahmstorf, S., Crucifix, M., Ganopolski, A., Goosse, H., Ka-menkovich, I., Knutti, R., Lohmann, G., Marsh, R., Mysak, A., Wang, Z., and Weaver, A. J.: Thermohaline circulation hysteresis: A model intercomparison, Geophys. Res. Lett., 32, 2005GL023655, https://doi.org/10.1029/2005GL023655, 2005.

Reinhard, C. T., Olson, S. L., Kirtland Turner, S., Pälike, C., Kanzaki, Y., and Ridgwell, A.: Oceanic and atmospheric methane cycling in the cGENIE Earth system model – release v0.9.14, Geosci. Model Dev., 13, 5687–5706, https://doi.org/10.5194/gmd-13-5687-2020, 2020.

Ridgwell, A.: derpycode/muffingen: v0.9.24 (v0.9.24), Zenodo [code], https://doi.org/10.5281/zenodo.7545809, 2023.

Ridgwell, A., Hargreaves, J. C., Edwards, N. R., Annan, J. D., Lenton, T. M., Marsh, R., Yool, A., and Watson, A.: Marine geochemical data assimilation in an efficient Earth System Model of global biogeochemical cycling, Biogeosciences, 4, 87–104, https://doi.org/
10.5194/bg-4-87-2007, 2007.

Ridgwell, A., Hülse, D., Peterson, C., Ward, B., Ted, evansmn, and Jones, R.: derpycode/muffindoc: v0.24-574-NSF (v0.24-574- NSF), Zenodo [code], https://doi.org/10.5281/zenodo.13377225, 2024.

Ridgwell, A., Reinhard, C., van de Velde, S., Adloff, M., Mon- teiro, F., Camillalxy98, Vervoort, P., Kanzaki, Y., Ward, B., Hülse, D., Wilson, J., Li, M., InkyANB, and Kirtland Turner, S.: derpycode/cgenie.muffin: v0.9.59 (v0.9.59), Zenodo [code], https://doi.org/10.5281/ zenodo.14680971, 2025.

Sakai, K. and Peltier, W. R.: Dansgaard–Oeschger Oscilla- tions in a Coupled Atmosphere–Ocean Climate Model, J. Climate, 10, 949–970, https://doi.org/10.1175/1520- 0442(1997)010<0949:DOOIAC>2.0.CO;2, 1997.

Sarr, A., Donnadieu, Y., Laugié, M., Ladant, J., Suchéras- Marx, B., and Raisson, F.: Ventilation Changes Drive Orbital-Scale Deoxygenation Trends in the Late Creta- ceous Ocean, Geophys. Res. Lett., 49, e2022GL099830, https://doi.org/ 10.1029/2022GL099830, 2022.

Scotese, C. R. and Wright, N.: PALEOMAP paleodig- ital elevation models (PaleoDEMS) for the Phanero- zoic, Paleomap Proj, https://www.earthbyte.org/ paleodem-resource-scotese-and-wright-2018/ (last access: 20 January 2025), 2018.

Semtner Jr., A. J.: A model for the thermodynamic growth of sea ice in numerical investigations of climate, J. Phys. Oceanogr., 6, 379–389, https://doi.org/10.1175/1520- 0485(1976)006<0379:AMFTTG>2.0.CO;2, 1976.

Sévellec, F. and Fedorov, A. V.: Millennial Variability in an Ide- alized Ocean Model: Predicting the AMOC Regime Shifts, J. Climate, 27, 3551–3564, https://doi.org/10.1175/JCLI-D-13- 00450.1, 2014.

Song, H., Song, H., Algeo, T. J., Tong, J., Romaniello, S. J., Zhu, Y., Chu, D., Gong, Y., and Anbar, A. D.: Uranium and carbon isotopes document global-ocean redox-productivity relationships linked to cooling during the Frasnian-Famennian mass extinc- tion, Geology, 45, 887–890, https://doi.org/10.1130/G39393.1, 2017.

Stein, W. E., Berry, C. M., Morris, J. L., Hernick, L. V., Man- nolini, F., Ver Straeten, C., Landing, E., Marshall, J. E., Wellman, C. H., Beerling, D. J., and Leake, J. R.: Mid-

Devonian Archaeopteris Roots Signal Revolutionary Change in Earliest Fossil Forests, Curr. Biol., 30, 421–431.e2, https://doi.org/ 10.1016/j.cub.2019.11.067, 2020.

Thornalley, D. J. R., Elderfield, H., and McCave, I. N.: Holocene oscillations in temperature and salinity of the surface subpolar North Atlantic, Nature, 457, 711–714, https://doi.org/10.1038/nature07717, 2009.

Turner, B. W. and Slatt, R. M.: Assessing bottom water anoxia within the Late Devonian Woodford Shale in the Arkoma Basin, southern Oklahoma, Mar. Petrol. Geol., 78, 536–546, https://doi.org/10.1016/j.marpetgeo.2016.10.009, 2016.

Valdes, P. J., Scotese, C. R., and Lunt, D. J.: Deep ocean temperatures through time, Clim. Past, 17, 1483–1506, https://doi.org/10.5194/ cp-17-1483-2021, 2021.

Vandenbroucke, T. R. A., Emsbo, P., Munnecke, A., Nuns, N., Duponchel, L., Lepot, K., Quijada, M., Paris, F., Servais, T., and Kiessling, W.: Metal-induced malformations in early Palaeozoic plankton are harbingers of mass extinction, Nat. Commun., 6, 7966, https://doi.org/ 10.1038/ncomms8966, 2015.

Vérard, C.: Panalesis: towards global synthetic palaeogeographies using integration and coupling of manifold models, Geol. Mag., 156, 320–330, https://doi.org/10.1017/S0016756817001042, 2019a.

Vérard, C.: Plate tectonic modelling: review and perspectives, Geol. Mag., 156, 208–241, https://doi.org/10.1017/ S0016756817001030, 2019b.

Vérard, C., Hochard, C., Baumgartner, P. O., Stampfli, G. M., and Liu, M.: 3D palaeogeographic recon- structions of the Phanerozoic versus sea-level and Sr- ratio variations, Journal of Palaeogeography, 4, 64–84, https://doi.org/10.3724/SP.J.1261.2015.00068, 2015.

Vervoort, P., Kirtland Turner, S., Rochholz, F., and Ridgwell, A.: 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, Paleooceanography and Paleoclimatology, 39, e2023PA004826, https://doi.org/10.1029/2023PA004826, 2024.

Vettoretti, G., Ditlevsen, P., Jochum, M., and Rasmussen, S. O.: Atmospheric CO₂ control of spontaneous millennial- scale ice age climate oscillations, Nat. Geosci., 15, 300–306, https://doi.org/10.1038/s41561-022-00920-7, 2022.

Von Der Heydt, A. and Dijkstra, H. A.: Effect of ocean gate- ways on the global ocean circulation in the late Oligocene and early Miocene, Paleooceanography, 21, 2005PA001149, https://doi.org/10.1029/2005PA001149, 2006.

Ward, B. A., Wilson, J. D., Death, R. M., Monteiro, F. M., Yool, A., and Ridgwell, A.: EcoGenIE 1.0: plankton ecology in the cGE- nIE Earth system model, Geosci. Model Dev., 11, 4241–4267, https://doi.org/10.5194/gmd-11-4241-2018, 2018.

Weaver, A. J., Eby, M., Wiebe, E. C., Bitz, C. M., Duffy, P., B., Ewen, T. L., Fanning, A. F., Holland, M. M., MacFadyen, A., Matthews, H. D., Meissner, K. J., Saenko, O., Schmittner, A., Wang, H., and Yoshimori, M.: The UVic earth system cli- mate model: Model description, climatology, and applications to past, present and future climates, Atmos. Ocean, 39, 361–428, https://doi.org/10.1080/07055900.2001.9649686, 2001.

Weijer, W., Cheng, W., Drijfhout, S. S., Fedorov, A. V., Hu, A., Jackson, L. C., Liu, W., McDonagh, E. L., Mecking, J. V., and Zhang, J.: Stability of the Atlantic Meridional Overturning Cir- culation: A Review and Synthesis, J. Geophys. Res.-Oceans, 124, 5336–5375, https://doi.org/10.1029/2019jc015083, 2019.

Wichern, N. M. A., Bialik, O. M., Nohl, T., Percival, L. M. E., Becker, R. T., Kaskes, P., Claeys, P., and De Vleeschouwer, D.: Astronomically paced climate and carbon cycle feedbacks in the lead-up to the Late Devonian Kellwasser Crisis, Clim. Past, 20, 415– 448, https://doi.org/10.5194/cp-20-415-2024, 2024.

Wong Hearing, T. W., Pohl, A., Williams, M., Donnadieu, Y., Harvey, T. H. P., Scotese, C. R., Sepulchre, P., Franc, A., and Vandenbroucke, T. R. A.: Quantitative comparison of geological data and model simulations constrains early Cambrian geography and climate, Nat. Commun., 12, 3868, https://doi.org/10.1038/s41467-021-24141-5, 2021.

Yang, S., Galbraith, E., and Palter, J.: Coupled climate impacts of the Drake Passage and the Panama Seaway, Clim. Dynam., 43, 37–52,