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

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 PO₄ concentration and atmospheric pO₂), climatic forcing (pCO₂), 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 PO₄ concentration, atmospheric pO₂ and pCO₂, 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 pO₂ and increased oceanic PO₄ concentration exacerbate ocean anoxia, consistent with established knowledge. The variation of pCO₂ 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 changes in the distribution and amount of incoming solar radiation received by the Earth through the modification of its position and orientation in space. Some of the Devonian OAEs have been linked to astronomical parameters through recent findings in different 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 Hangenberg 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, astronomical forcing may not have been enough on its own and would most probably have acted as a perturbation driving a preconditioned 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 factors preconditioning the ocean to anoxia, which is crucial for unravelling the mechanisms behind Devonian OAEs. A specific challenge for developing plausible anoxic scenarios involving ocean dynamics during the Devonian primarily comes from uncertainties associated with various factors, including atmospheric CO₂ concentration, orbital configuration, 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 analyses around those parameters. First, we explore the contribution of Devonian continental configuration and bathymetry to the global mean oceanic oxygen content through modification of ocean circulation patterns (Sect. 3.1). Of particular concern are the uncertainties surrounding palaeogeography, especially when it comes to bathymetry, as they exert 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 (lithosphere) has been destroyed in subduction zones due to plate tectonics activity (Scotese and Wright, 2018), preventing direct 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 dynamics 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 reconstructions: 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 palaeogeographic 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 continental configurations (Sect. 3.2). We perform sensitivity analyses of three biogeochemical quantities: atmospheric pCO₂ (partial pressure of CO₂), global mean oceanic PO₄ concentration ([PO₄]), and atmospheric pO₂ (partial pressure of O₂). 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 allow us to explore the physical contribution of essential components 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 (Crichton et al., 2021), to explore the contribution of ocean circulation reorganization to Late Ordovician anoxia (Pohl et al., 2021), and to understand the causes of the Atlantic Meridional Overturning Circulation slowdown under global warming (Gérard and Crucifix, 2024).

2 Methods

2.1 Model description

cGENIE is an Earth system model of intermediate complexity 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 biogeochemical components. The physical part of cGENIE consists 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 mixing. The framework relies on parameterized isoneutral diffusion and eddy-induced advection to efficiently transport heat, salinity, and biogeochemical tracers. The ocean surface 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 temperature and specific humidity to compute the heat and moisture 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 dynamic hydrological scheme for continents, setting land evaporation 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 variables 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 directly 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 Ridgwell et al. (2007), manages air–sea gas exchange along with the transformations and spatial redistribution of biogeochemical compounds at both surface and interior levels. Biogenically induced chemical fluxes are primarily constrained by nutrient limitation processes (Maier-Reimer, 1993; Ward et al., 2018). Concentration storage of atmospheric composition 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 sediment 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 interactions between physical and biogeochemical processes. To ensure simplicity and computer efficiency, the configuration of cGENIE used here includes the C-GOLDSTEIN framework together with the BIOGEM and ATCHEM biogeochemical modules.

2.2 Palaeogeographic reconstructions

