The future lifespan of Earth's oxygenated Atmosphere

Introduction

The field of exoplanet astronomy is on the verge of the first detailed characterization of the atmospheres of extrasolar planets, with heightened expectations for the potential detection of spectroscopic indicators of habitability, and possibly life, using the next generation of space and ground-based telescopes1–3. However, a robust framework for remote life detection requires a detailed understanding of the factors that lead to the emergence and long-term maintenance of atmospheric biosignatures, including a full understanding of potential atmospheric composition in the broader context of planetary biogeochemistry4. Although numerous potential atmospheric biosignatures have been proposed5, molecular oxygen ($\mathrm{O}_{2}$) and its photochemical by-product ozone ($\mathrm{O}_{3}$) remain at the forefront of approaches toward remote life detection due to their strong absorption at visible and near-infrared wavelengths, the importance of biology in controlling their atmospheric abundances, and the potential for homogeneity in abundance throughout the atmospheric column4,6. As a result, significant effort has been directed toward understanding atmospheric oxygenation on habitable worlds and the factors that might lead to abiotic oxygen-rich atmospheres7–9.

Earth’s photosynthetic biosphere currently supports large gross fluxes of $\mathrm{O}_{2}$ to the ocean-atmosphere system—roughly $9 \times10^{15}$ mol $\mathrm{O}_{2}$ $\mathrm{y}^{-1}$ (ref. 10)—with the result that the modern atmosphere is 20.5% $\mathrm{O}_{2}$ by volume. However, the presence of oxygenic phototrophs on its own may not be enough to maintain a strongly oxygenated atmosphere. Indeed, it is generally thought that the abundance of $\mathrm{O}_{2}$ in Earth’s atmosphere has been well below that of today for most of Earth’s history11, and that the modern atmospheric $\mathrm{O}_{2}$ abundance developed only after the emergence of land plants12,13, which have evolved to accelerate the geochemical cycles of bioessential elements (most importantly phosphorus) in Earth surface environments14,15. Thus, charting the history of planetary oxygenation after the evolution of oxygenic photosynthesis requires an understanding of how biospheric productivity (and the subsequent burial of reduced materials into sediments) changes in response to evolving environmental conditions.

Previous work on the future lifespan of Earth’s biosphere has focused principally on the links between secular changes in solar luminosity, the stability of the carbonate-silicate geochemical cycle, and the loss of surface water to space16–20. One robust result that emerges from this theoretical framework is a continuous decrease in the abundance of atmospheric $\mathrm{CO}_{2}$ through time, as steadily increasing solar luminosity inexorably drives down the abundance of atmospheric $\mathrm{CO}_{2}$ that is required to balance rates of crustal weathering with volcanic $\mathrm{CO}_{2}$ output16,17,21. However, a number of models suggest that in spite of this thermal buffering mechanism Earth will enter the “moist greenhouse” climate state at some point within the next ~2 billion years (Ga)—a regime in which elevated atmospheric temperatures allow Earth’s stratosphere to become moist, leading to permanent water loss from Earth’s surface16–19,22,23. In addition, this drawdown of atmospheric $\mathrm{CO}_{2}$ may result in a fundamental shift to a photosynthetic biosphere that is $\mathrm{CO}_{2}$-limited—for instance, some models predict that terrestrial $\mathrm{C}_{3}$ plants will cease to be viable at Earth’s surface less than ~500 Myr into Earth’s future16. If true, this might place a long-term physical limit on the ability of photosynthetic biospheres to maintain high levels of atmospheric oxygen, giving rise to a fundamental trade-off between long-term stellar evolution, the geologic carbon cycle, and the intrinsic timescale of atmospheric oxygenation.

Here, we examine the future lifespan of Earth’s oxygenated atmosphere using an Earth system model of biogeochemistry and climate (Figure 1; see Methods). Assessing the lifespan of Earth’s oxygenated atmosphere requires quantitative models that include a complete treatment of oxygen-related biogeochemical processes at the planetary surface over geologic time, as well as mechanistic linkages between surface environments and the planetary interior. Our model builds upon previous similar Earth system models12 that track the coupled carbon, oxygen, phosphorus, and sulphur cycles in the exogenic (i.e., atmosphere, ocean, and crust) system (Figure 1). However, we extend this framework by incorporating a global biogeochemical methane ($\mathrm{CH}_{4}$) cycle, including a range of key biological metabolisms, parameterized atmospheric photochemistry within the $\mathrm{O}_{2}$-$\mathrm{O}_{3}$-$\mathrm{CH}_{4}$ system24,25, and the radiative impact of $\mathrm{CH}_{4}$ on global energy balance26.

