Continental configuration controls ocean oxygenation during the Phanerozoic

The early evolutionary and much of the extinction history of marine animals is thought to be driven by changes in dissolved oxygen concentrations ([O2]) in the ocean1–3. In turn, [O2] is widely assumed to be dominated by the geological history of atmospheric oxygen (pO2)4,5. Here, by contrast, we show by means of a series of Earth system model experiments how continental rearrangement during the Phanerozoic Eon drives profound variations in ocean oxygenation and induces a fundamental decoupling in time between upper-ocean and benthic [O2]. We further identify the presence of state transitions in the global ocean circulation, which lead to extensive deep-ocean anoxia developing in the early Phanerozoic even under modern pO2. Our finding that ocean oxygenation oscillates over stable thousand-year (kyr) periods also provides a causal mechanism that might explain elevated rates of metazoan radiation and extinction during the early Palaeozoic Era6. The absence, in our modelling, of any simple correlation between global climate and ocean ventilation, and the occurrence of profound variations in ocean oxygenation independent of atmospheric pO2, presents a challenge to the interpretation of marine redox proxies, but also points to a hitherto unrecognized role for continental configuration in the evolution of the biosphere. Analysis of a series of Earth system model experiments shows that continental rearrangement during the Phanerozoic had a marked influence on variations in ocean oxygenation, independent of atmospheric pO2.