We consider the continental reconstructions produced by Scotese and Wright (2018) and Vérard (2019a). The palaeogeographic reconstructions of Scotese and Wright (2018) span the Devonian from 360 to 420 Ma at roughly 5 Myr intervals, encompassing a total of 13 distinct configurations. They are inferred from a local determination of the palaeoenvironment primarily using lithostratigraphic and palaeontological analyses, with a focus on shoreline regions. Local reconstructions are then extrapolated to address temporal uncertainties and geological gaps, ultimately contributing to a global palaeogeographic reconstruction through interpolation. This process is repeated across multiple sites worldwide to create a detailed understanding of past environmental conditions. The palaeogeographic reconstructions of Vérard (2019a), on the other hand, are obtained through the conversion of a full-plate tectonic model (global coverage with fully closed plate boundaries; see Vérard, 2019b) into a quantified topography (i.e. both orography and bathymetry) following Vérard et al. (2015). With the preliminary version of Panalesis (Panalesis.v0, Vérard, 2019a), these reconstructions are available at 370, 383, 393, 408, and 420 Ma. The approach of Scotese and Wright (2018) has spatial gaps, especially in the ocean because there is no constraining information. 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 reconstructions, 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 reconstructions 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 series of reconstructions side by side, we adjusted all Vérard (2019a) reconstructions to conform to the sea level specifications 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 feature 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 reason, consider a scheme of sinking organic matter remineralization independent of temperature. It allows us to better isolate the contribution of physical changes to the biogeochemical 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 dependence on the land–sea mask. Although variations in continental configurations are known, in reality, to affect albedo, our experiments are constrained by the absence of a land surface scheme (i.e. no snow cover, clouds, or ice sheets), limiting the realism, similarly to Cermeño et al. (2022) and Vervoort 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 indicates 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 Devonian 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).

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.

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.

2.3 Description of the experiments

We present three sets of experiments where all simulations were run for 10 000 years until the system reached equilibrium 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₂ (280 ppm) and pO₂ (20.95 %) as well as global mean oceanic [PO₄] (2.159 µmol kg⁻¹). Throughout this study, these preindustrial values are also referred to as 1×pCO₂, 1×pO₂, and 1×[PO₄], respectively. Additional simulations with the same solar constant across all 13 Scotese and Wright (2018) continental 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 continental 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 atmospheric pCO₂, with 1, 2, 4, 8, and 16 times the preindustrial 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 configurations. Values for pCO₂ and pO₂ were selected to capture a significant portion of the range and uncertainties associated 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 parameter 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 (removing any contribution from the precession) and obliquity of 23°. Table 1 provides a schematic overview of the three experiment sets.

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 combinations have been tested.

3 Results

3.1 Set 1: continental configuration

The impact of palaeogeographical evolution on oceanic oxygenation and ventilation using the continental reconstructions of Scotese and Wright (2018) (red dots) and Vérard (2019a) (blue dots) is shown in Fig. 3. The two sets of reconstructions generate fairly similar mean oxygen concentrations ([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 deepwater 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\text{ value}=\mathrm{4.3}\times {\mathrm{10}}^{-\mathrm{7}}$ and r²=0.93 with p value=0.008 for the set of reconstructions from Scotese and Wright, 2018, and Vérard, 2019a, respectively). This suggests that the variations in oxygen concentration are primarily explained by changes in ventilation, caused by modifications in ocean circulation dynamics induced by the prescribed differences in continental configuration. Second-order variations in oxygen observed in our simulations can be attributed to changes in oxygen solubility due to alterations in the solar constant and modifications in exported organic matter (biological pump), which are influenced by temperature adjustments and differences in continental configuration. All the simulations converge towards a stable steady state, except Scot390M. The latter is characterized 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 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⁻¹.

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.

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évellec 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, 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 reconstructions 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) reconstructions 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 exhibit 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).

Figure 4 Anoxic ([O₂]<6.5 µmol kg⁻¹) and hypoxic ([O₂]<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.

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, following 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 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 advection 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.

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.

3.2 Set 2: geochemical state and climate forcing

3.2.1 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₂, global mean oceanic [PO₄], and atmospheric pO₂. Here, atmospheric pCO₂ is used as a climate forcing, where a doubling of pCO₂ 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 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₂ 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.

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 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) Scotese 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).

3.2.2 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 doubling 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 reconstructions 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).

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

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 explain the small mismatch between Vera370M and Vera420M (blue and green curves in Fig. 7b) when increasing from 2 to 4×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.

3.2.3 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 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.