Our model is also the first of its type to explicitly evaluate redox exchange between the crust and mantle, allowing a comprehensive evaluation of the planetary processes controlling atmospheric $\mathrm{O}_{2}$ levels over geologic time. The model is designed to capture the major components of Earth system biogeochemistry and climate, but is abstracted enough to allow for large model ensembles that can be run for billions of years of model time. We employ a stochastic approach in which we draw randomly from a wide range of potential values for parameters controlling the major geophysical boundary conditions of the model, including variability (cycle and amplitude) in outgassing rates and erosional forcing (see Methods and Extended Data Figs. 1–3), and initialize the model at 600 Ma. We then repeatedly run the model forward in time until the present day ($n =~400{,}000$), and subsample model runs ($n = 4{,}787$) that roughly reproduce modern atmospheric composition ($p\mathrm{O}_2 = 21\pm1\%$ and $p\mathrm{CO}_2 = 300\pm150$ ppmv), global average surface temperature ($288\pm2$ K), and seawater sulphate ($\mathrm{SO}_4^{2-}$) levels ($\mathrm{SO}_4^{2-} = 25\pm15$ mM). The simulations from this subsampled ensemble are then run forward in time in order to estimate the future lifespan of Earth’s oxygenated atmosphere.

Solar brightening promotes atmospheric deoxygenation

Our default analysis (Figure 2 and Extended Data Figs. 4–8) indicates that atmospheric $\mathrm{O}_2$ levels will decrease substantially in Earth’s long-term future (Figure 2a). Whilst the magnitude and timing of this trajectory show a certain degree of uncertainty, the model-projected decrease in atmospheric $\mathrm{O}_2$ abundance is robust and suggests that the future lifespan of Earth’s oxygenated atmosphere will be less than 1.5 Gyr. Specifically, when we define the future lifespan as the period within which atmospheric $\mathrm{O}_2$ remains above a given threshold value, the lifespan of atmospheric $\mathrm{O}_2$ levels greater than 1% of the present atmospheric level (PAL)—$t_{1\%\mathrm{PAL}}$—is estimated at $1.08\pm0.14$ Gyr ($1\sigma$). Increasing the threshold value to 10% PAL results in a broadly similar lifespan of $\tau_{10\%\mathrm{PAL}} = 1.05\pm0.16$ Gyr ($1\sigma$), a consequence of the extremely short timescale over which atmospheric deoxygenation occurs. Several factors are responsible for the secular decrease in atmospheric $\mathrm{O}_2$ (see Supplementary Information), most notably increasing surface temperatures and atmospheric $\mathrm{CO}_2$ starvation of the terrestrial photosynthetic biosphere driven by solar brightening (Figure 2c and Extended Data Figs. 4–7). Indeed, when we repeat our analysis with a constant solar luminosity we find that no secular $\mathrm{O}_2$ trend remains (Figure 2a and Extended Data Fig. 4). As the geologic drivers of the decreasing $\mathrm{CO}_2$ trend (solar brightening and mantle cooling) are largely undisputed, the first-order features of our projected atmospheric $\mathrm{O}_2$ trajectory should be robust.

Our model further predicts that $\mathrm{CO}_2$ limitation of Earth’s photosynthetic biosphere should ultimately lead to photochemical destabilization of Earth’s oxygenated atmosphere and an abrupt shift in atmospheric $\mathrm{O}_2$ abundance to very low values (Figure 2a). Indeed, our model projects an atmospheric composition for Earth’s long-term future that is similar in many respects to that of the Archaean Earth prior to the so-called ‘Great Oxidation Event’27. In particular, atmospheric $\mathrm{O}_2$ drops to values that are many orders of magnitude below that of the modern atmosphere ($<10^{-6}$ PAL), while the atmospheric abundance of $\mathrm{CH}_4$ increases sharply (Figure 2b). However, a critical difference from the Archaean Earth system is the prediction of much lower atmospheric $\mathrm{CO}_2$ levels in Earth’s long-term future (Figure 2c) and thus the potential for greatly elevated $\mathrm{CH}_4/\mathrm{CO}_2$ ratios (Figure 2d). Such elevated $\mathrm{CH}_4/\mathrm{CO}_2$ ratios are generally predicted to result in the formation of ‘organic haze’ in the atmosphere28–30, with potentially important impacts on climate, atmospheric redox chemistry, and remote life detection. Our results thus suggest that organic haze will be a critical component of Earth’s long-term climate stability into the far future, as well as a potentially promising atmospheric biosignature on Earth-like worlds round late main-sequence stars such as the exoplanet candidate Kepler-452b (orbiting a G2 star with an age of $\sim\$6$\pm$2 Gyr and receiving$\sim$10$% more insolation than Earth receives from the Sun today)31,32.

Controls on the timing of atmospheric deoxygenation

Our default analysis provides the first systematic estimate of the potential lifetime of Earth’s oxygenated atmosphere into the far future.

