Articles | Volume 15, issue 4
Clim. Past, 15, 1463–1483, 2019
Clim. Past, 15, 1463–1483, 2019

Research article 05 Aug 2019

Research article | 05 Aug 2019

Simulating the climate response to atmospheric oxygen variability in the Phanerozoic: a focus on the Holocene, Cretaceous and Permian

Simulating the climate response to atmospheric oxygen variability in the Phanerozoic: a focus on the Holocene, Cretaceous and Permian
David C. Wade1, Nathan Luke Abraham1,2, Alexander Farnsworth3, Paul J. Valdes3, Fran Bragg3, and Alexander T. Archibald1,2 David C. Wade et al.
  • 1Department of Chemistry, Centre for Atmospheric Science, Cambridge, UK
  • 2Department of Chemistry, National Centre for Atmospheric Science, Cambridge, UK
  • 3School of Geographical Sciences, University of Bristol, Bristol, UK

Correspondence: David C. Wade ( and Alexander T. Archibald (


The amount of dioxygen (O2) in the atmosphere may have varied from as little as 5 % to as much as 35 % during the Phanerozoic eon (54 Ma–present). These changes in the amount of O2 are large enough to have led to changes in atmospheric mass, which may alter the radiative budget of the atmosphere, leading to this mechanism being invoked to explain discrepancies between climate model simulations and proxy reconstructions of past climates. Here, we present the first fully 3-D numerical model simulations to investigate the climate impacts of changes in O2 under different climate states using the coupled atmosphere–ocean Hadley Centre Global Environmental Model version 3 (HadGEM3-AO) and Hadley Centre Coupled Model version 3 (HadCM3-BL) models. We show that simulations with an increase in O2 content result in increased global-mean surface air temperature under conditions of a pre-industrial Holocene climate state, in agreement with idealised 1-D and 2-D modelling studies. We demonstrate the mechanism behind the warming is complex and involves a trade-off between a number of factors. Increasing atmospheric O2 leads to a reduction in incident shortwave radiation at the Earth's surface due to Rayleigh scattering, a cooling effect. However, there is a competing warming effect due to an increase in the pressure broadening of greenhouse gas absorption lines and dynamical feedbacks, which alter the meridional heat transport of the ocean, warming polar regions and cooling tropical regions.

Case studies from past climates are investigated using HadCM3-BL and show that, in the warmest climate states in the Maastrichtian (72.1–66.0 Ma), increasing oxygen may lead to a temperature decrease, as the equilibrium climate sensitivity is lower. For the Asselian (298.9–295.0 Ma), increasing oxygen content leads to a warmer global-mean surface temperature and reduced carbon storage on land, suggesting that high oxygen content may have been a contributing factor in preventing a “Snowball Earth” during this period of the early Permian. These climate model simulations reconcile the surface temperature response to oxygen content changes across the hierarchy of model complexity and highlight the broad range of Earth system feedbacks that need to be accounted for when considering the climate response to changes in atmospheric oxygen content.

1 Introduction

The primary driver of climate over the Phanerozoic is atmospheric CO2 (Royer et al.2004). However, atmospheric oxygen content may also have varied across the Phanerozoic. Atmospheric dioxygen (O2) plays a vital role in the Earth system (Catling et al.2005), regulating the biosphere through fire ignition (Watson et al.1978) and metabolism of aerobic biota. Hence, variability in the partial pressure of dioxygen (pO2, a measure of the mass of O2 in the atmosphere, assuming N2 and the volume of the atmosphere have been constant) over time has been invoked as an evolutionary trigger (Berner et al.2007) of both animals (Falkowski et al.2005) and plants (He et al.2012) at many points in the Phanerozoic (Saltzman et al.2011; Robinson1990; Beerling and Berner2000; Scott and Glasspool2006; Edwards et al.2017).

Figure 1Oxygen content reconstructions in the Phanerozoic from Algeo and Ingall (2007), Arvidson et al. (2013), Bergman et al. (2004), Berner (2009), Berner and Canfield (1989) and Glasspool and Scott (2010). The mean (black line) and range (grey shading) of the Arvidson et al. (2013), Bergman et al. (2004) and Berner (2009) is indicated, as these reconstructions were most consistent with ice core evidence (Stolper et al.2016). Present-day atmospheric oxygen content is indicated by the solid horizontal grey line. Timings of the palaeo case studies explored in this study are indicated by the dotted vertical lines (As: Asselian, Wu: Wuchiapingian, Ma: Maastrichtian).


While strong biological and geological feedbacks prevent rapid swings in atmospheric oxygen (Catling and Claire2005), reconstructions of past atmospheric oxygen content suggest that there have been substantial excursions from the 21 % oxygen content present in today's atmosphere at times in the Phanerozoic eon. These reconstruction methods can be divided into forward and inversion models. Forward models include nutrient/weathering models (Bergman et al.2004; Arvidson et al.2013; Hansen and Wallmann2003) and isotope mass balance models (Berner2009; Falkowski et al.2005), while inversion models infer oxygen content from proxies such as charcoal (Glasspool and Scott2010), organic-carbon-to-phosphorus ratios (Algeo and Ingall2007) and plant resin δ13C (Tappert et al.2013). Figure 1 shows the reconstructed oxygen content for a variety of these methods. There is general agreement in the trends in the reconstructions, in that oxygen content increased from 5 % to 25 % in the early Paleozoic to 20 %–35 % in the Permian and subsequently stabilised at levels around 15 %–30 % from the Middle Triassic onward. However, there is uncertainty in the absolute amount of O2 for the different reconstructions (grey shading in Fig. 1). Indeed, there is support for elevated O2 by carbon isotope measurements during the Permian (Beerling et al.2002). However, disagreement is particularly evident in the Mesozoic, with low values simulated by isotope mass balance approaches. Mills et al. (2016) have shown that this could be due to an inappropriate choice of δ13C and that adjusting this value with geological constraints leads to a higher reconstructed oxygen content in better agreement with wildfire records. At the time of writing, there are no direct geochemical proxies for atmospheric pO2 on the Phanerozoic timescale. However, there is isotopic evidence of oceanic oxygenation in steps at approximately 560 (Dahl et al.2010), 400 (Dahl et al.2010; Lu et al.2018) and 200 Ma (Lu et al.2018). pO2 in the last 800 000 years has been reconstructed using O2∕N2 ratios in ice cores (Stolper et al.2016). A roughly 7 ‰ decline in pO2 is consistent with the ability to change oxygen content by the order of a few percent in ∼10 Myr. The reconstructions of Bergman et al. (2004), Arvidson et al. (2013) and Berner (2009) are the most plausible based on ice core data (Stolper et al.2016). Considering these three models alone would still suggest a large uncertainty in oxygen content for most of the Phanerozoic, except for elevated levels in the late Carboniferous/early Permian and reduced levels in the late Devonian.

Phanerozoic means “visible life” and one of the marked changes to carbon cycling between the Proterozoic and Phanerozoic was caused by the emergence of land plants. The radiation of land plants has led to strong regulation of atmospheric CO2 and O2 which both play important roles in photosynthesis. Land plants likely led to a substantial sequestration of carbon in the terrestrial biosphere and possibly led to the Ordovician glaciation (Lenton et al.2012). Increases in organic carbon sequestration in the aftermath of the evolution of lignin production may also have contributed to the cooling (Robinson1989). This fundamental change to the Earth system may have constrained CO2 levels to between 10 and 200 Pa ever since (Franks et al.2014). Watson et al. (1978) have argued that strong fire feedbacks prevent large fluctuations in oxygen levels, due to runaway burning at high oxygen levels. However, subsequent experiments using natural fuels support the possibility of the Earth system to support higher oxygen levels (Wildman et al.2004). Charcoal appears in the fossil record continuously since the late Devonian (∼360 Ma; Algeo and Ingall2007; Scott and Glasspool2006). This suggests a floor on oxygen levels in the region of 12 % (Wildman et al.2004) to 16 % (Belcher and McElwain2008; Belcher et al.2010) since then due to limits on ignition.

Variations in pO2 also have important implications for photosynthesis and therefore the operation of the terrestrial carbon cycle. The primary CO2-fixing enzyme, RuBisCO, possesses a dual carboxylase–oxygenase function (Smith1976). A photosynthetic carboxylase pathway removes CO2 from the atmosphere, while oxygenation leads to photorespiration and CO2 evolution. Therefore, increases in pO2 ought to lead to O2 outcompeting CO2 for active sites on the RuBisCO enzyme and leading to a reduction in net primary productivity (less photosynthesis, more respiration). However, photorespiration is likely to be necessary for removal of harmful byproducts in the photosynthetic metabolic pathway (Hagemann et al.2016), and a recent study suggests that increases in photorespiration may actually promote photosynthesis (Timm et al.2015). Photosynthesis is itself sensitive to the background CO2 content (Beerling and Berner2000). In addition, temperature modifies the relative solubilities of CO2 and O2 (Jordan and Ogren1984). Temperature also affects the specificity of RuBisCO for CO2 (Long1991). Therefore, the coevolution of pO2, pCO2 and temperature across the Phanerozoic has the capacity to significantly impact the terrestrial carbon cycle.

This paper focuses on investigating the climate impacts of atmospheric mass variation as the result of altering the concentration of O2. Lower atmospheric mass leads to less Rayleigh scattering so more shortwave radiation reaches the Earth’s surface. This enhances atmospheric convection and the hydrological cycle, which leads to more tropospheric water vapour, further enhancing warming. However, lower atmospheric mass leads to a reduction in the pressure broadening of greenhouse gas absorption lines which should lead to a weaker greenhouse effect and lead to cooling. Previous modelling studies have investigated which factor dominates with conflicting results. Goldblatt et al. (2009) presented radiative–convective model simulations for the Archean (∼ 3 Ga), which suggested that a nitrogen inventory around 3 times larger than present would help to keep the early Earth warm at a time when solar input was only around 75 % of what it is today, potentially solving the “Faint Young Sun” paradox (Feulner2012). Charnay et al. (2013) investigated this using a general circulation model (GCM) coupled to a slab ocean and found that for their idealised early Earth simulations they achieved a strong warming (+7 C) in response to a doubling in atmospheric mass. Poulsen et al. (2015) simulated the climate impacts of changes in O2 content over a range of 5 %–35 % using the GENESIS climate model with a slab ocean and a continental configuration consistent with the Cenomanian (mid-Cretaceous, 95 Ma) and found the opposite response – lower atmospheric mass at low pO2 was associated with a strong warming. Subsequent 1-D calculations cast doubt on this result (Goldblatt2016; Payne et al.2016); however, it is plausible that other climate feedbacks such as changes to relative humidity and cloud changes may be important as atmospheric mass changes. These would not be accounted for in 1-D radiative–convective equilibrium simulations. Cloud feedbacks in particular are a good candidate for explaining the discrepancy as cloud feedbacks under CO2-driven climate change have strong model dependency (Bony et al.2015). Another feedback which has not been considered is the possible impact of changes in the mechanical forcing of wind on the ocean circulation. In the absence of this effect, Earth’s surface temperature would be 8.7 K cooler (Saenko2009) and the Equator-to-pole temperature gradient would be steeper. Wind stress (τ) is parameterised in GCMs as τ=ρuu, where ρ is the atmospheric density and u is the surface wind vector. So, as atmospheric density increases, the wind stress on the ocean and therefore ocean heat transport should increase accordingly. Increased meridional heat transport in high-density atmospheres is also supported by an idealised 2-D modelling study of the early Earth (Chemke et al.2016). As slab ocean models assume a constant or diffusive ocean heat transport, the Charnay et al. (2013) and Poulsen et al. (2015) simulations cannot account for these effects.

As pO2 variability may alter the radiative budget of the atmosphere, it may also have impacts on the sensitivity of the climate state to CO2 changes. The equilibrium climate sensitivity (ECS) is a metric for the sensitivity of a climate model to an abrupt doubling of atmospheric CO2. Understanding this value is important for predictions of both past and future climatic changes. As the radiative forcing of CO2 is approximately logarithmic with concentration, theoretically the ECS should be constant in time as carbon dioxide changes. However, there is growing evidence that ECS has not been constant over Earth's history (Caballero and Huber2013). Changes to the incoming solar radiation (Lunt et al.2016), palaeogeography (Lunt et al.2016), CO2 levels themselves (Meraner et al.2013) and tropical sea-surface temperatures (Caballero and Huber2013) may lead to changes in the sensitivity of a particular climate state to changes in CO2.

2 Methods and simulations

2.1 Models

The impact of oxygen content variability is investigated with two coupled atmosphere–ocean general circulation models (AOGCMs): Hadley Centre Coupled Model version 3 (HadCM3-BL) and Hadley Centre Global Environmental Model version 3 (HadGEM3-AO).

HadGEM3-AO is an AOGCM (Nowack et al.2014). The atmosphere component is the UK Met Office Unified Model version 7.3 (Davies et al.2005) in the HadGEM3-A r2.0 climate configuration (Hewitt et al.2011). It employs a regular Cartesian grid of 3.75 longitude by 2.5 latitude (N48). In the vertical, 60 hybrid height vertical levels are employed – “hybrid” indicating that the model levels are sigma levels near the surface, changing smoothly to pressure levels near the top of the atmosphere (Simmons and Strüfing1983). The model top is 84 km, which permits a detailed treatment of stratospheric dynamics. A 20 min time step is used. The model employs a non-hydrostatic and fully compressible dynamical core, using a semi-implicit semi-Lagrangian advection scheme on a staggered Arakawa C grid (Arakawa and Lamb1977). Radiation is represented using the Edwards and Slingo (1996) scheme with six shortwave and nine longwave bands, accounting for the radiative effects of water vapour, carbon dioxide, methane, nitrous oxide and ozone. The Met Office Surface Exchange Scheme 2 (MOSES2) land surface scheme is used (Cox et al.1999), which simulates atmosphere–land exchanges and hydrology. A fixed present-day vegetation distribution of plant functional types is employed. The ocean component of the model is NEMO-OPA (Nucleus for European Modelling of the Ocean - Océan Parallélisé; Madec2008) model version 3.0 (Hewitt et al.2011), run at a 96 min time step. In the vertical, 31 model levels are used, which increase in thickness steadily between 10 m in the shallowest to 500 m in the deepest layer at 5 km in depth. NEMO employs a tripolar, locally anisotropic grid (ORCA2; Madec2008) which permits a more detailed treatment of the north polar region and higher resolution in the tropics. This yields an approximate horizontal resolution of 2 in both longitude and latitude, with an increased resolution of up to 0.5 in the tropics. The sea-ice component of the model is CICE (Los Alamos Community Ice CodE) at version 4.0 (Hunke and Lipscomb2008), run at a 96 min time step. This treats sea ice in a five-layer model, allowing the simulation of different ice types. The atmosphere and ocean/sea-ice components exchange fields every 24 h, while NEMO and CICE exchange fields every time step. HadGEM3-AO can be thought of as a close relation to the newest-generation HadGEM3 coupled model that will be used to support the next Intergovernmental Panel on Climate Change (IPCC) assessment (Williams et al.2018) and so represents the state of science in numerical climate models.

HadCM3-BL (Valdes et al.2017) is an AOGCM coupled to an interactive vegetation model. The model was originally developed by the United Kingdom Met Office Hadley Centre (Pope et al.2000) but has since been substantially developed further by the University of Bristol. The atmosphere component of the model employs a regular Cartesian grid of 3.75 longitude by 2.5 latitude. In the vertical, 19 hybrid height vertical levels are employed. A 30 min time step is used. The primitive equation set of White and Bromley (1995) is solved to conserve energy and angular momentum, solved on a staggered Arakawa B grid (Arakawa and Lamb1977) in the horizontal. Radiation is represented using the Edwards and Slingo (1996) scheme with six shortwave and eight longwave bands, accounting for the radiative effects of water vapour, carbon dioxide and ozone, amongst other radiative active species. The ocean component of the model employs the same horizontal grid as the atmosphere component of the model, 3.75 longitude by 2.5 latitude. In the vertical, 20 model levels are used which increase in depth from 10 m in the shallowest layer to 616 m in the deepest layer. A time step of 60 min is employed and the ocean and atmosphere components exchange required fields once per day. The ocean component is based on the Cox (1984) model, solving the full primitive equation set in three dimensions. A staggered Arakawa B grid is employed in both atmosphere and ocean models. Sea ice is treated as a zero thickness layer on the surface of the ocean grid. Ice is assumed to form at the base at a freezing point of 1.8 C but can also form from freezing in ice leads and by falling snow. A simple parameterisation of sea-ice dynamics is also employed (Gordon et al.2000) and sea-ice formation due to convergence from drift being limited to 4 m thick. Sea-ice albedo is fixed at 0.8 for temperatures below 10 C, decreasing linearly to 0.5 at 0 C. The MOSES2.1 land surface model is employed to simulate the fluxes of energy and water between the land surface and the atmosphere (Cox et al.2000; Essery et al.2003). TRIFFID (Top-down Representation of Interactive Foliage and Flora Including Dynamics; Cox et al.1998) predicts the distribution of vegetation using a plant functional type (PFT) approach. TRIFFID is run in equilibrium mode with averaged fluxes calculated over a 5-year period. TRIFFID calculates vegetation properties for five PFTs: broadleaf trees, needleleaf tree, C3 grass, C4 grass and shrubs. Grid boxes can contain a mixture of PFTs based on a “fractional coverage co-existence approach” (Valdes et al.2017). Net primary productivity (NPP) is also calculated, using a photosynthesis-stomatal conductance model (Cox et al.1998) accounting for a number of factors including atmospheric oxygen content, which affects the photorespiration compensation point (Clark et al.2011). The predicted vegetation distribution impacts the atmosphere component by altering surface albedo, evapotranspiration and surface roughness.

Figure 2Annual average surface air temperature simulated in (a) PI-CM21, (b) Ma-CM21, (c) Wu-CM21 and (d) As-CM21. Global-mean values are indicated in the top right. The 0 C isoline is indicated in pink. Continental outlines are indicated with a solid black line.

Figure 3Annual average precipitation simulated in (a) PI-CM21, (b) Ma-CM21, (c) Wu-CM21 and (d) As-CM21. Global-mean values are indicated in the top right. Continental outlines are indicated with a solid black line.

2.2 Boundary conditions

Both models simulated the climate response to oxygen variability in a pre-industrial Holocene (PIH) climate. HadCM3-BL was additionally run for three time periods across the Phanerozoic: the Maastrichtian (late Cretaceous, 66.0–72.1 Ma), Wuchiapingian (late Permian, 254.14–259.1 Ma) and the Asselian (early Permian, 295.0–298.9 Ma). The continental reconstructions employed were developed by and are from ©Getech. These reconstructions have been widely employed in a number of previous studies using the HadCM3-BL climate model (e.g. Lunt et al.2016). All three are periods of time in which models have suggested that atmospheric oxygen may have deviated significantly from the present level of 21 % (see Fig. 1). Modifications were made to alter the oxygen content of the atmosphere by adjusting the mass mixing ratios of major and minor gases, the surface pressure and other physical characteristics of the atmosphere such as the specific gas constant in an analogous way to Poulsen et al. (2015). Figure 2 shows the annual average surface temperatures and Fig. 3 shows the annual average precipitation for the (a) PI-CM21, (b) Ma-CM21, (c) Wu-CM21 and (d) As-CM21 simulations.

Table 1Summary of experiments. Experiment names AA-BBB include the continental configuration (AA) and model used (BBB). Experiment names N × AA-BBB indicate a multiplier of CO2 with respect to AA-BBB. A star (*) indicates that the CO2 multiplier was applied instantaneously and the transient adjustment to climate was analysed for the purpose of a Gregory et al. (2004) analysis.

Download Print Version | Download XLSX

A summary of the experiments performed can be found in Table 1. When an experiment with a particular oxygen content is referred to, it will be indicated in superscript; e.g. EXP21 indicates a 21 % oxygen simulation. The 21 % simulations (PI-CM21, As-CM21, Ma-CM21 and Wu-CM21) were integrated for 50 model years, as these simulations had already been spun up at that CO2 content. For 10 % and 35 % pO2, the model was spun off the 21 % pO2 simulation and iterated for at least 1000 model years. The 2 × PI-CM*, 2 × Ma-CM* and 2 × As-CM* experiments were spun off from the end of the PI-CM, Ma-CM and As-CM experiments and iterated for 100 years in order to perform a Gregory et al. (2004) analysis. For HadGEM3-AO, model integrations (PI-GEM35, PI-GEM10, 4 × PI-GEM35 and 4 × PI-GEM10) were performed for 300 model years with a 10× acceleration of the deep ocean to reduce the time for equilibrium, then integrated for a further 500 years to spin up the shallow ocean without acceleration. The last 50 years were used for model analysis.

Pre-Quaternary pCO2 is poorly constrained due to the absence of glacial ice; however, there is growing evidence that CO2 is unlikely to have been significantly higher than the order of hundreds of Pa since the radiation of land plants (Breecker et al.2010; Franks et al.2014). For the Maastrichtian, 56 Pa is used in agreement with stomatal proxy-based reconstructions (Steinthorsdottir and Pole2016). For the Asselian, 28 Pa is used in agreement with carbonate and fossil plant reconstructions (Montañez et al.2007). For the Wuchiapingian, 112 Pa is used (Brand et al.2012).

1-D atmospheric chemistry simulations have simulated higher O3 column with increasing pO2 (Kasting et al.1979; Payne et al.2016). More detailed 2-D model simulations, which capture critical latitudinal gradients in photolysis and zonal-mean transport (Hadjinicolaou and Pyle2004; Haigh and Pyle1982), support a monotonically increasing ozone column with increasing pO2 (Harfoot et al.2007). However, simulated ozone column was more sensitive to N2O levels than pO2 (Harfoot et al.2007). In addition, while column ozone reduces at low pO2 in Harfoot et al. (2007), there are increases in ozone concentration in the tropical tropopause region where the radiative effect of O3 is stronger (Forster and Shine1997). Changes in lightning are important for understanding future changes in tropospheric ozone (Banerjee et al.2014); however, they are subject to considerable uncertainty (Finney et al.2018). There may be more lightning at high pO2 due to a higher pO2pN2 ratio or less due to reduced convection (Goldblatt et al.2009). Low pO2 may also enhance isoprene emissions (Rasulov et al.2009), which could enhance tropospheric ozone and alter cloud properties (Kiehl and Shields2013). Ozone is also sensitive to changes in CH4 and N2O; the changes to inventories of these chemically active species on the Phanerozoic timescale (Beerling et al.2009, 2011) are highly uncertain. Ozone is also sensitive to dynamical changes. Given these large uncertainties in possible changes to chemical sources, reactivity and transport, we neglected including changes in atmospheric ozone concentration in these simulations. However, we recommend that follow-up work should focus on this specific question in detail. In HadGEM3-AO, the mass of tropospheric and stratospheric ozone is fixed at PIH values simulated by Nowack et al. (2014) using a tropopause height-matching scheme. This prevents a rising tropopause leading to stratospheric levels of ozone existing in the troposphere, particularly in the 4 × PI-GEM experiments. Not accounting for a rising tropopause has been found to artificially increase climate sensitivity (Heinemann2009) and initial tests not accounting for this led to instability for 4 × PI-GEM10. In HadCM3-BL, tropospheric ozone is set to 6 ppbv and stratospheric ozone is set to 1.66 ppmv for the 21 % simulations. These values are adjusted to conserve total ozone mass in the alternative pO2 scenarios.

2.3 Data

Data for the Cenomanian (Poulsen et al.2015) 21 % O2 and 10 % O2 simulations were obtained from (last access: 21 February 2018). At the time of writing, the 35 % simulation contained missing data, so it was not used for analysis.

2.4 1-D energy balance model

A 1-D energy balance model (EBM) has been used to deconvolve the contributions from changes in different parts of the climate system. This 1-D EBM approach has been applied to zonal-mean quantities for climate simulations of the Eocene by Heinemann et al. (2009) following Budyko (1969) and Sellers (1969):

(1) SW t ( ϕ ) 1 - α ( ϕ ) - 1 2 π R 2 cos ( ϕ ) F ( ϕ ) ϕ = ϵ ( ϕ ) σ T s , ebm 4 ( ϕ ) ,

where SWt is the incident shortwave radiation at the top of the atmosphere, ϕ is the latitude, α is the surface albedo, R is the radius of Earth, ϵ is the effective surface emissivity, and Ts, ebm is the EBM surface temperature (Heinemann et al.2009). F(ϕ)ϕ is the divergence of total meridional heat transport and is given by

(2) F ( ϕ ) ϕ = - 2 π R 2 cos ( ϕ ) SW t net ( ϕ ) + LW t net ( ϕ ) ,

where SWtnet and LWtnet are the net top-of-atmosphere shortwave and longwave radiative fluxes, respectively (positive downward; Heinemann et al.2009). Solving for Ts, ebm, the EBM surface temperature for each latitude can be calculated using zonal- and annual-mean radiative fluxes from the GCM. Where clear-sky radiative fluxes are also available, cloud radiative effects can be deconvolved from clear-sky radiative effects. The clear-sky albedo αc and clear-sky effective surface emissivity ϵc can be calculated by

(3) α c = SW t , c SW t , ϵ c = LW t , c LW s ,

where SWt,c is the upward top-of-atmosphere clear-sky shortwave radiative flux and LWt,c is the upward top-of-atmosphere clear-sky longwave radiative flux. When considering the temperature change between two experiments, the contributions from different components can be quantified by calculating Ts, ebm with different combinations of components from each experiment (Heinemann et al.2009).

2.5 Climate sensitivity

To estimate the climate sensitivity to CO2 changes, the linear regression methodology of Gregory et al. (2004) is employed. This assumes a linear relationship between the changes in global, annual-mean radiative imbalance at the top of the atmosphere (N, W m−2) and surface temperature anomalies (Ts, C):

(4) N = F + ξ Δ T s ,

where ξ is the effective climate feedback parameter (Wm-2C-1) and F is the effective forcing (W m−2) accounting for fast climate adjustments and effective radiative forcing. The effective ECS is then ΔTs when N=0. While there are weaknesses of this approach, particularly due to non-linearities in ξ as ΔTs changes (Armour et al.2013; Li et al.2013), the climate response when simulations are continued to equilibrium shows an accuracy to within 10 % (Li and Sharma2013). Furthermore, the contributions to ξ and F from longwave (LW) and shortwave (SW), clear-sky (CS) and cloudy-sky (CRE) components can be decomposed by a linear decomposition as

(5) F = F CS , SW + F CS , LW + F CRE , SW + F CRE , LW

for the effective forcing and

(6) ξ = ξ CS , SW + ξ CS , LW + ξ CRE , SW + ξ CRE , LW

for the effective climate feedback parameter.

3 Results

Where results are presented from a single simulation, the oxygen content for that run is written in superscript, i.e. EXP21 indicates a 21 % oxygen simulation. Where results are presented as an anomaly between simulations with different oxygen content, EXP021 indicates that the quantity presented is EXP21 minus EXP0. A summary of results is shown in Table 2.

Table 2Summary of results for EXP21, then EXP1035. Where applicable, results calculated for the Poulsen et al. (2015) Cenomanian 21 %–10 % oxygen simulation are also presented. Abbreviations: Teq-pole (Equator-to-pole surface air temperature gradient), Teq-pole,cold month (Equator-to-pole surface air temperature gradient for cold months), EBM (quantities obtained using a Budyko–Sellers 1-D energy balance model following Heinemann, 2009), Ts,ebm (EBM surface temperature), Ts,csky,ebm (EBM surface temperature change accounting for changes in clear-sky radiative fluxes), Ts,cre,ebm (EBM surface temperature change accounting for changes in cloudy-sky radiative fluxes), Ts,mht,ebm (EBM surface temperature change accounting for changes in meridional heat flux divergence).

Download Print Version | Download XLSX

3.1 Surface climate

Figure 4 (left) shows the annual-mean surface air temperature differences between the 35 % and 10 % runs. For the pre-industrial Holocene, PI-GEM1035 (Fig. 4a) shows a global-mean surface temperature response of +1.50 C, while PI-CM1035 (Fig. 4b) shows a similar global-mean surface temperature response of +1.22 C. It is worth noting that HadCM3-BL and HadGEM3-AO are not completely distinct climate models, for instance, sharing the Edwards and Slingo (1996) radiation scheme, so this is unlikely to capture the full variability in possible climate model responses. That the results are in reasonable agreement with the 1-D results of Payne et al. (2016), who simulated a temperature response between +1.05 and +2.21 C depending on assumptions about atmospheric ozone, gives some confidence in the HadCM3-BL and HadGEM3-AO results. Similarly, the As-CM1035 (Fig. 4c) case exhibits a global-mean surface temperature response of +1.29 C. For the warmest climates, a response of 0.82C is simulated for Wu-CM1035 (Fig. 4e) and +0.17 C for 4 × PI-GEM1035 (Fig. 4f). In the Ma-CM1035 case (Fig. 4d), a global-mean surface temperature response of +0.70 C is simulated. This suggests that the climate response to pO2 variability depends on the background climate state.

Figure 4Surface air temperature change for (a) PI-GEM1035, (b) PI-CM1035, (c) As-CM1035, (d) Ma-CM1035, (e) Wu-CM1035 and (f) 4 × PI-GEM1035 in the annual mean (left), cold-month mean (change in the mean grid box temperature of the coldest month in the monthly mean climatology, middle) and warm-month mean (change in the mean grid box temperature of the warmest month in the monthly mean climatology, right). The change in global-mean values (C) are offset to the top right of each plot. Note the strong high-latitude warming in the cold-month mean and tropical cooling in the warm-month mean.

There is a strong seasonal dependence in the surface air temperature response. Considering the changes in coolest average monthly temperature in each grid box (Fig. 4, middle column), the change in cool months dominates the warming response, particularly at high latitudes. By contrast, the warm-month mean is smaller/less negative in all cases (Fig. 4, right). A cooling of continental land masses is evident in the tropics and particularly in the Wu-CM (Fig. 4e) and 4 × PI-GEM (Fig. 4f) cases. These could be in part due to free-air lapse rate changes which should be stronger at high pO2; as for a given topographic height, the change in pressure is higher for high pO2, which should lead to a larger temperature reduction with height. The changes to the seasonal cycle are consistent with the radiative changes associated with changing oxygen content. The reduction in incident surface shortwave radiation should have its strongest effect on extratropical temperatures in the summer; therefore, the Rayleigh scattering component will most strongly affect the warm-month temperature. Warming from pressure broadening of greenhouse gas absorption lines as atmospheric mass increases will be most evident in extratropical winter, as with anthropogenic climate change, due to sea-ice and surface heat flux changes (Dwyer et al.2012). The reduction in the amplitude of the seasonal cycle in temperature simulated by both HadGEM3-AO and HadCM3-BL is therefore supported by a consideration of the changes to atmospheric radiation.

Figure 5Zonally and annually averaged surface air temperature difference (solid lines) from 10 % to 35 % oxygen content for PI-GEM (blue), 4 × PI-GEM (red), As-CM (pink), Ma-CM (green) and Wu-CM (purple). The difference from the annual mean to cold-month mean for each run is indicated by the shading. Values are smoothed by a Savitzky–Golay filter (Savitzky and Golay1964).


The zonal- and annual-mean surface air temperature changes are shown in Fig. 5. The Northern Hemisphere Equator-to-pole temperature gradient is reduced by 6.6 C in the PI-GEM1035 case (blue line) and 4.0 C in the PI-CM1035 case (not shown). The zonal structure of the surface temperature change is similar in the palaeoclimate case studies. In the Maastrichtian, the Equator-to-pole temperature gradient is reduced by 2.0 C (Ma-CM1035) and in the Asselian the Equator-to-pole temperature gradient is reduced by 2.3 C (As-CM1035). The Equator-to-pole temperature gradient reduces even in the Wu1035 case despite the reduction in global-mean surface temperatures.

Figure 6Annually averaged total precipitation change from 10 % to 35 % oxygen content for (a) PI-GEM1035, (b) PI-CM1035, (c) As-CM1035, (d) Ma-CM1035, (e) Wu-CM1035 and (f) 4 × PI-GEM1035. Global-mean values (mm d−1) are offset to the top right of each plot.

The hydrological cycle is also affected by changing oxygen content. Increases in Rayleigh scattering at high pO2 ought to reduce incident shortwave at the Earth's surface (Poulsen et al.2015) and inhibit convection (Goldblatt et al.2009), which should lead to reductions in precipitation. This is analogous to stratospheric sulfate or solar radiation management geoengineering where precipitation is reduced in geoengineering experiments with respect to an unperturbed climate with the same global-mean surface temperature (Irvine et al.2016). Poulsen et al. (2015) simulated large reductions in precipitation as pO2 increased in the GENESIS climate model; however, much of this could be explained by the surface temperature changes. Annually averaged precipitation change between the 10 % and 35 % oxygen content runs is shown in Fig. 6. In all cases, increasing oxygen content leads to a decline in global-mean total precipitation, despite the increase in surface temperatures, however, with strong regional differences. For PI-GEM (Fig. 6a), PI-CM (Fig. 6b) and 4 × PI-GEM (Fig. 6f), there is a clear northward shift in the tropical rainbelts. A northward shift in the Intertropical Convergence Zone (ITCZ) would be consistent with stronger warming in the Northern Hemisphere due to Bjerknes compensation (Bjerknes1964; Broccoli et al.2006). Heat transport is more hemispherically symmetric in the Maastrichtian, Asselian and Wuchiapingian cases so latitudinal ITCZ shifts are not evident. While global precipitation is reduced in Wu-CM35, the increase in ocean–land temperature contrast leads to a significant increase in tropical land precipitation. Despite the increases in global-mean surface temperatures simulated for most cases, precipitation is still reduced in all simulations.

Comparing the surface temperature and precipitation response between HadCM3-BL and HadGEM3-AO suggests that the model responses are broadly consistent. A grid-box-by-grid-box comparison of annual-mean surface air temperature and precipitation anomalies for PI-GEM1035 vs. PI-CM1035 is presented in Fig. S1 in the Supplement. The largest discrepancy in surface air temperature response between the two models occurs for the largest temperature changes simulated by HadGEM, which are strongest in Northern Hemisphere polar regions. This could be linked to differences in the representation of polar climate processes and amplification by polar ice feedbacks between the two models. There is broad consistency in cold- and warm-month means (Fig. 4a and b) with stronger warming in the cold-month mean and terrestrial cooling in the warm-month mean.

3.2 Energy balance decomposition

The drivers of the changes in surface temperature can be understood by decomposing the terms which contribute to surface temperature change in a 1-D energy balance model following Heinemann et al. (2009). For PI-CM1035, these results are shown in Fig. 7. These show that the 1-D EBM can reasonably capture the temperature response in the HadCM3-BL simulations, with slight errors (where the black and grey lines are not overlapping) evident in the polar regions. This could be due to averaging over the polar rows in the HadCM3-BL model. There are positive contributions to the surface temperature change in the clear-sky emissivity and albedo at the poles. This is consistent with the increase in pressure broadening of absorption lines and the simulated reduction in sea-ice extent. By contrast, extrapolar contributions to clear-sky albedo provide a negative contribution to the temperature change which is consistent with an increase in Rayleigh scattering which would be expected to be strongest in the tropics where the maximum in incoming solar radiation is located. Combined, the clear-sky component of the temperature change is +1.45C and the cloudy-sky component is 0.35C. This suggests that HadCM3-BL supports a cloud feedback which acts to cool the climate at high pO2 and partially offset the clear-sky temperature changes.

Figure 7The 1-D EBM decomposition for PI-CM1035. (a) EBM results (grey) vs. GCM results (black). (b) Decomposition of EBM into the emissivity (purple), albedo (green) and heat transport (orange) components of the temperature change. (c) Clear-sky emissivity (dark purple) and clear-sky albedo (dark green) components of the EBM. The all-sky components are included for comparison. (d) Decomposition of EBM into the total clear-sky (blue), cloudy-sky (red) and all-sky (grey) components.


Figure 81-D EBM decomposition for PI-GEM1035. (a) EBM results (grey) vs. GCM results (black). (b) Decomposition of EBM into the emissivity (purple), albedo (green) and heat transport (orange) components of the temperature change. (c) Clear-sky emissivity (dark purple) and clear-sky albedo (dark green) components of the EBM. The all-sky components are included for comparison. (d) Decomposition of EBM into the total clear-sky (blue), cloudy-sky (red) and all-sky (grey) components.


The same analysis was performed for the HadGEM3-AO PIH simulations. Figure 8 shows that a somewhat weaker cloud feedback is simulated by HadGEM3-AO (0.21 C). Clear-sky contributions are slightly stronger (+1.51 C). The largest differences between the simulations appear in the all-sky albedo and emissivity changes, where there appear to be competing factors which lead to a similar climate response possibly related to partitioning between the longwave and shortwave contributions to the cloud response. This is perhaps unsurprising, as cloud feedbacks to CO2 changes represent a large uncertainty in future climate change projections, and given the relatively small global-mean temperature changes, a relatively small change in cloud radiative effects has the power to considerably mediate the climate response. However, the qualitative agreement in latitudinal structure of the clear-sky albedo and emissivity changes between these structurally different models gives some confidence that the relevant climate feedbacks are well captured in these simulations.

Analysis of the palaeo case studies (As-CM, Fig. S2; Ma-CM, Fig. S3; Wu-CM, Fig. S4) shows a similar pattern. In all simulations, irrespective of surface temperature response, the clear-sky emissivity is a positive contribution to global-mean surface temperature change, while clear-sky albedo is a more negative contribution. The emissivity contribution becomes less positive as pCO2 increases from As-CM to Ma-CM to Wu-CM. By contrast, the albedo contribution becomes more negative as pCO2 increases. This is consistent with the reduction in planetary albedo as sea-ice extent is reduced on the ocean and dark vegetated surfaces increase on the land.

Figure 9Gregory analysis of HadCM3-BL: regression of top-of-atmosphere radiative imbalance against surface air temperature change (solid lines) for the first 100 years of 2 × PI-CM10 (pink) and 2 × PI-CM35 (blue) cases. Annual averages are indicated by crosses and decadal averages are indicated by filled circles. The regression was performed on the decadal averages.


Figure 10Gregory analysis of HadCM3-BL: regression of top-of-atmosphere radiative imbalance against surface air temperature change (solid lines) for the first 100 years of 2 × As-CM10 (pink) and 2 × As-CM35 (blue) cases. Annual averages are indicated by crosses and decadal averages are indicated by filled circles. The regression was performed on the decadal averages.


3.3 Climate sensitivity

Here, we investigate the impact of oxygen variability on climate sensitivity. The HadGEM3-AO and HadCM3-BL results suggest that increasing CO2 content leads to a reduction in the surface temperature change on increasing pO2 (compare 4 × PI-GEM and PI-GEM in Fig. 4). For reference, HadGEM3 has a climate sensitivity of +3.6 C (Nowack et al.2014), while HadCM3 has a climate sensitivity of +3.1 C (Johns et al.2006). From the 4 × PI-GEM and PI-GEM experiments, a reduction in climate sensitivity of 0.65 C can be inferred based on the changes in surface temperatures. For HadCM3, CO2-doubling experiments were performed, and a regression of the change in top-of-atmosphere radiative imbalance against change in surface temperature following Gregory et al. (2004) (see also Sect. 2.4) is shown in Fig. 9. The PI-CM35 climate state has a smaller ECS than PI-CM10 by 0.7 C. While the changes in total radiative forcing F are very similar, ξ is less negative (1.08 vs. 1.37 Wm-2C-1) at low pO2. The decomposition of these changes into their longwave and shortwave components, clear-sky and cloudy-sky components is also shown in Fig. 9. The clear-sky longwave radiative flux changes are slightly higher in 2 × PI-CM35 (4.0 W m−2) than 2 × PI-CM10 (3.8 W m−2) as would be expected due to the pressure broadening of CO2. The clear driver for the less negative ξ value are from the longwave cloud radiative effect changes, which is much steeper for 2 × PI-CM10 (+0.62 Wm-2C-1) than 2 × PI-CM35 (+0.17 Wm-2C-1). This is somewhat offset by stronger clear-sky shortwave radiative feedbacks in 2 × PI-CM35 (+1.00 Wm-2C-1) than 2 × PI-CM10 (+0.57 Wm-2C-1). This highlights the important role that cloud radiative feedbacks play in determining the climate sensitivity.

Figure 11Gregory analysis of HadCM3-BL: regression of top-of-atmosphere radiative imbalance against surface air temperature change (solid lines) for the first 100 years of 2 × Ma-CM10 (pink) and 2 × Ma-CM35 (blue) cases. Annual averages are indicated by crosses and decadal averages are indicated by filled circles. The regression was performed on the decadal averages.


An increase in ECS appears to be robust across the HadCM3-BL experiments. For As-CM, ECS is 0.8 C lower at 35 % O2 than 10 % O2 (Fig. 10). For Ma-CM, this value is much larger. A 3.3 C reduction in ECS is simulated, which is also driven by the longwave cloud radiative effects in conjunction with a weaker clear-sky shortwave radiative effect which tended to cool the low pO2 (Fig. 11). Unlike the clear-sky shortwave effects, the longwave cloud radiative effects seem consistent across the three experiments. It should be noted that attempts were made to simulate 2× experiments for the Wuchiapingian; however, what would have been the 2 × Wu-CM10, in the nomenclature used here, was numerically unstable.

Figure 12Change in column water vapour in (a) PI-CM1035, (b) As-CM1035, (c) Ma-CM1035 and (d) Wu-CM1035. Global sum values (petagrams) are offset to the top right of each plot. Note the atmospheric drying at high pO2 is enhanced in the warmer climate states of the Wuchiapingian and Maastrichtian and more subdued in the cooler climate states of the Asselian and Holocene.

The increase in climate sensitivity appears to be linked to the reduction in temperature anomaly in a warmer climate state. We propose that this is due to more vigorous convection at low pO2 (Goldblatt et al.2009), leading to an atmospheric moistening (Fig. 12) which causes warming analogously to Rose and Ferreira (2013). This is consistent with the increases in climate sensitivity observed – in a warmer climate the atmosphere can hold more water vapour, so any changes to water vapour will be amplified in their impacts on the radiative budget of the atmosphere. This water vapour feedback is also consistent with the weaker clear-sky shortwave radiative effect observed in 2 × Ma-CM and the temperature response observed in the Wuchiapingian simulations.

Figure 13Dominant surface type for each oxygen level simulation for (a) As-CM and (b) Wu-CM. BLT: broadleaf tree; NLT: needleleaf tree.

3.4 Response of Permian vegetation to pO2

Changes in pO2 and surface temperatures have the potential to impact the terrestrial carbon cycle by altering the competition between the oxidative and photosynthetic metabolic pathways for RuBisCO. Beerling and Berner (2000) simulated significant changes to vegetation productivity in the Permian due to changes in oxygen content. The modelled changes to vegetation in the final 50 years of the Asselian and Wuchiapingian experiments are investigated. Focusing on changes to vegetation across the Permian, the dominant vegetation fractions for As-CM35, As-CM10, Wu-CM35 and Wu-CM10 are shown in Fig. 13. For low pCO2 in the Asselian, increasing pO2 leads to a reduction in the extent of broadleaf trees and greater proliferation of grasses and shrubs. This would be consistent with increases in photorespiration at high pO2. The reverse is true in the Wuchiapingian simulations with increases in the extent of tropical broadleaf forests. It should be noted that the simulation of plant functional types is carefully tuned to present-day vegetation which was likely considerably different in the past. Therefore, caution should be exercised when extrapolating to past vegetation changes.

Figure 14As-CM1035 (left) and Wu-CM1035 (right) anomalies for (a) net primary productivity, (b) total carbon storage and (c) water use efficiency.

Figure 14a shows the change in net primary productivity (NPP) for As-CM1035 and Wu-CM1035. The Asselian simulations shows a large reduction in net primary productivity (NPP) as pO2 is increased (Fig. 14a, 59 Pg C yr−1), while the reverse is true in the Wuchiapingian simulations (+33 Pg C yr−1). At low pCO2, it is expected that competition for RuBisCO will be won out by O2 and therefore that rates of photorespiration should lead to a decline in photosynthesis. This is reflected in the gross primary productivity (GPP, 34 %) and NPP (52 %) response for the Asselian. During the Wuchiapingian, there may be sufficient CO2 that competition is much less sensitive to the pO2 so changes to NPP are much less significant. In fact, NPP is increased by 14 % (GPP +18 %). Tropical water use efficiency is also higher in Wu-CM35 (Fig. 14c), which suggests that water economy of plants could alter to adapt to a higher pO2 (Beerling and Berner2000).

These net primary productivity changes are reflected in the total carbon storage (Fig. 14b) which is lower as pO2 is increased in As-CM (338 Pg C) and higher as pO2 is increased in Wu-CM (+379 Pg C). This is dominated by changes in the tropics (in agreement with Beerling and Berner2000), where broadleaf trees cover more area in Wu-CM35. Cooler terrestrial tropical temperatures, particularly in the warm months (Fig. 4e), reduce the pO2 inhibition of RuBisCO and reduce the rate of respiration by vegetation and soils (Long1991; Beerling and Berner2000).

As these simulations are fully coupled and changes to oxygen content affect temperatures, radiation and precipitation, it is challenging to explore all the possible contributions to differences between these results and the more idealised Beerling and Berner (2000) simulations. However, there is general agreement that changes occur in the signs of the response of NPP and total carbon storage. This supports the conclusions of Beerling and Berner (2000) that high pO2 in the early Permian may have played an important role in the evolution of plants. Note that while the atmosphere and vegetation are coupled in the physical sense, the carbon cycle is not interactive (atmospheric CO2 is fixed), so determining the impacts of atmospheric pO2 on the carbon cycle remains an outstanding problem.

Figure 15The 21 %–10 % O2 anomalies for Poulsen et al. (2015) simulations. (a) Annual-mean, cold-month mean and warm-month mean surface air temperature difference. (b) Change to diurnal cycle and (c) annual-mean precipitation (mm d−1).

Figure 161-D energy balance decomposition analogous to Fig. 7 for the 21 %–10 % O2 Poulsen et al. (2015) simulations.


4 Discussion

Through its impact on atmospheric mass, oxygen content has the capacity to alter the radiative budget of the atmosphere and therefore has implications for the Earth's climate. These simulations suggest that the interactions between radiative and dynamical feedbacks lead to some consistent climatic changes in HadCM3-BL with increasing pO2:

  • reduction in the seasonal cycle in surface air temperature,

  • reduction in Equator-to-pole temperature gradient and

  • reduction in global precipitation.

HadCM3-BL simulates a reduced equilibrium climate sensitivity mainly due to changes in longwave cloud feedbacks. HadGEM3-AO results also support a reduced sensitivity to CO2 content at high pO2. The pre-industrial Holocene results are supported by 1-D radiative convective simulations (Payne et al.2016), 2-D model simulations (Chemke et al.2016) and slab ocean 3-D model simulations of the Archean (Charnay et al.2013). This raises a discrepancy with the Poulsen et al. (2015) study, which simulated a reduction in global-mean surface temperature when increasing oxygen content in the GENESIS model. Figure 15a shows the surface air temperature change between the 10 % and 21 % Cenomanian (100.5–93.9 Ma) simulations from the Poulsen et al. (2015) study. These show a 2.04 C change in the annual mean. To understand the mechanisms behind this, we performed the 1-D energy balance decomposition on the Poulsen et al. (2015) Cenomanian 21 %–10 % model output. The results are shown in Fig. 16. This shows that the cloudy-sky contribution to the temperature change dominates the climate response, contributing 1.45 C. However, the clear-sky contribution is also negative (0.60 C), including both clear-sky emissivity (0.58 C) and clear-sky albedo (0.11 C). This appears to support the argument that tropical cloud feedbacks explain the discrepancy between the Poulsen et al. (2015) simulations and results of 1-D radiative convective models (Goldblatt2016); however, this cannot be the only factor. An increase in pressure broadening of absorption lines would be expected to lead to a positive contribution from the clear-sky emissivity. This suggests that cloud feedbacks alone cannot explain the discrepancy and that the implementation of pressure broadening may play a role in the anomalous Poulsen et al. (2015) response. In addition, changes to the seasonal cycle (Fig. 15a) simulated by Poulsen et al. (2015) are also inconsistent with the HadGEM-AO and HadCM3-BL results, in which all simulations led to a reduced seasonal cycle as pO2 increases. The Poulsen et al. (2015) Cenomanian simulations actually simulates a larger seasonal cycle at high pO2, which is challenging to reconcile with the radiative and physical processes. Note that the Poulsen et al. (2015) simulations were for an earlier Cretaceous period (Cenomanian) than those performed in HadCM3-BL (Maastrichtian); however, the continental configurations and the global-mean temperatures are reasonably similar (22.2 C in HadCM3-BL vs. 20.5 C in Poulsen et al.2015).

The simulations presented in here suggest that perturbations to the wind-driven ocean circulation by increasing atmospheric mass lead to warmer temperatures, particularly at high latitudes. The magnitude of the results varies depending on the precise continental configuration and background climate state. Gyre circulations vary between the pre-industrial and the Maastrichtian and Asselian case studies. Given the importance of the wind-driven ocean circulation response, this suggests that a 3-D representation of ocean circulation is necessary in order to capture the temperature response to atmospheric mass changes. It should be noted, however, that Charnay et al. (2013) simulated higher surface temperatures for the early Earth at high atmospheric mass with a slab ocean model.

The use of 3-D oceans is now widespread in the palaeoclimate community; however, this is not widely used in the exoplanet/early Earth community (e.g. Kilic et al.2017) and for early Earth studies such as the Archean (e.g. Charnay et al.2013). While boundary conditions for these studies are sparse or in some cases non-existent the additional uncertainty associated with using a slab ocean should be considered. AOGCM studies remain the best way to assess the complex coupling between potentially competing radiative and dynamical effects.

One criticism of high oxygen variability in the Phanerozoic is the possibility of runaway fire at high oxygen content (Watson et al.1978). While subsequent experiments have put this in doubt (Wildman et al.2004), fire is undoubtedly a negative feedback on oxygen content. However, the cooling of warmest-month temperatures over tropical and midlatitude continents in Wu-CM35 may provide somewhat of a protective mechanism against runaway fire regimes taking hold. Lightning is a major cause of palaeofire (Scott and Jones1994), so the reduction in convection at high pO2 would also lead to fewer lightning strikes, which would reduce fire initiation. In addition, higher fire risk could have favoured the evolution and spread of more fire-resistant species (Robinson1990).

The simulations of Permian climate (As-CM and Wu-CM) also suggest a strong role for pO2 variability in the terrestrial carbon cycle. However, there are many limitations to the modelling approach employed here. The plant functional types employed here are the same as those for the present day. In particular, C4 photosynthesis likely evolved in the Oligocene (Sage2004), although there is evidence of vegetation which causes C4-like fractionation in the Mississippian, suggesting different vegetation adaptations operating in the past (Jones1994). In addition, angiosperms did not evolve until the Cretaceous, so gymnosperms such as cycads were more widespread in the Permian (Taylor et al.2009). TRIFFID and other dynamic plant models were not developed with these changes in plant types in mind, so simulating past vegetation changes is still a considerable challenge. However, scientific understanding of the role of plants in the climate in the Paleozoic is still immature. While early evidence suggested that late Paleozoic vegetation was unproductive based on analysis of the closest modern relatives, this perspective is increasingly being challenged (Wilson et al.2017). Other approaches such as trait-based methods (Van Bodegom et al.2012; Porada et al.2016) may be able to achieve more insights into the role of pO2 in the Earth system. We also have not accounted for changes to the ocean carbon cycle. A biogeochemical model study suggests that pervasive oceanic anoxia and euxinia only occur below an oxygen level of around 10 % (Ozaki and Tajika2013), which may be below the fire threshold (Belcher et al.2010) and therefore not of relevance to many periods in the Phanerozoic. However, the extent of oceanic anoxic events may be sensitive to atmospheric pO2 (Clarkson et al.2018).

Figure 17(a) Reconstructed CO2 (doubling from Pleistocene values, blue) and O2 content (red) and 95 % confidence intervals (shading) from Royer et al. (2014) Geocarb simulations. (b) GMST reconstructed using Geocarb pCO2 and climate sensitivity values (purple) and the uncertainty in GMST from pCO2 uncertainty (purple shading). GMST reconstructed, accounting for pO2 according to Poulsen et al. (2015) global-mean temperature sensitivities (solid orange) and the uncertainty due to Geocarb pO2 (dashed orange).


Given the small changes in global-mean surface temperature (GMST, 1.5 C maximum) compared to ECS (∼3C), this raises the question of how much pO2 variability contributes to uncertainty in Phanerozoic surface temperature even with such large uncertainties in pO2 reconstructions. Figure 17 shows reconstructed Phanerozoic surface temperatures based on CO2 content and climate sensitivity from the Geocarb model (purple; Royer et al.2014). The uncertainty associated with the 95 % confidence interval in simulated pCO2 is also indicated (purple shading). Analysis of the Poulsen et al. (2015) simulations suggests a global-mean surface temperature reduction of 0.21 C per percentage increase in O2. Accounting for the pO2 simulated by Royer et al. (2014) leads to a mean absolute difference in global-mean surface temperature of 0.80 C and maximum absolute difference of 2.59 C (Fig. 17 orange line). The largest deviations from the Geocarb values occur during the largest deviations from present atmospheric levels of O2 during the Permian. However, pO2 contributes little to the uncertainty in reconstruction of global-mean surface temperature compared to pCO2 (Fig. 17, dashed orange lines), even if the temperature changes simulated by Poulsen et al. (2015) are reasonable. The HadGEM3-AO and HadCM3-BL simulations show even less sensitivity of global-mean surface temperature to pO2 changes, which suggests this is likely an overestimate.

pO2 therefore remains a secondary contribution to climatic variability in the Phanerozoic in agreement with Payne et al. (2016) but is most likely to be important during the Permian. The Artinskian (early Permian, 283.5–290.1 Ma) is associated with a rapid increase in CO2 content from ∼500 to ∼3500 ppmv, which is associated with considerable restructuring of tropical vegetation (Montañez et al.2007). The results in this study suggest that pO2 variability could have modulated the climate and terrestrial vegetation response to this increase in CO2 content. Feulner (2017) suggested that Earth was close to entering a Snowball Earth in the late Carboniferous, when pO2 was higher than today. We hypothesise that the carbon cycle and physical climate feedbacks described in this paper would strongly mitigate against this. If pCO2 and pO2 are intimately linked such that cooler climates tends to increase pO2, this would suggest that pO2 responses have helped to prevent Snowball Earth initiation in the Phanerozoic.

5 Conclusions

The numerical simulations performed in this study reconcile the surface temperature response to oxygen content changes across the hierarchy of model complexity:

  • Under pre-industrial Holocene conditions, increasing atmospheric pO2 leads to an increase in global-mean surface temperature in agreement with 1-D radiative–convective model simulations. This increase is greater in the cold-month mean than in the warm-month mean. The Equator-to-pole temperature gradient is reduced, particularly in the cold-month mean, consistent with a stronger greenhouse effect at high atmospheric pressure.

  • Lower incident surface shortwave radiation leads to a slowdown of the hydrological cycle. Precipitation decreases globally under high pO2, with regional variations.

  • The climate sensitivity is lower at high pO2, particularly in the Maastrichtian. This appears to reconcile the results of the 1-D and 3-D modelling approaches.

  • The climate response simulated by Poulsen et al. (2015) is inconsistent with the radiative changes when considering a 1-D energy balance model decomposition of the surface temperature changes. Tropical cloud feedbacks alone were not sufficient to explain the discrepancy.

  • The climate response to oxygen content variability is state-dependent, so it should be considered on a case-by-case basis. However, the changes are relatively small compared to the role of CO2 in the Phanerozoic (Royer et al.2014).

Code and data availability

Processed model output and analysis scripts will be made available on the NERC data centre. The Met Office Unified Model is available for use under licence. Please see (last access: 27 July 2019) for more information. The ocean bathymetry and land orography reconstructions are ©Getech. Readers who would like advice on how to implement alterations to pO2 in their climate model are encouraged to contact the corresponding authors. UM users can obtain the code changes for these particular versions from the corresponding authors.


The supplement related to this article is available online at:

Author contributions

The study was performed by DCW, conceived by ATA and DCW and refined with input from PJV and the co-authors. DCW and ATA led the preparation of the manuscript, and all co-authors helped in proofreading and checking of the manuscript.

Competing interests

The authors declare that they have no conflict of interest.


This work was carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol – (last access: 27 July 2019). We acknowledge use of the MONSooN system, a collaborative facility supplied under the Joint Weather and Climate Research Programme, a strategic partnership between the Met Office and the Natural Environment Research Council. David C. Wade acknowledges Eric Wolff and David Stevenson for their comments on the PhD thesis of which this paper largely forms a part.

Financial support

This research has been supported by the NERC (grant no. DTP-1502139).

Review statement

This paper was edited by Yannick Donnadieu and reviewed by Jim Kasting and two anonymous referees.


Algeo, T. J. and Ingall, E.: Sedimentary Corg : P ratios, paleocean ventilation, and Phanerozoic atmospheric pO2, Palaeogeogr. Palaeocl., 256, 130–155,, 2007. a, b, c

Armour, K. C., Bitz, C. M., and Roe, G. H.: Time-Varying Climate Sensitivity from Regional Feedbacks, J. Clim., 26, 4518–4534,, 2013. a

Arvidson, R. S., Mackenzie, F. T., and Guidry, M. W.: Geologic history of seawater: A MAGic approach to carbon chemistry and ocean ventilation, Chem. Geol., 362, 287–304,, 2013. a, b, c, d

Arakawa, A. and Lamb, V. R.: Computational Design of the Basic Dynamical Processes of the UCLA General Circulation Model, Methods in Computational Physics: Advances in Research and Applications, 17, 173–265, 1977. a, b

Banerjee, A., Archibald, A. T., Maycock, A. C., Telford, P., Abraham, N. L., Yang, X., Braesicke, P., and Pyle, J. A.: Lightning NOx, a key chemistry-climate interaction: impacts of future climate change and consequences for tropospheric oxidising capacity, Atmos. Chem. Phys., 14, 9871–9881,, 2014. a

Beerling, D., Lake, J., Berner, R., Hickey, L., Taylor, D., and Royer, D.: Carbon isotope evidence implying high O2∕CO2 ratios in the Permo-Carboniferous atmosphere, Geochim. Cosmochim. Ac., 66, 3757–3767,, 2002. a

Beerling, D., Berner, R. A., Mackenzie, F. T., Harfoot, M. B., and Pyle, J. A.: Methane and the CH4-related greenhouse effect over the past 400 million years, Am. J. Sci., 309, 97–113,, 2009. a

Beerling, D. J. and Berner, R. A.: Impact of a Permo-Carboniferous high O2 event on the terrestrial carbon cycle, P. Natl. Acad. Sci. USA, 97, 12428–12432,, 2000. a, b, c, d, e, f, g, h

Beerling, D. J., Fox, A., Stevenson, D. S., and Valdes, P. J.: Enhanced chemistry-climate feedbacks in past greenhouse worlds, P. Natl. Acad. Sci. USA, 108, 9770–9775,, 2011. a

Belcher, C. M. and McElwain, J. C.: Limits for Combustion in Low O2 Redefine Paleoatmospheric Predictions for the Mesozoic, Science, 321, 1197–200, 2008. a

Belcher, C. M., Yearsley, J. M., Hadden, R. M., McElwain, J. C., and Rein, G.: Baseline intrinsic flammability of Earth's ecosystems estimated from paleoatmospheric oxygen over the past 350 million years, P. Natl. Acad. Sci. USA, 107, 22448–22453,, 2010. a, b

Bergman, N. M., Lenton, T. M., and Watson, A. J.: COPSE: A new model of biogeochemical cycling over Phanerozoic time, Am. J. Sci., 304, 397–437,, 2004. a, b, c, d

Berner, R. A.: Phanerozoic atmospheric oxygen: New results using the GEOCARBSULF model, Am. J. Sci., 309, 603–606,, 2009. a, b, c, d

Berner, R. A. and Canfield, D. E.: A new model for atmospheric oxygen over Phanerozoic time, Am. J. Sci., 289, 333–361,, 1989. a

Berner, R. A., Vandenbrooks, J. M., and Ward, P. D.: Evolution, Oxygen and evolution, Science, 316, 557–558,, 2007. a

Bjerknes, J.: Atlantic Air-Sea Interaction, Adv. Geophys., 10, 1–82,, 1964. a

Bony, S., Stevens, B., Frierson, D. M. W., Jakob, C., Kageyama, M., Pincus, R., Shepherd, T. G., Sherwood, S. C., Siebesma, A. P., Sobel, A. H., Watanabe, M., and Webb, M. J.: Clouds, circulation and climate sensitivity, Nat. Geosc., 8, 261–268,, 2015. a

Brand, U., Posenato, R., Came, R., Affek, H., Angiolini, L., Azmy, K., and Farabegoli, E.: The end-Permian mass extinction: A rapid volcanic CO2 and CH4-climatic catastrophe, Chem. Geol., 322/323, 121–144, 2012. a

Breecker, D. O., Sharp, Z. D., and McFadden, L. D.: Atmospheric CO2 concentrations during ancient greenhouse climates were similar to those predicted for AD 2100, P. Natl. Acad. Sci. USA, 107, 576–80,, 2010. a

Broccoli, A. J., Dahl, K. A., and Stouffer, R. J.: Response of the ITCZ to Northern Hemisphere cooling, Geophys. Res. Lett., 33, 1–4,, 2006. a

Budyko, M. I.: The effect of solar radiation variations on the climate of the Earth, Tellus, 21, 611–619,, 1969. a

Caballero, R. and Huber, M.: State-dependent climate sensitivity in past warm climates and its implications for future climate projections, P. Natl. Acad. Sci. USA, 110, 14162–14167,, 2013. a, b

Catling, D. C. and Claire, M. W.: How Earth's atmosphere evolved to an oxic state: A status report, Earth Planet. Sc. Lett., 237, 1–20,, 2005. a

Catling, D. C., Glein, C. R., Zahnle, K. J., and McKay, C. P.: Why O2 Is Required by Complex Life on Habitable Planets and the Concept of Planetary “Oxygenation Time”, Astrobiology, 5, 415–438,, 2005. a

Charnay, B., Forget, F., Wordsworth, R., Leconte, J., Millour, E., Codron, F., and Spiga, A.: Exploring the faint young Sun problem and the possible climates of the Archean Earth with a 3-D GCM, J. Geophys. Res.-Atmos., 118, 10414–10431,, 2013. a, b, c, d, e

Chemke, R., Kaspi, Y., and Halevy, I.: The thermodynamic effect of atmospheric mass on early Earth's temperature, Geophys. Res. Lett., 43, 11414–11422,, 2016. a, b

Clark, D. B., Mercado, L. M., Sitch, S., Jones, C. D., Gedney, N., Best, M. J., Pryor, M., Rooney, G. G., Essery, R. L. H., Blyth, E., Boucher, O., Harding, R. J., Huntingford, C., and Cox, P. M.: The Joint UK Land Environment Simulator (JULES), model description – Part 2: Carbon fluxes and vegetation dynamics, Geosci. Model Dev., 4, 701–722,, 2011. a

Clarkson, M. O., Stirling, C. H., Jenkyns, H. C., Dickson, A. J., Porcelli, D., Moy, C. M., Pogge von Strandmann, P. A. E., Cooke, I. R., and Lenton, T. M.: Uranium isotope evidence for two episodes of deoxygenation during Oceanic Anoxic Event 2, P. Natl. Acad. Sci. USA, 115, 2918–2923,, 2018. a

Cox, P.: A primitive equation, 3-dimensional model of the ocean, GFDL Ocean Group Technical Report No. 1, Tech. Rep., Geophysical Fluid Dynamics Laboratory, Princeton, New Jersey, 1984. a

Cox, P., Huntingford, C., and Harding, R.: A canopy conductance and photosynthesis model for use in a GCM land surface scheme, J. Hydrol., 212/213, 79–94,, 1998. a, b

Cox, P., Betts, R. A., Jones, C., A. Spall, S., and J. Totterdell, I.: Acceleration of global warming due to carbon-cycle feedbacks in a coupled model, Nature, 408, 184–187, 2000. a

Cox, P. M., Betts, R. A., Bunton, C. B., Essery, R. L. H., Rowntree, P. R., and Smith, J.: The impact of new land surface physics on the GCM simulation of climate and climate sensitivity, Clim. Dynam., 15, 183–203,, 1999. a

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

Davies, T., Cullen, M. J. P., Malcolm, A. J., Mawson, M. H., Staniforth, A., White, A. A., and Wood, N.: A new dynamical core for the Met Office's global and regional modelling of the atmosphere, Q. J. Roy. Meteorol. Soc., 131, 1759–1782,, 2005. a

Dwyer, J. G., Biasutti, M., and Sobel, A. H.: Projected Changes in the Seasonal Cycle of Surface Temperature, J. Clim., 25, 6359–6374,, 2012. a

Edwards, C. T., Saltzman, M. R., Royer, D. L., and Fike, D. A.: Oxygenation as a driver of the Great Ordovician Biodiversification Event, Nat. Geosci., 10, 925–929,, 2017. a

Edwards, J. M. and Slingo, A.: Studies with a flexible new radiation code, I: Choosing a configuration for a large-scale model, Q. J. Roy. Meteorol. Soc., 122, 689–719,, 1996. a, b, c

Essery, R. L. H., Best, M. J., Betts, R. A., Cox, P. M., Taylor, C. M., Essery, R. L. H., Best, M. J., Betts, R. A., Cox, P. M., and Taylor, C. M.: Explicit Representation of Subgrid Heterogeneity in a GCM Land Surface Scheme, J. Hydrometeorol., 4, 530–543,<0530:EROSHI>2.0.CO;2, 2003. a

Falkowski, P. G., Katz, M. E., Milligan, A. J., Fennel, K., Cramer, B. S., Aubry, M. P., Berner, R. A., Novacek, M. J., and Zapol, W. M.: The Rise of Oxygen over the Past 205 Million Years and the Evolution of Large Placental Mammals, Science, 309, 1–3, 2005. a, b

Feulner, G.: The faint young Sun problem, Rev. Geophys., 50, RG2006,, 2012. a

Feulner, G.: Formation of most of our coal brought Earth close to global glaciation, P. Natl. Acad. Sci. USA, 114, 11333–11337,, 2017. a

Finney, D. L., Doherty, R. M., Wild, O., Stevenson, D. S., MacKenzie, I. A., and Blyth, A. M.: A projected decrease in lightning under climate change, Nat. Clim. Change, 8, 210–213,, 2018. a

Forster, P. M. and Shine, K. P.: Radiative forcing and temperature trends from stratospheric ozone changes, J. Geophys. Res.-Atmos., 102, 10841–10855,, 1997. a

Franks, P. J., Royer, D. L., Beerling, D. J., Van de Water, P. K., Cantrill, D. J., Barbour, M. M., and Berry, J. A.: New constraints on atmospheric CO2 concentration for the Phanerozoic, Geophys. Res. Lett., 41, 4685–4694, 2014. a, b

Glasspool, I. J. and Scott, A. C.: Phanerozoic concentrations of atmospheric oxygen reconstructed from sedimentary charcoal, Nat. Geosci., 3, 627–630,, 2010. a, b

Goldblatt, C.: Comment on “Long-term climate forcing by atmospheric oxygen concentrations”, Science, 353, 1–2, 2016. a, b

Goldblatt, C., Claire, M. W., Lenton, T. M., Matthews, A. J., Watson, A. J., and Zahnle, K. J.: Nitrogen-enhanced greenhouse warming on early Earth, Nat. Geosci., 2, 891–896,, 2009. a, b, c, d

Gordon, C., Cooper, C., Senior, C. A., Banks, H., Gregory, J. M., Johns, T. C., Mitchell, J. F. B., and Wood, R. A.: The simulation of SST, sea ice extents and ocean heat transports in a version of the Hadley Centre coupled model without flux adjustments, Clim. Dynam., 16, 147–168,, 2000. a

Gregory, J. M., Ingram, W. J., Palmer, M. A., Jones, G. S., Stott, P. A., Thorpe, R. B., Lowe, J. A., Johns, T. C., and Williams, K. D.: A new method for diagnosing radiative forcing and climate sensitivity, Geophys. Res. Lett., 31, L03205,, 2004. a, b, c, d

Hadjinicolaou, P. and Pyle, J. A.: The Impact of Arctic Ozone Depletion on Northern Middle Latitudes: Interannual Variability and Dynamical Control, J. Atmos. Chem., 47, 25–43,, 2004. a

Hagemann, M., Weber, A. P., and Eisenhut, M.: Photorespiration: origins and metabolic integration in interacting compartments, J. Exp. Bot., 67, 2915–2918,, 2016. a

Haigh, J. D. and Pyle, J. A.: Ozone perturbation experiments in a two-dimensional circulation model, Q. J. Roy. Meteorol. Soc., 108, 551–574,, 1982. a

Hansen, K. W. and Wallmann, K.: Cretaceous and Cenozoic evolution of seawater composition, atmospheric O2 and CO2: A model perspective, Am. J. Sci., 303, 94–148,, 2003. a

Harfoot, M. B. J., Beerling, D. J., Lomax, B. H., and Pyle, J. a.: A two-dimensional atmospheric chemistry modeling investigation of Earth's Phanerozoic O3 and near-surface ultraviolet radiation history, J. Geophys. Res.-Atmos., 112, D07308,, 2007. a, b, c

He, T., Pausas, J. G., Belcher, C. M., Schwilk, D. W., and Lamont, B. B.: Fire-adapted traits of Pinus arose in the fiery Cretaceous, New Phytol., 194, 751–759, 2012. a

Heinemann, M.: Warm and sensitive Paleocene-Eocene Climate, Ph.D. thesis, Universität Hamburg, 119 pp., 2009. a

Heinemann, M., Jungclaus, J. H., and Marotzke, J.: Warm Paleocene/Eocene climate as simulated in ECHAM5/MPI-OM, Clim. Past, 5, 785–802,, 2009. a, b, c, d, e

Hewitt, H. T., Copsey, D., Culverwell, I. D., Harris, C. M., Hill, R. S. R., Keen, A. B., McLaren, A. J., and Hunke, E. C.: Design and implementation of the infrastructure of HadGEM3: The next-generation Met Office climate modelling system, Geosci. Model Dev., 4, 223–253,, 2011. a, b

Hunke, E. and Lipscomb, W.: The Los Alamos sea ice model documentation and software user's manual, Version 4.0, LA-CC-06-012., Tech. Rep., Los Alamos National Laboratory, 76 pp., 2008. a

Irvine, P. J., Kravitz, B., Lawrence, M. G., and Muri, H.: An overview of the Earth system science of solar geoengineering, Wires Clim. Change, 7, 815–833,, 2016. a

Johns, T. C., Durman, C. F., Banks, H. T., Roberts, M. J., McLaren, A. J., Ridley, J. K., Senior, C. A., Williams, K. D., Jones, A., Rickard, G. J., Cusack, S., Ingram, W. J., Crucifix, M., Sexton, D. M. H., Joshi, M. M., Dong, B.-W., Spencer, H., Hill, R. S. R., Gregory, J. M., Keen, A. B., Pardaens, A. K., Lowe, J. A., Bodas-Salcedo, A., Stark, S., and Searl, Y.: The New Hadley Centre Climate Model (HadGEM1): Evaluation of Coupled Simulations, J. Clim., 19, 1327–1353,, 2006. a

Jones, T. P.: 13C enriched Lower Carboniferous fossil plants from Donegal, Ireland: carbon isotope constraints on taphonomy, diagenesis and palaeoenvironment, Rev. Palaeobot. Palyno., 81, 53–64,, 1994. a

Jordan, D. B. and Ogren, W. L.: The CO2/O2 specificity of ribulose 1,5-bisphosphate carboxylase/oxygenase, Planta, 161, 308–313,, 1984. a

Kasting, J. F., Liu, S. C., and Donahue, T. M.: Oxygen levels in the prebiological atmosphere, J. Geophys. Res., 84, 3097,, 1979. a

Kiehl, J. T. and Shields, C. A.: Sensitivity of the Palaeocene-Eocene Thermal Maximum climate to cloud properties, Philos. T. R. Soc. A, 371, 20130093,, 2013. a

Kilic, C., Raible, C. C., and Stocker, T. F.: Multiple Climate States of Habitable Exoplanets: The Role of Obliquity and Irradiance, Astrophys. J., 844, 13 pp.,, 2017. a

Lenton, T. M., Crouch, M., Johnson, M., Pires, N., and Dolan, L.: First plants cooled the Ordovician, Nat. Geosci., 5, 86–89,, 2012. a

Li, C., von Storch, J.-S., and Marotzke, J.: Deep-ocean heat uptake and equilibrium climate response, Clim. Dynam., 40, 1071–1086,, 2013. a

Li, J. and Sharma, A.: Evaluation of volcanic aerosol impacts on atmospheric water vapor using CMIP3 and CMIP5 simulations, J. Geophys. Res.-Atmos., 118, 4448–4457,, 2013. a

Long, S. P.: Modification of the response of photosynthetic productivity to rising temperature by atmospheric CO2 concentrations: Has its importance been underestimated?, Plant Cell Environ., 14, 729–739,, 1991. a, b

Lu, W., Ridgwell, A., Thomas, E., Hardisty, D. S., Luo, G., Algeo, T. J., Saltzman, M. R., Gill, B. C., Shen, Y., Ling, H.-F., Edwards, C. T., Whalen, M. T., Zhou, X., Gutchess, K. M., Jin, L., Rickaby, R. E. M., Jenkyns, H. C., Lyons, T. W., Lenton, T. M., Kump, L. R., and Lu, Z.: Late inception of a resiliently oxygenated upper ocean, Science, 361, 174–177,, 2018. a, b

Lunt, D. J., Farnsworth, A., Loptson, C., Foster, G. L., Markwick, P., O'Brien, C. L., Pancost, R. D., Robinson, S. A., and Wrobel, N.: Palaeogeographic controls on climate and proxy interpretation, Clim. Past, 12, 1181–1198,, 2016. a, b, c

Madec, G.: NEMO ocean engine, Note du Pôle de modélisation, Institut Pierre-Simon Laplace (IPSL), France, No 27, ISSN No 1288-1619, 2008. a, b

Meraner, K., Mauritsen, T., and Voigt, A.: Robust increase in equilibrium climate sensitivity under global warming, Geophys. Res. Lett., 40, 5944–5948, 2013. a

Mills, B. J., Belcher, C. M., Lenton, T. M., and Newton, R. J.: A modeling case for high atmospheric oxygen concentrations during the Mesozoic and Cenozoic, Geology, 44, 1023–1026,, 2016. a

Montañez, I. P., Tabor, N. J., Niemeier, D., Dimichele, W. A., Frank, T. D., Fielding, C. R., Isbell, J. L., Birgenheier, L. P., and Rygel, M. C.: CO2-forced climate and vegetation instability during Late Paleozoic deglaciation, Science, 315, 87–91,, 2007. a, b

Nowack, P. J., Luke Abraham, N., Maycock, A. C., Braesicke, P., Gregory, J. M., Joshi, M. M., Osprey, A., and Pyle, J. A.: A large ozone-circulation feedback and its implications for global warming assessments, Nat. Clim. Change, 5, 41–45,, 2014. a, b, c

Ozaki, K. and Tajika, E.: Biogeochemical effects of atmospheric oxygen concentration, phosphorus weathering, and sea-level stand on oceanic redox chemistry: Implications for greenhouse climates, Earth Planet. Sc. Lett., 373, 129–139,, 2013. a

Payne, R. C., Britt, A. V., Chen, H., Kasting, J. F., and Catling, D. C.: The Response of Phanerozoic Surface Temperature to Variations in Atmospheric Oxygen Concentration, J. Geophys. Res.-Atmos., 121, 10089–10096,, 2016. a, b, c, d, e

Pope, V. D., Gallani, M. L., Rowntree, P. R., and Stratton, R. A.: The impact of new physical parametrizations in the Hadley Centre climate model: HadAM3, Clim. Dynam., 16, 123–146,, 2000. a

Porada, P., Lenton, T. M., Pohl, A., Weber, B., Mander, L., Donnadieu, Y., Beer, C., Pöschl, U., and Kleidon, A.: High potential for weathering and climate effects of non-vascular vegetation in the Late Ordovician, Nat. Commun., 7, 12113,, 2016. a

Poulsen, C. J., Tabor, C., and White, J. D.: Long-term climate forcing by atmospheric oxygen concentrations, Science, 348, 1238–41, 2015. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v

Rasulov, B., Hüve, K., Välbe, M., Laisk, A., and Niinemets, U.: Evidence that light, carbon dioxide, and oxygen dependencies of leaf isoprene emission are driven by energy status in hybrid aspen, Plant Physiol., 151, 448–60,, 2009. a

Robinson, J.: Phanerozoic O2 variation, fire, and terrestrial ecology, Palaeogeogr. Palaeocl., 75, 223–240,, 1989. a

Robinson, J. M.: Lignin, land plants, and fungi: Biological evolution affecting Phanerozoic oxygen balance, Geology, 18, 607–610,<0607:LLPAFB>2.3.CO;2, 1990. a, b

Rose, B. E. J. and Ferreira, D.: Ocean Heat Transport and Water Vapor Greenhouse in a Warm Equable Climate: A New Look at the Low Gradient Paradox, J. Clim., 26, 2117–2136,, 2013. a

Royer, D., Berner, R., Montañez, I., and Tabor, N.: CO2 as a primary driver of Phanerozoic climate, GSA today, 14, 4–10, 2004. a

Royer, D. L., Donnadieu, Y., Park, J., Kowalczyk, J., and Godderis, Y.: Error analysis of CO2 and O2 estimates from the long-term geochemical model GEOCARBSULF, Am. J. Sci., 314, 1259–1283,, 2014. a, b, c, d

Saenko, O. A.: On the Climatic Impact of Wind Stress, J. Phys. Oceanogr., 39, 89–106,, 2009. a

Sage, R. F.: The evolution of C4 photosynthesis, New Phytol., 161, 341–370,, 2004. a

Saltzman, M. R., Young, S. A., Kump, L. R., Gill, B. C., Lyons, T. W., and Runnegar, B.: Pulse of atmospheric oxygen during the late Cambrian, P. Natl. Acad. Sci. USA, 108, 3876–81,, 2011. a

Savitzky, A. and Golay, M. J. E.: Smoothing and Differentiation of Data by Simplified Least Squares Procedures, Anal. Chem., 36, 1627–1639,, 1964. a

Scott, A. C. and Glasspool, I. J.: The diversification of Paleozoic fire systems and fluctuations in atmospheric oxygen concentration, P. Natl. Acad. Sci. USA, 103, 10861–10865,, 2006. a, b

Scott, A. C. and Jones, T. P.: The nature and influence of fire in Carboniferous ecosystems, Palaeogeogr. Palaeocl., 106, 91–112,, 1994. a

Sellers, W. D.: A Global Climatic Model Based on the Energy Balance of the Earth-Atmosphere System, J. Appl. Meteorol., 8, 392–400,<0392:AGCMBO>2.0.CO;2, 1969. a

Simmons, A. J. and Strüfing, R.: Numerical forecasts of stratospheric warming events using a model with a hybrid vertical coordinate, Q. J. Roy. Meteorol. Soc., 109, 81–111,, 1983. a

Smith, B. N.: Evolution of C4 photosynthesis in response to changes in carbon and oxygen concentrations in the atmosphere through time, Biosystems, 8, 24–32,, 1976. a

Steinthorsdottir, M. and Pole, M.: Global trends of pCO2 across the Cretaceous–Paleogene boundary supported by the first Southern Hemisphere stomatal proxy-based pCO2 reconstruction, Palaeogeogr. Palaeocl., 464, 143–152,, 2016. a

Stolper, D. A., Bender, M. L., Dreyfus, G. B., Yan, Y., and Higgins, J. A.: A Pleistocene ice core record of atmospheric O2 concentrations, Science, 353, 1427–1430,, 2016. a, b, c

Tappert, R., McKellar, R. C., Wolfe, A. P., Tappert, M. C., Ortega-Blanco, J., and Muehlenbachs, K.: Stable carbon isotopes of C3 plant resins and ambers record changes in atmospheric oxygen since the Triassic, Geochim. Cosmochim. Ac., 121, 240–262,, 2013. a

Taylor, T. N., Taylor, E. L., and Krings, M.: Paleobotany : the biology and evolution of fossil plants, Academic, 1–1252, 2009. a

Timm, S., Wittmiß, M., Gamlien, S., Ewald, R., Florian, A., Frank, M., Wirtz, M., Hell, R., Fernie, A. R., and Bauwe, H.: Mitochondrial Dihydrolipoyl Dehydrogenase Activity Shapes Photosynthesis and Photorespiration of Arabidopsis thaliana, Plant Cell, 27, 1968–84,, 2015. a

Valdes, P. J., Armstrong, E., Badger, M. P. S., Bradshaw, C. D., Bragg, F., Crucifix, M., Davies-Barnard, T., Day, J. J., Farnsworth, A., Gordon, C., Hopcroft, P. O., Kennedy, A. T., Lord, N. S., Lunt, D. J., Marzocchi, A., Parry, L. M., Pope, V., Roberts, W. H. G., Stone, E. J., Tourte, G. J. L., and Williams, J. H. T.: The BRIDGE HadCM3 family of climate models: HadCM3@Bristol v1.0, Geosci. Model Dev., 10, 3715–3743,, 2017.  a, b

Van Bodegom, P. M., Douma, J. C., Witte, J. P. M., Ordoñez, J. C., Bartholomeus, R. P., and Aerts, R.: Going beyond limitations of plant functional types when predicting global ecosystem-atmosphere fluxes: exploring the merits of traits-based approaches, Global Ecol. Biogeogr., 21, 625–636,, 2012. a

Watson, A., Lovelock, J. E., and Margulis, L.: Methanogenesis, fires and the regulation of atmospheric oxygen, Biosystems, 10, 293–298,, 1978. a, b, c

White, A. A. and Bromley, R. A.: Dynamically consistent, quasi-hydrostatic equations for global models with a complete representation of the Coriolis force, Q. J. Roy. Meteorol. Soc., 121, 399–418,, 1995. a

Wildman, R. A., Hickey, L. J., Dickinson, M. B., Berner, R. A., Robinson, J. M., Dietrich, M., Essenhigh, R. H., and Wildman, C. B.: Burning of forest materials under late Paleozoic high atmospheric oxygen levels, Geology, 32, 457–460,, 2004. a, b, c

Williams, K., Copsey, D., Blockley, E., Bodas-Salcedo, A., Calvert, D., Comer, R., Davis, P., Graham, T., Hewitt, H., Hill, R., Hyder, P., Ineson, S., Johns, T. C., Keen, A. B., Lee, R. W., Megann, A., Milton, S. F., Rae, J. G. L., Roberts, M. J., Scaife, A. A., Schiemann, R., Storkey, D., Thorpe, L., Watterson, I. G., Walters, D. N., West, A., Wood, R. A., Woollings, T., and Xavier, P. K.: The Met Office global coupled model 3.0 and 3.1 (GC3. 0 and GC3. 1) configurations, J. Adv. Model. Earth Sy., 10, 357–380, 2018. a

Wilson, J. P., Montañez, I. P., White, J. D., DiMichele, W. A., McElwain, J. C., Poulsen, C. J., and Hren, M. T.: Dynamic Carboniferous tropical forests: new views of plant function and potential for physiological forcing of climate, New Phytol., 215, 1333–1353,, 2017. a

Short summary
The amount of O2 in the atmosphere may have varied from as little as 10 % to as much as 35 % during the last 541 Myr. These changes are large enough to have led to changes in atmospheric mass, which may alter the radiative budget of the atmosphere. We present the first fully 3-D numerical model simulations to investigate the climate impacts of changes in O2 during different climate states. We identify a complex new mechanism causing increases in surface temperature when O2 levels were higher.