Figure 9 Representation of the limit cycles in (a) benthic [O₂] (µmol kg⁻¹), benthic temperature (°C), (b) maximum overturning strength (Sv), and northern sea ice extent (%) for different continental configurations and atmospheric pCO₂. 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₂.

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 that the anoxic extent obtained with 8×pCO₂ could vary up to more than 20 % against less than 10 % with 16×pCO₂. 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×pCO₂ (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 pCO₂ 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×pCO₂, the warming is linked to an increase in the benthic anoxic extent, whereas going from 8 to 16×pCO₂ 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 pCO₂ increases, the northern circulation cell intensifies, gradually taking over the southern cell because of the seesaw effect (see Fig. 9d–f). Although the northern overturning cell strengthens with each increase in pCO₂, 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×pCO₂ do deep waters of northern origin efficiently propagate into the benthic ocean, leading to increased oxygen levels and the reduction of benthic anoxic extent.

3.2.4 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 behaviours, potential signs of bistability, and their implications for the ocean time response. Bistability arises when a system, 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 increasing (blue) and decreasing (red) pCO₂.

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.

Given that the simulation restarts from a stable state (previous spin-up of 10 000 years) and considering the slow evolution 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 oscillations, 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 discrepancies 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 surrounding 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, indicating that it only acts as a boundary when the limit cycle is already present. This feature of the oscillation interval accounts 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.

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 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₂]<62.5 µmol kg⁻¹. 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 remove any contribution from precession and 0.06, respectively. In HE simulations, precession effects are maximized, with the different months indicating when in the year the perihelion is reached. All simulations were performed using 2×pCO₂, 1.25×[PO₄], and 0.7×pO₂.

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 tend to disappear, although the threshold depends on the precise 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.

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₂.

4 Discussion

4.1 Implications for current hypotheses about Devonian anoxic events