However, it is important to note that our default model ensemble includes a robust terrestrial biosphere, which may not be a universal feature of habitable Earth-like planets. To examine the impact of the terrestrial biosphere on the longevity of the $\mathrm{O}_2/\mathrm{O}_3$-based biosignatures, we performed an additional ensemble analysis in which we remove the terrestrial biosphere from our default successful runs (Extended Data Figs. 4 and 10). As expected, the absence of land plants results in somewhat lower atmospheric $\mathrm{O}_2$ levels throughout the course of planetary evolution relative to our default analysis. However, we find that the future lifespan of Earth’s remotely-detectable oxygenated atmosphere is similar to that of our default analysis (e.g., $\tau_{1\%\mathrm{PAL}} = 1.04^{+0.17}_{-0.16}$ Gyr; cf. $1.08\pm0.14$ Gyr for the default simulations) (Figure 3a). This is an important result, as it strongly suggests that the presence or absence of terrestrial biosphere exerts a secondary control on the timescale of the planetary deoxygenation and long-term atmospheric redox state. Although it is difficult to explicitly model uncertainties associated with biological evolution and adaptation in photosynthetic lineages, our statistical approach explores a wide range of parameter space for key biological response functions (see Methods and Extended Data Figs. 1–3), with results suggesting that uncertainty in the parameters governing the response of the biosphere plays a secondary role in controlling the future lifespan of oxygenated atmospheric conditions on Earth (Extended Data Fig. 9).

Our analysis also delineates an important relationship between the future lifespan of oxygenated atmospheric conditions and the long-term average flux of reducing power from the mantle to the surface system, $\Phi_{1\%}(\mathrm{Red})$ (Figure 3b). In particular, our model suggests that larger $\Phi_{1\%}(\mathrm{Red})$ values tend to yield shorter oxygenated lifespans. Improvements in quantitative estimates of subduction and mantle degassing over time would reduce uncertainty in estimates of the future lifespan of Earth’s oxygenated atmosphere, and will be critical for extending our model framework to very different tectonic regimes. In addition, our analysis suggests that oceanic redox chemistry—which ultimately controls oceanic P levels and the redox state and subduction fluxes of major C, S, and Fe species in the exogenic system—would be critical for the long-term trajectory of atmospheric redox state on habitable Earth-like planets (Figure 3b and Extended Data Fig. 4). The lifespans of oxygen-rich atmospheric states on Earth-like planets are thus likely to be modulated by mantle redox

state, rates of surface-interior exchange, and their evolution through time.

It is important to note that there are multiple biogeochemical and climate processes that are not considered in our model that may play a role in constraining the future lifespan of Earth’s biosphere and the timing/mode of transition to more reducing atmospheric conditions. In particular, “reverse weathering” (the formation of authigenic silicates in marine sediments, resulting in net $\mathrm{CO}_2$ release to the ocean-atmosphere system)33 could potentially extend the lifespan of oxygenated atmospheric conditions under certain scenarios by prolonging the timescale over which atmospheric $\mathrm{CO}_2$ is above the levels expected to result in $\mathrm{CO}_2$ limitation of the photosynthetic biosphere34. On the other hand, this process would commensurately elevate global temperature in the face of increasing solar luminosity, with a corresponding potential for limitation of biospheric $\mathrm{O}_2$ fluxes (Extended Data Fig. 1). The net effect on the future lifespan of Earth’s oxygenated atmosphere is uncertain, but is likely of secondary importance to the overall trend presented here (see also Extended Data Fig. 5 for the sensitivity to terrestrial weatherability). In addition, we hypothesize that haze-induced climate cooling28 could potentially act as a brake on the overall magnitude of atmospheric deoxygenation, or result in the inception of oxygenation/deoxygenation cycles during Earth’s terminal habitability. Both of these topics represent important areas of future research.

Implications for the search for life on exoplanets

Our results also have important implications for the search for life on Earth-like planets beyond our solar system (e.g., habitable planets with abundant liquid water at the surface, exposed silicate crust, and a biosphere with oxygenic photosynthesis). From the perspective of planetary evolution, our results imply that atmospheric oxygenation is not a permanent condition on habitable worlds hosting oxygenic photosynthesis, and that only a fraction of Earth’s history will be characterized by robustly detectable levels of atmospheric $\mathrm{O}_2$ (Figure 4). For example, it is estimated that a G-dwarf stellar flux $\sim20\%$ greater than that of the present Sun ($S/S_0 \sim1.2$) is required to achieve a water loss timescale of $\sim1$ Gyr in the moist greenhouse regime35. Assuming it will take $\sim2$ Gyr to achieve $S/S_0 \sim1.2$ (ref.36), together with the assumption of surface liquid water by 4.4 Ga37, yields an estimated habitable lifetime on Earth of $\sim7.4$ Gyr. Our results indicate that direct detection of $\mathrm{O}_2$ at visible wavelengths would be challenging for all but $\sim1.5$-2.0 Gyr of this history — or roughly 20-30% of Earth’s lifetime as a habitable world38. On the other hand, observations of $\mathrm{O}_3$ at UV wavelengths39 could potentially significantly extend this observability timescale (Figure 4). In any case, our model suggests that a fundamental feature of the carbonate-silicate cycle, considered central to the long-term maintenance of habitable climate conditions on terrestrial planets, may ultimately act to drive atmospheric composition toward an anoxic, hazy state similar to that of the Archaean Earth but with much lower atmospheric $\mathrm{CO}_2$. These predictions emphasize the need to better understand the factors regulating redox balance between the planetary interior and surface environments40–44, and provide additional motivation for implementing a wide wavelength range in future exoplanet characterization missions as a contingency against “false negatives” for remote life detection3,38.