The early evolutionary and much of the extinction history of marine animals is thought to be driven by changes in dissolved oxygen concentrations ([O 2 ]) in the ocean [1][2][3] . In turn, [O 2 ] is widely assumed to be dominated by the geological history of atmospheric oxygen (pO 2 ) 4,5 . Here, by contrast, we show by means of a series of Earth system model experiments how continental rearrangement during the Phanerozoic Eon drives profound variations in ocean oxygenation and induces a fundamental decoupling in time between upper-ocean and benthic [O 2 ]. We further identify the presence of state transitions in the global ocean circulation, which lead to extensive deep-ocean anoxia developing in the early Phanerozoic even under modern pO 2 . Our finding that ocean oxygenation oscillates over stable thousand-year (kyr) periods also provides a causal mechanism that might explain elevated rates of metazoan radiation and extinction during the early Palaeozoic Era 6 . The absence, in our modelling, of any simple correlation between global climate and ocean ventilation, and the occurrence of profound variations in ocean oxygenation independent of atmospheric pO 2 , presents a challenge to the interpretation of marine redox proxies, but also points to a hitherto unrecognized role for continental configuration in the evolution of the biosphere.
The availability of dissolved oxygen in the ocean exerts a critical control on habitability for marine animals 1,7 and is thought to have strongly affected their evolution from the late Neoproterozoic and through the Phanerozoic Eon (541-0 million years ago (Ma)) [1][2][3] . Indeed, the emergence of the first unambiguous metazoans during the Ediacaran Period (635-541 Ma) may have been triggered by the occurrence of ambient dissolved oxygen concentrations ([O 2 ]) sufficient to support the metabolism of a large body size 1 . Further increases in ocean oxygenation may have contributed to the Great Ordovician Biodiversification Event 8 as well as the rise of large predatory fish in the Devonian 9 . Conversely, episodic deoxygenation is considered the primary kill mechanism during many of the main Phanerozoic mass extinction events, including the Late Ordovician 10 , Late Devonian 11 and end Permian 7 .
The past few decades have seen the development and application of a variety of geochemical redox proxies, such as δ 238 U (ref. 12 ), δ 98 Mo (ref. 9 ), I/Ca (ref. 13 ) and iron speciation 14 , which have provided us with unprecedented insights into the Phanerozoic history of ocean oxygenation 5,13 . In parallel, numerical models have been instrumental in helping interpret the proxy trends and test hypotheses for the underlying driving mechanisms 4,15 . Principal among the ideas that have arisen is that changes in atmospheric pO 2 can be inferred from the ocean 4,5,15 . However, for computational reasons, this link between ocean [O 2 ] and atmospheric pO 2 has invariably been made without considering how ocean circulation may have changed over time 4,5,15 . Although the importance of continental configuration for ocean circulation and surface climate has started to be studied systematically [16][17][18] , spatial patterns of ocean geochemistry and redox have only been considered in a limited number of continental configurations (such as refs. 7,19 ) and often analysed in temporal isolation. Here we explicitly address this gap using an Earth system model of intermediate complexity to quantify how continental configuration can modulate the distribution of [O 2 ] in the ocean, thereby leading us to a radically different conclusion about the inferences that can be drawn from marine redox proxy records.
We base our analysis on a series of past ocean circulation scenarios generated using cGENIE 20 , an Earth system model of intermediate complexity designed to simulate the large-scale biogeochemical cycles and patterns in the ocean, and one that has previously been shown to be capable of simulating regional-scale distributions of [O 2 ] in the present-day ocean 20 as well as in the geological past 19,21 (see Methods and Extended Data Fig. 1). In these experiments, we consider both the potential role of changing ecological structure in the ocean (using the size-structured plankton model of ref. 22 ) as well as the influence of temperature on metabolic rates (following ref. 23 ) in creating 3D realizations of the potential distribution of [O 2 ] in the ocean. On the basis of the continental reconstructions of Scotese and Wright 24 , we conducted one simulation every 20 Myr through the Phanerozoic, for a total of 28 simulated time slices. To help isolate the impact of continental configuration from other changes occurring through the Phanerozoic, only the continental configuration 24 (plus physical atmospheric boundary Article conditions-see Methods) was varied from one time slice to the other. Solar luminosity (1,368 W m −2 ), atmospheric oxygen concentration (20.95%) and ocean nutrient inventory (2.1 μmol kg −1 PO 4 ) were kept identical in every model run, and in a first series of simulations (#1), we also kept atmospheric pCO 2 fixed (at 2,240 ppm; see Methods). Hence, we are not aiming to reconstruct the Phanerozoic evolution of climate here (see, for example, refs. 16,18 ) nor necessarily provide a reconstruction of past [O 2 ] in the ocean or even recover modern distributions with high fidelity, but, rather, expose the specific impact of changes in continental configuration on gross ocean oxygenation.
Despite assuming both invariant solar constant and pCO 2 , substantive variability in climate through the Phanerozoic occurs in the model (Fig. 1b). This is an expected result of sea-level-driven changes in the proportion of absorptive ocean surface versus more reflective land surface 17,25 , as well as of continental latitudinal distribution, and, in a principal component analysis conducted across all 28 time slices, we find a strong positive correlation between global mean sea-surface temperature (SST) and global ocean surface area ( Fig. 1e; R 2 = 0.585). Given our assumption of invariant pO 2 (deliberately chosen to expose the role of the continental configuration), we find a relatively straightforward and expected inverse relationship between mean ocean-surface [O 2 ] and SST (Fig. 1b,e)-driven primarily by changes in oxygen solubility and modulated by the latitudinal distribution of ocean surface area. Global mean subsurface (around 90-190 m depth in the model) [O 2 ] exhibits a more exaggerated temporal evolution (Fig. 1c). Respiration of organic matter exported from the ocean surface now strongly amplifies subsurface oxygen variability (Fig. 1c,e) and introduces new features unconnected to SST changes, such as the approximately 35 μmol kg −1 decrease from 460 to 400 Ma. This is our first hint of the importance of changing continental distribution and, hence, large-scale pattern of ocean circulation, which-in the case of subsurface [O 2 ]-is by means of changes in the resupply of nutrients to the surface. Even without attempting to simulate a realistic sequence of Phanerozoic climate states (and, hence, changing ocean-surface oxygen solubility), our modelling provides hints to the timing of oceanic anoxic events (OAEs) during the Phanerozoic ( Fig. 2 and Extended Data Fig. 2). For instance, seafloor ocean deoxygenation develops in the Palaeo-Tethys during the Permian-Triassic transition (ca. 260 Ma) 26 and in the central Atlantic during the Early Cretaceous (ca. 120 Ma) 27 and Late Cretaceous (ca. 100 Ma) 28 . That the Late Devonian OAE (380 Ma) 29 and Toarcian OAE (ca. 180 Ma) 30 are not captured by our simulations may be a result of the coarse model grid (for the Tethys Ocean during the Toarcian OAE in particular 30 ; Fig. 2 and Extended Data Fig. 2). This could also be interpreted as reflecting the need for a strong perturbation (warming and nutrient-driven increase in primary productivity) at that time to trigger an OAE, as the continental configuration itself is not particularly prone to ocean deoxygenation in our model simulations.
Our unexpected discovery from the numerical modelling is how variable global deep-ocean [O 2 ] (calculated as the mean of all benthic grid points in the model below 1,000 m depth) is (Fig. 1d), ranging from 30 μmol kg −1 (and near fully anoxic) to 178 μmol kg −1 (a little higher than modern). The early Palaeozoic (541 to 460 Ma) stands out in particular and is characterized by anomalously poor deep seafloor oxygenation (generally less than 45 μmol kg −1 ) in the model, with an abrupt increase to 171 μmol kg −1 at 440 Ma and mostly remaining above 100 μmol kg −1 thereafter (Fig. 1d). The step change in benthic oxygenation simulated between 460 Ma and 440 Ma is similarly observed using an alternative set of continental reconstructions 31 (Extended   The role of continental configuration in creating conditions of extreme deep-ocean deoxygenation and weakened ventilation (old water mass ages) during the early Palaeozoic (before 440 Ma in series #1 (Fig. 1d) and before 420 Ma in series #2 (Fig. 3d)), arises from the existence of state transitions in the large-scale circulation of the ocean in the model. Using series #1 as an example, the oldest time slices (540-500 Ma) are characterized by a poorly ventilated seafloor with deep waters forming over the South Pole (Figs. 1d and 2 and Extended Data Fig. 8a). At 480 and 460 Ma, the ocean circulation oscillates between the previous poorly ventilated state and a better-ventilated state, the latter state subsequently becoming stabilized at 440 Ma by further tectonic rearrangement (Figs. 1d and 2 and Extended Data Fig. 8a). Of course, the differences in the two simulation series illustrate that the climate state also influences the ocean circulation regimes (contrast Figs. 1 and 3; see also Extended Data Figs. 8 and 9). Several lines of evidence suggest that these state transitions are a robust characteristic of the early Palaeozoic, at least in the cGENIE model. Indeed, further simulations conducted using idealized land-sea masks show that these state transitions may be a characteristic feature of continental configurations with one pole free of land (Extended Data Fig. 10 and Methods), such as the one of the early Palaeozoic (Fig. 2), and a similar reorganization of the global ocean circulation has previously been shown to best explain the expansion of seafloor anoxia during Late Ordovician glaciation (Extended Data Fig. 10 and ref. 21 ).
The occurrence of climate-driven state transitions in ocean circulation and oxygenation is well established in the context of changes in Atlantic and Pacific basin meridional overturning since the last glacial maximum (21 ka) 32,33 . By contrast, the potential for ocean state transitions in deeper time has received much less attention but is no less important (see, for example, ref. 21 ). We also have good reason to  Article conclude that the occurrence of oscillatory ocean circulation modes is not an artefact of the particular model used here and may occur as a transitional regime between different stable steady states (Extended Data Fig. 10). For instance, oscillations have previously been identified in both (ocean-only) 3D ocean circulation models 34,35 and 3D ocean circulation models coupled to simplified representations of atmospheric dynamics (Earth system models of intermediate complexity) 36,37 . Non-steady solutions in volume-integrated ocean temperature are also hinted at in several recent fully coupled general circulation model (GCM) simulations under specific continental configurations and climate states 16 . Further work conducted using higher-resolution Earth system models, such as the Coupled Model Intercomparison Project phase 6 generation, will be important to distinguish robust from model-dependent modes of ocean circulation, although identifying the boundary conditions giving rise to oscillatory modes will be computationally challenging. Our finding of a fundamental decoupling of subsurface and benthic [O 2 ] during the Phanerozoic (Fig. 1e; R 2 = 0.001) implies that no simple relationship exists during the Phanerozoic among the concentrations of oxygen in the atmosphere, the subsurface ocean and at the seafloor. Furthermore, the highly non-linear nature of the system means that even very minor changes in continental configuration (for example, contrast 460 versus 440 Ma; Fig. 2) can induce pronounced reorganizations of large-scale circulation and, hence, ocean oxygenation and, similarly, confusticate any predictable relationship between changes in atmospheric pCO 2 (and climate) and ocean oxygenation. This has implications not only for how we interpret deep-ocean redox proxy records and specifically proxies for the global areal extent of seafloor anoxia 5 (such as, δ 238 U (ref. 12 ) and δ 98 Mo (ref. 9 )) but also how we understand the underlying controls on ocean anoxia and infer changes in atmospheric composition. For instance, Stolper and Keller 38 analysed the ratio of Fe 3+ to total Fe (Fe 3+ /ΣFe) in hydrothermally altered basalts formed in ocean basins to quantitatively reconstruct deep-ocean [O 2 ] over the past 3,500 Myr. Their results indicate that the deep ocean became oxygenated only in the Phanerozoic and probably not until the late Palaeozoic (<420 Myr). Our simulations support the vision of a poorly oxygenated early Palaeozoic ocean 4,5,38 (until ca. 440 Ma in the model; see Figs. 1d and 3d). However, although Fe 3+ /ΣFe data are interpreted as reflecting an increase in atmospheric oxygen concentrations 38 , our Earth system simulations, by contrast, suggest that deep-ocean oxygenation around 440-420 Ma could have been a consequence of continental rearrangement and concomitant changes in ocean circulation and ventilation. The implication is that early Phanerozoic atmospheric pO 2 need not have been appreciably lower than modern, which aligns with the measurements of atmospheric oxygen trapped in fluid-gas inclusions of halite 39 . Furthermore, the step changes in early Palaeozoic seafloor oxygenation are indicative of bistability between deep-ocean circulation states; that is, for certain continental configurations, multiple advection and convection patterns will co-exist (what is observed depends on the initial conditions of the simulation). Switches between such states could provide an explanation for the swings 5 in early Palaeozoic global anoxia extent documented on the grounds of δ 238 U (refs. [40][41][42] and δ 98 Mo (refs. 41,43 ). Uncertainties in both the model and geological data preclude us from determining if our simulated scenario may have actually taken place in the geological past, and we recognize that continental rearrangement is also not exclusive of secular variation in other drivers of ocean anoxia, such as atmospheric oxygen, ocean nutrient inventory and global warming 13,19,21,26 . That said, and despite our idealized boundary conditions and invariant atmospheric chemistry, there are shared temporal trends (albeit with different timings and interval durations) with recent I/Ca-based oxygen proxy compilations 13 , namely, low oxygenation in the early and late Palaeozoic with a high during the Devonian, and peak values occurring before late Cenozoic cooling.
Finally, our results have implications for the evolution of marine animal ecosystems through the Phanerozoic. If kyr-scale oscillations in ocean oxygenation were a characteristic feature of early Palaeozoic oceans, as supported by the high-resolution Mo record of ref. 44 and our modelling (Extended Data Fig. 8), this would lend support to arguments linking the progressive (here, palaeogeographical) stabilization of a steady-state Earth system to decreasing extinction rates through the Phanerozoic 6,45 .

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-05018-z.

Description of the model cGENIE 20 -an 'Earth system model of intermediate complexity'-is based
around a 3D ocean circulation model, which-for speed-is coupled to a 2D energy-moisture balance atmospheric component. We configured the model on a 36 × 36 equal-area grid with 16 unevenly spaced vertical levels to a maximum of depth of 5,500 m in the ocean. The cycling of carbon and associated tracers in the ocean is based on a single (phosphate) nutrient limitation of biological productivity accounting for plankton ecology based on refs. 22,46 , but adopts the Arrhenius-type temperature-dependent scheme for the remineralization of organic matter exported to the ocean interior of ref. 23 . We further modify the ecological model of ref. 22 to account for limitation of productivity under sea ice, rescale the global mean annual C/P ratio of exported particulate organic matter to match that of the standard model (106) 47 , and only diagnose the ocean-surface mixed-layer depth for calculating light attenuation rather than entrain this into the physical circulation. We include a sulphate (SO 4 2− ) tracer as an electron acceptor in the model in addition to dissolved oxygen, such that when [O 2 ] starts to approach zero, sulphate reduction occurs in association with the remineralization of organic matter. The hydrogen sulphide (H 2 S) created is transported by means of ocean circulation and oxidized in the presence of O 2 to re-form sulphate. Note that we do not include a nitrogen cycle (and, hence, nitrate reduction) nor methanogenesis in our model tracer configuration. See refs. 48,49 for details of the overall ocean redox scheme.
Despite its low spatial resolution, cGENIE satisfactorily simulates first-order ocean [O 2 ] spatial patterns and values for the present-day ocean when configured with pre-industrial boundary conditions such as 278 ppm CO 2 (Extended Data Fig. 1 and ref. 20 ), as well as in the geological past 19,21 . For instance, simulations for Oceanic Anoxic Event 2 conducted using cGENIE (ref. 19 , their Fig. 2f,g) compare well with simulations conducted using the more complex model IPSL-CM5A2 (ref. 28 , their Fig. 10), in terms of both global extent of anoxia and response to gateway alteration.

Description of the numerical experiments
We adopted the (flat-bottomed) Phanerozoic continental reconstructions of ref. 24 but substituted the deep-ocean bathymetry of ref. 50 when available (140-0 Ma) to account for mid-ocean ridges. To specifically quantify the impact of palaeogeographical evolution, we kept solar luminosity (1,368 W m −2 ), atmospheric oxygen concentrations (20.95%) and ocean nutrient inventory (2.1 μmol kg −1 PO 4 ) invariant. We used a null eccentricity-minimum obliquity orbital configuration, which provides an equal mean annual insolation to both hemispheres with minimum seasonal contrasts. Atmospheric CO 2 concentration was kept fixed in series #1 (2,240 ppm) but varied in series #2 to approximately 'correct' for the palaeogeographical impacts on climate. These forcing combinations ensured that simulated ocean temperatures are on the right order of magnitude compared with the Phanerozoic proxy-derived trend 18 . The mean Phanerozoic tropical (25° S-25° N) SST is 29.67 °C (standard deviation 1.88 °C) in series #1 (mean of 29.39 °C and standard deviation of 0.45 °C in series #2) compared with a mean of 28.1 °C (standard deviation 2.51 °C) in the data compilation 18 .
To generate the physical atmospheric boundary conditions required by cGENIE for each different continental configuration, we ran FOAM 51 GCM experiments for 2-3 kyr, until deep-ocean temperature equilibrium. We then derived the 2D wind speed and wind stress, and 1D zonally averaged albedo forcing fields required by the cGENIE model, using the 'muffingen' open-source software (https://doi.org/10.5281/ zenodo.5500687), following the methods used in ref. 52 .
Simulations were initialized with a sea-ice-free ocean and homogeneous temperature and salinity in the ocean (5 °C and 33.9‰, respectively) and integrated for a total of 20,000 years. Data are presented reflecting either the average over the last 5,000 years of the simulation (for example, Fig. 2) or the mean and range of any oscillation over the last 5,000 years (for example, Figs. 1 and 3), as detailed in the figure captions.

Principal component analysis
Principal component analyses of Fig. 1 and Extended Data Fig. 5  , global water age, benthic water age, sea-surface temperature ('SST'), benthic temperature, pCO 2 , meridional overturning circulation maximum ('MOC max') and minimum ('MOC min'), export production ('export'; that is, the flux of organic carbon produced in the ocean photic zone by primary producers that is not recycled (remineralized) before sinking deeper in the water column) and total ocean area. The principal component analysis was conducted using the 'PCA' function of the R package 'FactoMineR' 53 and the first two axes/principal components were represented (capturing about 68% of the variance in both cases). The objective here is not to provide a thorough principal component analysis but rather to (1) confirm the robustness of the correlations illustrated in the form of time series and R 2 in the other panels of Fig. 1 and Extended Data Fig. 5 and (2) provide information on the correlation (or not) of further key variables in a synthetic and graphical manner.

More discussion on ocean circulation regimes
Most of the simulations in Extended Data Fig. 8a,b show that ocean dynamics obtain stable steady-state solutions. However, among these are cases in which the system overshoots (for example, compare 0 Ma and 120 Ma in Extended Data Fig. 8a) before approaching a steady-state solution. This may imply the presence of an oscillatory mode, in agreement with the linear stability analysis of low-resolution GCM steady states under present-day conditions 34,35 . The simulations in Extended Data Fig. 8 show that, depending on changes to continental configuration, the steady-state solutions can lose stability, resulting in stable oscillations. It is known that, under certain climate conditions, advective feedbacks in the ocean circulation can destabilize the steady state 54 . Simulations for 540 Ma and 500 Ma in Extended Data Fig. 8a,b, respectively, are cases of steady-state solutions with oscillatory modes that are linearly stable, but only weakly stable, resulting in clear examples of damped oscillations. Because the frequencies of the damped and self-sustained oscillations (500 and 520 Ma in Extended Data Fig. 8b, respectively) are very similar, the transition could represent a Hopf bifurcation, resulting in self-sustained oscillations (characterized by fluctuations in circulation pattern and strength) around an unstable steady state.
To better understand the model behaviour and, ultimately, the nature of the state transitions in the global ocean circulation, we conducted further cGENIE simulations using two contrasting, idealized continental configurations featuring either a latitudinal strip of land from the North Pole to the Southern Hemisphere mid-latitudes ('Drake world' 55 ) or a strip of land extending from the North Pole to the South Pole ('ridge world' 55 ) (Extended Data Fig. 10). Drake world, with a pole free of land (that is, purely oceanic at all longitudes), schematically represents the continental configuration of the early Palaeozoic (or a high-latitude circumpolar path such as in the late Cenozoic), whereas ridge world represents an early-Mesozoic-like continental arrangement with land masses present from the North Pole to the South Pole (see Fig. 2).
With Drake world, at low pCO 2 levels, deep-water formation takes place at the North Pole, where a slightly higher hemispherical fraction of land cover with lower thermal inertia favours winter cooling over the surrounding ocean and convection (label 1 in Extended Data Fig. 10a,b). As atmospheric CO 2 is steadily increased across a series of experiments, the onset of a first set of oscillations, between ×5.5 and ×7.5 CO 2 (label 2), marks the development of a new equilibrium that stabilizes in the subsequent stable steady state at around ×8.25 CO 2 (label 3). Then, a new set of oscillations, from ×9 to ×9.75 CO 2 , marks another step towards the weakening of the deep-water formation over the North Pole and strengthening over the South Pole (label 4). A third stable equilibrium finally stabilizes above ×10.75 CO 2 (label 5). Each set of oscillations thus marks the shift to a new regime of steady-state circulation and the ultimate result is a shift in the locus of deep-water formation from the North Pole (with land barriers) to the South Pole (oceanic) in response to global warming from ×4 to ×16 CO 2 . A comparable (but latitudinally reversed) behaviour was previously observed by Pohl et al. 21 for the latest Ordovician. In response to global warming (their Fig. 4e-i; see also Extended Data Fig. 10c), the locus of deep-water formation shifts in their simulations from the South Pole (with lands) to the North Pole (oceanic). This is very similar to the change in ocean circulation simulated between 460 Ma and 440 Ma in simulation series #1 (Fig. 2), although the state transition here involves continental rearrangement in addition to global warming (global SST increases from 24.1 °C at 460 Ma to 24.38 °C at 440 Ma). No similar state transition is observed using the ridge world configuration (Extended Data Fig. 10a).

Code availability
The version of the cGENIE code used in this paper is tagged as release v0.9.31 and is available at https://doi.org/10.5281/zenodo.6823664. Necessary boundary condition files are included as part of the code release. Configuration files for the specific experiments presented in the paper can be found in the installation subdirectory: genie-userconfigs/PUBS/ published/Pohl_et_al.2022. Details of the experiments, plus the command line needed to run each one, are given in the readme.txt file in that directory. A manual describing code installation, basic model configuration and an extensive series of tutorials is provided (https:// doi.org/10.5281/zenodo.5500696). The FOAM output is hosted on Zenodo (https://doi.org/10.5281/zenodo.5780096).