We used cGENIE to investigate the contribution of continental configuration, geochemical states, and astronomical forcing to Devonian anoxia. First, we found that the continental 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₂, oceanic [PO₄], pO₂, and astronomical forcings. The differences between reconstructions from Scotese and Wright (2018) and Vérard (2019a) underscore the importance of bathymetry for circulation dynamics and anoxic extent. Notably, the presence of self-sustained oscillations in the Scot390M configuration emphasizes the nonlinear nature 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 reliability of different palaeogeographic reconstructions in the deep geological past. Our findings also reveal that continental 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 experimental design, the continental configuration offers an explanation 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₂] 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 emphasized the potentially important role of palaeogeography as a preconditioning factor to the development of OAEs. Pohl et al. (2022) showed, with cGENIE, how continental rearrangements, even small, during the Phanerozoic eon led to profound variations in ocean oxygenation through reorganizations of large-scale ocean circulation. Recent studies on nitrogen 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 factors. They suggest that the Hangenberg black shales do not represent a single global anoxic event but rather multiple, diachronous basins similar to the Black Sea, occurring around the world. To explain this, they propose that the tectonic configuration 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₂ and global mean oceanic [PO₄] on the benthic anoxic extent. As a general rule, both an increase in [PO₄] and a decrease in atmospheric pO₂ 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 terrestrial phosphorus from enhanced weathering. Our analysis also reveals that continental configuration significantly influences the amplitude and even sometimes the sign of the contribution of a pCO₂ increase to benthic oxygen levels. This is due to the relative importance of two competing mechanisms: surface solubility loss and overturning strengthening. Studies generally suggest that OAEs of the Devonian occurred under warm conditions; however, some research indicates that these events could have also taken place during periods of a relatively cold global climate (Percival et al., 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 dynamics with substantial effects on ocean oxygen levels including self-sustained oscillations, variations in overturning cell dominance, modification of sea ice cover, emergence of multistability, 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 relevant 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 prolonged duration of certain OAEs. Applying constraints on atmospheric pCO₂ and surface temperature, such as those presented 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 extensive oxygen-depleted zones, with eccentricity and precession 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 astronomical forcing on ocean anoxia. Additionally, astronomical forcing can shift the system's equilibrium, transitioning between stable oscillations and steady states, with obliquity and eccentricity being the primary determinants of these equilibrium states.

4.2 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. Reassuringly, cGENIE has repeatedly demonstrated its ability to capture large-scale ocean circulation patterns (Rahmstorf et al., 2005; Lunt et al., 2006; Cao et al., 2009; Marsh et al., 2013; Pohl et al., 2022; Gérard and Crucifix, 2024), providing confidence in the robustness of the simulated global circulation. Furthermore, even with its coarse spatial resolution, the model effectively captured the variations between continental reconstructions to produce very distinct ocean dynamics regimes. Yet, model limitations have some consequences regarding 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 active 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 detailed 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 considerations from Weijer et al. (2019), multistability and limit cycles 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 include 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, potentially impacting the circulation. Finally, cGENIE's resolution is too coarse to properly capture the connection between the deep ocean and epicontinental margins. Yet, the development of deep-ocean anoxia would help explain the synchronous character of anoxia observed in distant geological 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 simulations 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 reconstructions, which would combine continental locations and shallow waters from Scotese and Wright (2018) with the deep ocean from Vérard (2019a). Such hybrid configurations were explored by Cermeño et al. (2022), merging 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, particularly EMIC and Coupled Model Intercomparison Project-class models to elucidate similarities and potential disparities between cGENIE and finer models.

5 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 systematic 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.

1. Continental configuration consistently modulates the signal of all other considered components of the Devonian climate (oceanic [PO₄], atmospheric pO₂, and pCO₂ as well as orbital forcing) on anoxia, highlighting its crucial contribution to ocean anoxic events. Furthermore, the continental configuration greatly influences ocean circulation and global oxygen content, with, at one point, the development of stable oceanic oscillations 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 during the Middle and Late Devonian periods by creating more extensive oxygen-depleted zones. Our investigation also revealed the presence of pronounced multistability (more than one stable state of the system exists for the same climatic configuration) for a specific continental configuration, providing a potential mechanism to account for the prolonged duration of these events (“lockdown” of the system in a state with strong oxygen-depleted conditions).

2. The influence of global mean oceanic [PO₄] and atmospheric [O₂] on ocean oxygen content is in line with current knowledge. Simulations indicate that both decreased atmospheric [O₂] and increased oceanic [PO₄] exacerbate ocean anoxia. The variation of atmospheric pCO₂ (here used as a climate forcing to modify temperatures) revealed important changes in ocean anoxia through modification of oxygen solubility, ocean circulation, 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 sensitive to continental configuration. The hysteresis feature also provides a plausible mechanism to explain the prolonged duration of some ocean anoxic events. Furthermore, these processes and behaviours lead to important modifications of the ocean oxygen levels, therefore strongly impacting the preconditioning of the ocean to anoxia.

3. Astronomical forcing substantially impacts ocean anoxia by altering ocean circulation and oxygen solubility. Obliquity consistently emerges as the primary orbital parameter driving the vast majority of observed oceanic oxygen variations due to its strong annual average influence on high latitudes. In our experiments, high obliquity persistently led to more extensive oxygen-depleted zones. Simulations demonstrate that astronomical 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 cycles, 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, influenced by continental configuration, CO₂ levels, and orbital forcing, exert significant impacts on ocean anoxia. Ocean circulation revealed surprisingly complex dynamics, with response times ranging from years to tens of thousands of years, which are relevant for understanding anoxia in the Devonian 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 cycle observed in the Scot420M continental setup. Similar to Scot390M, we find that these oscillations arise from the interplay 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 composite of a similar limit cycle as in Scot390M and a deep-decoupling oscillation. This additional deep decoupling oscillation involves the gradual accumulation of heat below the surface, as shown in Fig. A1. Mack (2021) have also observed analogous deep decoupling oscillations, manifesting as two quasi-stable states – an unventilated phase and a ventilated phase – using the MITgcm ocean model for Devonian palaeoclimate conditions. Furthermore, the close similarity between Fig. A1 and Fig. 3.8 from Mack (2021) further demonstrates the ability of cGENIE to capture complex large-scale ocean dynamics patterns. The combination of these two limit cycles generates this “two-step” pattern of the oscillations present for the Scot420M continental configuration (see Fig. 9a and b).

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₂.