In sum, our stochastic analysis suggests that the eventual deoxygenation of Earth’s atmosphere is a robust outcome of increasing solar luminosity, irrespective of large uncertainties in geophysical/biological boundary conditions, and that the collapse of Earth’s oxygenated atmosphere is likely to occur before the inception of moist greenhouse conditions in Earth’s climate system although its precise timing is modulated by net input flux of reducing power from the mantle to Earth’s surface system. In addition, our model projections of atmospheric chemistry implicate organic haze as a critical factor in regulating Earth’s long-term future climate. Evaluating these factors using more sophisticated biogeochemistry and climate models is an important task for future research. Nevertheless, our analysis clearly identifies the first-order response of atmospheric O$_2$ to declining CO$_2$ and biospheric productivity, opening a new perspective on the future evolution of Earth’s climate and atmospheric redox chemistry.

Corresponding author:

Correspondence to: Kazumi Ozaki, E-mail: [email protected]

Acknowledgements: The authors thank E. Tajika, Y. Sekine, S. Kadoya, Y. Watanabe, and E. Schwieterman for helpful discussions. K.O. acknowledges support from the NASA Postdoctoral Program at the NASA Astrobiology Program, administered by Universities Space Research Association under contact with NASA. This work was supported by JSPS KAKENHI Grant Number JP20K04066. C.T.R. acknowledges support from the NASA Astrobiology Institute (Grant 13-13NAI7_2-0027) and the NASA Nexus for Exoplanet System Science (NExSS) (Grant 80NSSC19KO461).

Author Contributions: K.O. and C.T.R. designed the study. K.O. constructed the model and performed experiments. K.O. and C.T.R. analysed the results and wrote the paper. All authors discussed and contributed intellectually to the interpretation of the results.

Competing interests: The authors declare no competing interests.

Methods

Model description. Our model is based on the Carbon Oxygen Phosphorus Sulphur Evolution (COPSE) model12,48, the default version of which is described in full elsewhere12,48. The model has been extensively tested and validated against geologic records from the most recent $\sim500$ million years of Earth’s history and can predict the coupled histories and controls on atmospheric $\mathrm{CO_2}$, $\mathrm{O_2}$ and ocean composition ($\mathrm{PO_4^{3-}}$, $\mathrm{NO_3^-}$, $\mathrm{SO_4^{2-}}$) on geologic timescales. The COPSE model was recently revised by ref.[51], but the version used here is based mainly on a previous version of the model described in ref.12,15, because we have constructed the model with the aim of keeping the design as simple as possible. We do not anticipate our overarching conclusions to be altered by the differences. Our aim here is to produce a quantitative model that comprehensively encapsulates a range of key biogeochemical processes in oxic/anoxic biospheres and makes resulting predictions for the future lifespan of Earth’s oxygenated atmosphere. Specifically, we made the following additions to the basic structure of the model: (1) metabolic $\mathrm{CH_4}$ production/consumption processes in the ocean, (2) parameterized atmospheric photochemistry within the $\mathrm{O_2}$-$\mathrm{O_3}$-$\mathrm{CH_4}$ system24,25, (3) the radiative impact of $\mathrm{CH_4}$ on global energy balance26, (4) a dependence of marine phytoplankton growth rate on surface temperature24 and marine inorganic carbon availability, (5) exchange of C, S, and Fe between the exogenic (ocean-atmosphere-crust) system and the mantle, and (6) a full redox ($\mathrm{O_2}$) budget in the exogenic system.

Primary production and decomposition. In the COPSE model, the activity levels of terrestrial and marine ecosystems are evaluated separately. The global net primary production (NPP), $J_{\mathrm{NPP}}$ (in terms of organic carbon), is given by:

where $J_{\mathrm{NPP}}^{\mathrm{ocn}}$ and $J_{\mathrm{NPP}}^{\mathrm{land}}$ denote marine and terrestrial biospheric productivity, respectively. $J_{\mathrm{NPP}}^{\mathrm{ocn}}$ is assessed as follows:

where * denotes the modern value of $\sim45$ Gt C yr$^{-1}$, or 3750 Tmol C yr$^{-1}$ (ref.10), and $f_n$ and $f_T^{\mathrm{ocn}}$ denote the effect of nutrient availability and surface temperature on oceanic productivity, respectively. A Liebig’s minimum law is adopted for $f_n$:

where N, P, and $M_{\mathrm{CO_2+HCO_3^-}}$ denote reservoir sizes of seawater nitrate, phosphate and inorganic carbon ($\mathrm{CO_{2aq}}$ and $\mathrm{HCO_3^-}$) available for phytoplankton growth, and $r_{\mathrm{C/P}}$ (= 117) and $r_{\mathrm{N/P}}$ (= 16) denote the “Redfield ratios” of photosynthetic biomass. In the original COPSE model the growth rate of marine phytoplankton is limited by the availability of nutrients (N or P). However, even in the modern ocean, $\mathrm{CO_2}$-limitation can reduce growth rates among some phytoplankton taxa[52, 53]. Given a secular decrease in atmospheric $\mathrm{CO_2}$ levels (Figure 2c), marine phytoplankton could be limited by the availability of dissolved inorganic carbon ($\mathrm{CO_{2aq}}$ and $\mathrm{HCO_3^-}$) in the remote future. The dissolved concentrations of $\mathrm{CO_{2aq}}$ and $\mathrm{HCO_3^-}$ are calculated in the model by solving the seawater carbonate system (see Supplementary Table 3). For the temperature response of marine phytoplankton growth, we adopt the following general function based on 24:

where $k$ denotes a normalizing factor (= 1.2748), and $T_{\mathrm{opt}}$ and $\Delta T_{\mathrm{opt}}$ are assumed to be 301 K and 301 K$^3$, respectively. The maximum temperature for phytoplankton growth is controlled by $T_{\mathrm{ref}}$, whose uncertainty is explored by the statistical analysis (see below and Extended Data Figs. 1 and 3).

Almost all organic matter produced in the surface ocean is decomposed before burial. Fractions of organic matter decomposed by aerobic respiration (combined with denitrification/nitrification and microbial sulphate reduction/sulphide oxidation) and methanogenesis are expressed as follows:

where O denotes the oxygen reservoir size in the combined ocean-atmosphere system and $K_{\mathrm{O}_2}=13.6\times10^{18}$ mol. The total organic carbon decomposition rate is given by $J_{\mathrm{NPP}}^{\mathrm{ocn}}-J_{\mathrm{org}}^{\mathrm{b,ocn}}\equiv J_{\mathrm{org}}^{\mathrm{r,ocn}}$, where $J_{\mathrm{org}}^{\mathrm{b,ocn}}$ denotes the burial rate of organic carbon in marine sediments that is proportional to the square of the normalized primary production as in the original COPSE model48. The global methane generation rate via methanogenesis ($2\mathrm{CH}_2\mathrm{O}\rightarrow\mathrm{CH}_4+\mathrm{CO}_2$) is given by $0.5\times f_{\mathrm{CH}_4}\times J_{\mathrm{org}}^{\mathrm{r,ocn}}$. A fraction of this methane is consumed by aerobic/anaerobic methane oxidation processes (aerobic methanotrophy and AOM), given by:

where $K_{\mathrm{O}_2}^{\prime}=0.273\times10^{18}$ mol ([54]). Given the stoichiometric relationships, the net production rate of CH$_4$, CO$_2$, and O$_2$ in the ocean interior are given by:

For simplicity, the uncertainty associated with partitioning of $\mathrm{CH}_4$ consumption between aerobic methanotrophy and AOM is accounted for implicitly by prescribing $K_{\mathrm{O}_2}'$.

The NPP of Earth’s terrestrial biosphere is scaled by a factor $V$ which represents global land biomass normalized by the modern value48:

where the NPP of the present terrestrial biosphere, $J^{\mathrm{land},*}_{\mathrm{NPP}}$, is $\sim$60 Gt C yr$^{-1}$, or 5000 Tmol C yr$^{-1}$ (10), and $V$ is calculated, as follows:

where $V_{\mathrm{npp}}$ denotes the effects of atmospheric $\mathrm{O}_2$ and $\mathrm{CO}_2$ concentrations and surface temperature on terrestrial biomass as in the original model (Extended Data Fig. 1a,b)48, $f_{\mathrm{fire}}$ denotes the effect of the fire feedback12,48, and $f_{\mathrm{P}}$ denotes a forcing factor representing selective-P weathering by vascular plant activity12,15. Our model introduces an additional factor, $f_{\mathrm{UV}}$, in order to capture the effect of UV on the terrestrial biosphere as a function of atmospheric $\mathrm{O}_2$ levels:

where $p\mathrm{O}_2(\mathrm{PAL})$ denotes atmospheric $\mathrm{O}_2$ levels in terms of PAL, and $c_{\mathrm{UV}}$ is treated as one of the model parameters explored during our stochastic analysis (see below). The minimum $\mathrm{CO}_2$ concentration for the terrestrial plant activity is also treated as one of the uncertainties, because the ecological $\mathrm{CO}_2$ threshold is not precisely known and its absolute value is likely to change with global mean temperature via its influence on photosynthesizer respiration rate[55, 56].

The $\mathrm{CH}_4$ flux from the terrestrial biosphere is assumed to be proportional to the burial rate of organic matter on land, as follows:

where $f^{\mathrm{land}}_{T}$ denotes the temperature effect on the $\mathrm{CH}_4$ flux[57]:

where $T_C$ denotes the surface temperature in terms of °C. On the other hand, the net production rate of $\mathrm{CH}_4$, $\mathrm{CO}_2$, and $\mathrm{O}_2$ in the terrestrial biosphere can be written, as follows:

where $J^{\mathrm{r,land}}$ (=$J_{\mathrm{NPP}}^{\mathrm{land}}-J_{\mathrm{org}}^{\mathrm{b,land}}$) denotes the decomposition rate of organic matter in the terrestrial biosphere, and the fraction of methane oxidation is calculated by equation (7) above. A fraction of organic matter decomposed by methanogenesis, $g_{\mathrm{CH}_4}$, is determined from equations (14) and (16). Then $g_{\mathrm{O}_2}$ is determined from $g_{\mathrm{O}_2}=1-g_{\mathrm{CH}_4}$.

Photochemistry. The revised model includes parameterized photochemistry that allows calculation of the concentrations of atmospheric $\mathrm{O}_2$ and $\mathrm{CH}_4$. We incorporate a scheme for parameterized $\mathrm{O}_2$-$\mathrm{O}_3$-$\mathrm{CH}_4$ photochemistry based on the previous study of a 1-D photochemical model24. Specifically, the $\mathrm{CH}_4$ oxidation rate $k_{\mathrm{CH}_4\mathrm{ox}}$ (in mol$^{-1}$ y$^{-1}$) is expressed as a polynomial function of the reservoir sizes of $\mathrm{O}_2$ and $\mathrm{CH}_4$ (O and M)25.

where $\alpha_i^j$ is fitting coefficients for a given atmospheric $\mathrm{CH}_4$ levels and $\varphi_{\mathrm{O}_2}$ is given by $\log p_{\mathrm{O}_2}(\mathrm{bar})$ (Supplementary Table 4). The oxidation rate is evaluated based on Figure 3 of ref. 24 showing the oxidation rate as a function of $p_{\mathrm{O}_2}$ and $p_{\mathrm{CH}_4}$. We took the relationship between $k_{\mathrm{CH}_4\mathrm{ox}}$ and $p_{\mathrm{O}_2}$ for $p_{\mathrm{CH}_4}$ of $10^{-6}$, $10^{-5}$, $10^{-4}$, $10^{-3}$, $2\times10^{-3}$ bar, and $k_{\mathrm{CH}_4\mathrm{ox}}$ is calculated as a function of $p_{\mathrm{O}_2}$ and $p_{\mathrm{CH}_4}$ with a log-linear interpolation method.

For the atmospheres examined in this study, $\mathrm{CH}_4$ and $\mathrm{H}_2\mathrm{O}$ are the dominant hydrogen-bearing species in the stratosphere, allowing the hydrogen escape flux, $J_{\mathrm{esc}}$, to be calculated as:

where the proportional coefficient, $f_{\mathrm{esc}}$, is $6.68\times10^{15}$ mol yr$^{-1}$ and $f_T(\mathrm{H}_2)$ denotes the total hydrogen molecule mixing ratio:

Therefore, fluxes of hydrogen escape via $\mathrm{CH}_4$ and $\mathrm{H}_2\mathrm{O}$ are written in terms of the abundances in the stratosphere as:

We fit the empirical relationship between the stratospheric water vapour concentration and global average surface temperature obtained from a suite of 3-D GCM simulations of ref. 35 with a 3rd order polynomial in surface temperature:

where $f_i$ denotes the fitting coefficients ($f_0 = 198.34$, $f_1 = -1.6856$, $f_2 = 4.147\times10^{-3}$, and $f_3 = -2.7947\times10^{-6}$) (Supplementary Fig. 2).

Subduction flux. Despite its crucial role in controlling the redox budget in the exogenic system, a process not typically represented in large-scale biogeochemical models is the subduction of redox sensitive elements (C, S, and Fe species) into the mantle. The exchange flux between Earth surface and the mantle influences not only the sizes of the surficial reservoirs but the oxygen balance of the atmosphere on the long timescales considered here. For the modern Earth the subduction flux of organic C is estimated at $<0.2$-0.7 Tmol $\mathrm{O_2}$ equiv. yr$^{-1}$ (40,42), and that of sedimentary pyrite would also be small ($<1$ Tmol S yr$^{-1}$)[58] [42]. On the other hand, the rate of subduction of ferric iron could be 3 Tmol Fe yr$^{-1}$ (0.75 Tmol $\mathrm{O_2}$ equiv. yr$^{-1}$)40, or as high as 7 Tmol Fe yr$^{-1}$ (1.75 Tmol $\mathrm{O_2}$ equiv. yr$^{-1}$)[59] to 19 Tmol Fe yr$^{-1}$ (4.8 Tmol $\mathrm{O_2}$ equiv. yr$^{-1}$)42. However, there is a large uncertainty in the major oxidants ($\mathrm{O_2}$ vs. $\mathrm{SO_4^{2-}}$) for $\mathrm{Fe^{2+}}$ during aqueous oxidation of seafloor[60]. These redox-dependent subduction fluxes probably vary in intensity in response to changes in tectonic activity and oceanic redox state. In this study, the subduction fluxes of C, S, and Fe are parameterized as follows:

where * denotes a reference value, which is treated as a free parameter. The rate of subduction is scaled to the tectonic activity represented by a forcing $f_{\mathrm{G}}$ (D in the original COPSE model). The subduction fluxes of organic carbon and carbonate carbon are also proportional to their burial fluxes. We also assume that the subduction flux of Fe-oxides is a function of the degree of anoxia (DOA) which represents oceanic anoxicity. The rate of subduction of ferric iron is thought to be sensitive to the redox landscape of the global ocean (e.g., sulfidic or Fe-rich anoxia), but this is outside the scope of the present work. Even if subduction of iron oxides is assumed to continue in an anoxic, Fe-rich ocean, this would not change our overarching conclusion that the future lifespan of the oxygenated atmosphere is less than 1.5 Gyr because the subduction of ferric iron represents the loss of oxidizing power from the surface system. Given the large uncertainty in the magnitude of subduction fluxes, we explore a wide range of reference values in our stochastic analysis (Extended Data Fig. 3).

Global O$_2$ budget. In this study, redox balance in the model is tracked in units scaled to the oxidizing power of O$_2$. The model explicitly accounts for fluxes of oxidants and reductants into and out of the surface system (Supplementary Fig. 3). The global redox budget (GRB) for the ocean-atmosphere system (OAS) is written as follows:

where $J_{\mathrm{Hesc}}$ denotes the irreversible hydrogen escape to space, and second and third terms represent the influence of organic C and pyrite S subcycles on the GRB. $J_{\mathrm{H_2S}}^{\mathrm{d}}$ and $J_{\mathrm{red}}^{\mathrm{d}}$ are the mantle degassing of H$_2$S and other reducing gases (e.g., H$_2$). $J_{\mathrm{Fe-ox}}^{\mathrm{s}}$ is the subduction flux of ferrous iron. In this study, we also revised the baseline COPSE model to capture the O$_2$ budget in the exogenic (ocean-atmosphere-crust) system on geological timescales, according to:

The primary sink of reducing power from the exogenic system is the irreversible escape of hydrogen to space ($J_{\mathrm{Hesc}}$) and the subduction of organic C and pyrite S to the mantle ($J_{\mathrm{org}}^{\mathrm{s}}$ and $J_{\mathrm{py}}^{\mathrm{s}}$), while the primary source of reducing power is outgassing of reduced gases from the Earth’s interior ($J_{\mathrm{H_2S}}^{\mathrm{d}}$ and $J_{\mathrm{red}}^{\mathrm{d}}$). The subduction of ferric iron ($J_{\mathrm{Fe-ox}}^{\mathrm{s}}$) is also a loss of oxidizing power (i.e., source of reducing power). The net input flux of reducing power from the mantle to the exogenic system is defined as:

In Figure 3b of the main text, we assess the impact of long-term $\Phi(\mathrm{Red})$ on the future lifespan of oxygenated atmosphere. To do this, we average $\Phi(\mathrm{Red})$ over the future lifespan:

where $\tau_i$ denotes the future lifespan for the different values of O$_2$ threshold i (in PAL). Existing empirical estimates of the mantle degassing of reducing gases via submarine volcanism, $0.7 \pm0.3$ Tmol O$_2$ equiv. yr$^{-1}$ ((ref. 60)), the limited O$_2$ consumption through seafloor Fe$^{2+}$ oxidation, 0.14 Tmol O$_2$ equiv. yr$^{-1}$ ((ref. 60)), and the subduction flux of organic C of <0.2-0.7 Tmol O$_2$ equiv. yr $^{1}$, yield a present net input flux of reducing power from the mantle to the surface system, $(\mathrm{Red})^*$, of $>-0.16$ Tmol $\mathrm{O}_2$ equiv. yr$^{-1}$. A resampling analysis of $\Phi(\mathrm{Red})^*$ from our successful model simulations predicts a slightly lower value of $-0.24 \pm0.27$ Tmol $\mathrm{O}_2$ equiv. yr$^{-1}$ (Supplementary Fig. 4).

Climate model. Global mean surface temperature $T_{\mathrm{s}}$ (in K) is estimated with a zero-dimensional energy balance model (0-D EBM) in which the absorption of solar radiation is balanced by the heat loss by outgoing longwave (i.e., infrared) radiation (OLR) to space ($J_{\mathrm{OLR}}$):

where $\alpha$ and $S_{\odot}$ denote the top-of-atmosphere (TOA) albedo and the value of the solar constant at the time under consideration, respectively. For $J_{\mathrm{OLR}}$ and $\alpha$, we adopted the parametric expressions proposed by ref.61, with the following variables: $p = p_{\mathrm{CO}_2}$ (in bar); $\phi= \ln(p/3.3 \times10^{-4})$; $\mu= \cos(z)$, where $z$ is the zenith angle; and surface albedo $\alpha_{\mathrm{s}}$. $\alpha_{\mathrm{s}}$ is written as a function of $T_{\mathrm{s}}$:

where it is assumed that Earth is globally ice-covered for $T_{\mathrm{s}} < T_{\mathrm{ice}} (= 263$ K) with an albedo $\alpha_{\mathrm{ice}} (= 0.65$), and that Earth is characterized by a relatively low albedo $\alpha_0 (= 0.35)$ for $T_{\mathrm{s}} > T_0 (= 273$ K). For $T_{\mathrm{ice}} < T_{\mathrm{s}} < T_0$, the albedo is interpolated with a quadratic interpolation method[62]. $RF_{\mathrm{CH}_4}$ represents the radiative forcing by methane, which is quantified based on ref.26. Recent 1-D modelling results suggest that relatively small increases in solar luminosity could conceivably trigger a runaway greenhouse on Earth[63, 64]. However, 3-D climate models indicate that Earth’s climate is likely to be somewhat more stable than this as a consequence of expansion of unsaturated regions in the upper atmosphere attendant to changes in Hadley circulation that occur at higher solar fluxes[65] [23]. Given the quantitative consistency between the above 0-D EBM and recent 3-D climate models35 (Supplementary Fig. 1), we consider our utilization of 0-D EBM a defensible approximation of planetary climate for implementation on $10^6$-$10^9$-year timescales. Nevertheless, it is important to point out that future projections of Earth system behavior will likely be modulated by uncertain cloud feedbacks and continental configurations. As a result, in our MC analysis we explore uncertainties in OLR and TOA (Extended Data Fig. 3).

Stochastic analysis. The computational efficiency of the model enables us to examine the future lifespan of the oxygenated atmosphere via a stochastic analysis, in which we randomly sample from a prior distribution for poorly constrained parameters and quantify whole-model uncertainty (Extended Data Fig. 3). We draw poorly constrained initial values for crustal C and S reservoir sizes at 600 Ma (million years ago) over a wide range of values. We also vary the reference values of fluxes between the exogenic system and the mantle (subduction and mantle degassing) to capture uncertainties in model parameterisation. After running the model forwards in time until the present, we only sample model runs that allow for acceptable atmospheric composition ($\mathrm{pO}_2 = 21 \pm1\%$ and $\mathrm{pCO}_2 = 300 \pm150$ ppmv), global surface temperature ($288 \pm2$K) and seawater sulphate ($\mathrm{SO}_4^{2-}$) levels ($\mathrm{SO}_4^{2-} = 25 \pm15$ mM). The model runs that fail to obtain reasonable values have been abandoned. For our default stochastic analysis, we performed nearly 400,000 simulations to span the whole parameter space, and the resultant successful runs were $n = 4,787$ which were continued to be run for estimating the future lifespan of the oxygenated atmosphere.

For the future evolution of the degassing factor ($f_\mathrm{G}$) and erosion factor ($f_\mathrm{R}$; $U$ in the original model), we assumed following functions:

where the first term on the right-hand side denotes the periodic changes with an amplitude $a$ and a period of $\tau$ (in Myr), and the second term represents the secular decrease in the activity of tectonics[66]. These relationships have been adopted for examining to the geological past, and we consider that this assumption is a reasonable first-order approximation as there is no a priori reason for thinking that this style of decay will not continue into Earth’s far future. Nevertheless, we also emphasize that this is very under-studied and it is possible that non-linearity in tectonic regime will be important to consider for Earth and other habitable rocky planets in future work.

We draw model parameters from wide distributions that represent the uncertainty in their values (Extended Data Fig. 3). For example, the uncertainty in the biological responses is assessed by changing biological parameters, such as FERT. The strength of terrestrial weathering feedback is also a topic of debate[67, 68, 69], and we explore its uncertainty by changing the apparent activation energy (ACT) and a runoff factor (RUNsil)[70]. The amplification of terrestrial weatherability by land plants is also an important uncertainty[71, 72], which we explore here by changing LIFE[70]. To account for the uncertainty in the value of $c_\mathrm{UV}$, we draw it from a log-uniform distribution between $10^{-2.5}$ PAL and $10^{-1.5}$ PAL (Extended Data Fig. 1c). At the higher end of this range, terrestrial plant activity is suppressed when atmospheric $\mathrm{O}_2$ is $<10\%$ PAL. At the lower end, the threshold is lowered to $\sim1\%$PAL. We consider these two end-member cases to cover the full range of uncertainty of UV impacts on the terrestrial ecosystems. We also explore uncertainty in the strength of redox-dependent P recycling in marine sediments by changing $\mathrm{C}_{\mathrm{org}}/\mathrm{P}_{\mathrm{org}}$ ratios for anoxic sediments.

The sampling procedure is run for long enough to obtain stationary distributions for posterior probability densities. However, it is possible that our model predictions may underestimate the inherent variability in the future Earth’s evolution. For instance, the possibility of random catastrophe, such as the asteroid impacts, could also introduce additional uncertainty in the estimate of future lifespan of Earth’s oxygenated atmosphere. We also note that the uncertainty in the longevity of human impacts on Earth’s future environment is beyond the scope of this study.

Data availability. The numerical data for Figure 3 in the main text and Extended Data Figures 4 and 5 are presented in the Source data. The big data obtained by the statistical analysis is available to download at https://doi.org/10.6084/m9.figshare.13487487.v1.

Code availability. Our Fortran source code is available at https://github.com/kazumi-ozaki/lifespan.