# Quantification of ocean heat uptake from changes in atmospheric O2 and CO2 composition

L. Resplandy, R. F. Keeling, Y. Eddebbar, M. Brooks, R. Wang, L. Bopp, M.

C. Long, J. P. Dunne, W. Koeve, A. Oschlies

► To cite this version:

L. Resplandy, R. F. Keeling, Y. Eddebbar, M. Brooks, R. Wang, et al. Quantification of ocean heat uptake from changes in atmospheric O2 and CO2 composition. Scientific Reports, 2019, 9, 10.1038/s41598-019-56490-z. insu-03727000

HAL Id: insu-03727000

https://insu.hal.science/insu-03727000

Submitted on 28 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
OPEN

# Quantification of ocean heat uptake from changes in atmospheric O₂ and CO₂ composition

L. Resplandy<fnref n="1" />*, R. F. Keeling², Y. Eddebbar², M. Brooks², R. Wang³, L. Bopp⁴, M. C. Long⁵,
J. P. Dunne⁶, W. Koeve⁷ & A. Oschlies⁷

The ocean is the main source of thermal inertia in the climate system. Ocean heat uptake during recent decades has been quantified using ocean temperature measurements. However, these estimates all use the same imperfect ocean dataset and share additional uncertainty due to sparse coverage, especially before 2007. Here, we provide an independent estimate by using measurements of atmospheric oxygen (O₂) and carbon dioxide (CO₂) – levels of which increase as the ocean warms and releases gases – as a whole ocean thermometer. We show that the ocean gained 1.29 ± 0.79 × 10²² Joules of heat per year between 1991 and 2016, equivalent to a planetary energy imbalance of 0.80 ± 0.49 W watts per square metre of Earth's surface. We also find that the ocean-warming effect that led to the outgassing of O₂ and CO₂ can be isolated from the direct effects of anthropogenic emissions and CO₂ sinks. Our result – which relies on high-precision O₂ atmospheric measurements dating back to 1991 – leverages an integrative Earth system approach and provides much needed independent confirmation of heat uptake estimated from ocean data.

A fundamental measure of global warming is the heat uptake by the ocean, which represents more than 90% of the excess energy gained by the Earth<fnref n="1" />. This ocean warming has been quantified using hydrographic temperature measurements, including data from the Argo float program, which expanded coverage after 2007²⁻⁴. As shown in Fig. 1, the most recent temperature-based estimates of ocean warming⁵⁻⁸ show good agreement for 2007-2016 (1.09 ± 0.10 to 1.16 ± 0.2 × 10²² J yr⁻<fnref n="1" />), but a larger spread when extending back to include the sparser 1990s data (0.90 ± 0.09 to 1.36 ± 0.10 × 10²² J yr⁻<fnref n="1" /> for 1993–2015). The spread is mostly caused by gap-filling methods and systematic errors⁵,⁸,⁹, which together introduce uncertainties up to 25–50% in warming trends¹⁰. Because temperature-based estimates use the same upper-ocean observations and linear warming trend for depths below 2000 m (ref. ⁶.), they may share additional unknown systematic errors⁸. An alternative method based on the top of the atmosphere energy balance¹¹ is also not truly independent, because it is subject to large systematic errors when estimating long-term trends and therefore depends on the same hydrographic measurements for calibration¹¹⁻¹⁴. Here we introduce a third method, based on changes in the abundances of gases in the atmosphere, which respond to whole-ocean warming through the temperature dependence of gas solubility in seawater. This method is not limited by data sparseness, because fast mixing in the atmosphere efficiently integrates the global ocean signal.

Changes in ocean heat content on seasonal¹⁵ and glacial-interglacial¹⁶ time-scales have been reconstructed using measurements of noble gases in modern or ancient air. Our method is similar, but instead of relying on noble gases (e.g. Ar/N₂), which lack sufficient accuracy as yet¹⁵, we rely on measurements of atmospheric O₂ and CO₂, which can be summed to yield a tracer “atmospheric potential oxygen” (APO) that responds to warming similar to a noble gas¹⁷. When the ocean warms, the solubility of O₂ and CO₂ drops, and the amount of gas lost by

<footnote n="1">¹Department of Geosciences and Princeton Environmental Institute, Princeton University, Princeton, USA. ²Scripps Institution of Oceanography, University of California San Diego, La Jolla, USA. ³Department of Environmental Science and Engineering, Fudan University, Shanghai, 200433, China. ⁴Laboratoire de Météorologie Dynamique/ Institut Pierre Simon Laplace, CNRS/ENS/X/UPMC, Département de Géosciences, Ecole Normale Supérieure, Paris, France. ⁵National Center for Atmospheric Research, Boulder, USA. ⁶NOAA, Geophysical Fluid Dynamics Laboratory, Princeton, USA. ⁷GEOMAR Helmholtz Centre for Ocean Research, Kiel, Germany. *email: laurer@princeton.edu</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Change in global ocean heat content (ΔOHC). (a) ΔOHC derived from hydrographic and atmospheric observations (normalized to zero in 2007, ±1-σ uncertainty). (b) Linear least-squares trends for 1991–2016, 1993–2016 and 2007–2016 (±1-σ uncertainty). Hydrography-based ΔOHC estimates combine warming rates at ocean depths of 0 to 2,000 m warming rates (Cheng and co-authors (CHEN)⁸, Pacific Marine Environmental Laboratory (PMEL)⁶, Meteorological Research Institute (MRI)⁵ and National Center for Environmental Information (NCEI)⁴⁷) with revised deep ocean (depth of more than 2000 m) of ref. ⁷ (Tables S1, S2). The atmospheric-based estimate which uses observed atmospheric potential oxygen trends (ΔAPOₐₖₐₘ) and model-based ΔAPOₐₖₐₘ-to-ΔOHC ratio, does not resolve interannual variations.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Processes contributing to observed changes in atmospheric potential oxygen (ΔAPOₐₗ). Industrial processes (fossil fuel burning and cement production, ΔAPOₐₗ) and ocean sink for anthropogenic carbon (ΔAPOₐₖₐₘ) remove APO from the atmosphere. The fertilization effect of anthropogenic aerosol deposition which promotes marine photosynthesis (ΔAPOₐₗ) and the changes in solubility, biology and ocean circulation due to warming (ΔAPOₐₖₐₘ) release APO into the atmosphere. This study shows that ΔAPOₐₖₐₘ can be used to estimate long term changes in global ocean warming.</figcaption></figure>

the ocean can be quantified with the complementary change observed in the atmosphere. Precise atmospheric O₂ measurements began in 1991 (CO₂ in 1958), enabling APO-based ocean heat content reconstructions that span nearly three decades¹⁸.

## Results

**Atmospheric potential oxygen trend and components.** APO (=O₂ + 1.05 × CO₂) is computed using observed atmospheric O₂/N₂ mole ratios and CO₂ mole fractions (see Methods)¹⁸⁻²⁰. By design, APO is insensitive to exchanges with land ecosystems, which produce changes in O₂ and CO₂ that largely cancel in APO owing to their approximate 1.05 O₂:C oxidative ratio. Time-series measurements at remote sites show a global long-term decline in APO ΔAPOₐₗ = 255.93 ± 13.74 per meg (units defined in Methods) between 1991 and 2016. ΔAPOₐₗ is driven by four primary contributors, illustrated in Fig. 2:

$$ \Delta APO_{OBS} = \Delta APO_{FF} + \Delta APO_{Cant} + \Delta APO_{AtmD} + \Delta APO_{Climate} \quad (1) $$
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Data-based estimates of global ΔAPO<sub>Climate</sub>. (a) ΔAPO<sub>Climate</sub> estimated from observed APO (ΔAPO<sub>OBS</sub>) from the Scripps Institution of Oceanography network (1991–2016), and corrected from fossil fuel burning, ocean anthropogenic carbon uptake and anthropogenic aerosols deposition (ΔAPO<sub>Climate</sub> = ΔAPO<sub>OBS</sub> –ΔAPO<sub>FF</sub> –ΔAPO<sub>Cant</sub> –ΔAPO<sub>AtmD</sub>) and their 1σ uncertainty ranges. (b) The increase in global ΔAPO<sub>Climate</sub> (±1σ interval) exceeds the range of 26-year trends expected from natural variations in four earth system models (CESM, GFDL, IPSL and UVic, shown in grey). Uncertainties and contributions to ΔAPO<sub>Climate</sub> are given in Tables S3, S4 and S6.</figcaption></figure>

where ΔAPO<sub>FF</sub> is the decrease in APO caused by industrial processes (fossil fuel burning and cement production) which in aggregate consume more than 1.05 moles of O<sub>2</sub> for each mole of CO<sub>2</sub> released; ΔAPO<sub>Cant</sub> accounts for the oceanic uptake of excess anthropogenic atmospheric CO<sub>2</sub>; APO<sub>AtmD</sub> accounts for for air-sea exchanges driven by ocean fertilization from anthropogenic aerosol deposition (increased fertilization leads to increased photo-synthesis, with a concomitant release of O<sub>2</sub> and uptake of CO<sub>2</sub>); and ΔAPO<sub>Climate</sub> accounts for air-sea fluxes of O<sub>2</sub>, CO<sub>2</sub> and N<sub>2</sub> driven by ocean processes, including warming-induced changes in solubility, in ocean circulation, and in photosynthesis and respiration (N<sub>2</sub> influences O<sub>2</sub>/N<sub>2</sub> ratios). Here, we derive ΔAPO<sub>Climate</sub> from Eq. (1) and show that it tracks ocean warming.

We estimate ΔAPO<sub>FF</sub> using fossil fuel and cement inventories<sup>21</sup>, finding ΔAPO<sub>FF</sub> = -142.38 ± 7.65 per meg (Fig. 3). ΔAPO<sub>Cant</sub> is controlled by the increase in atmospheric CO<sub>2</sub> and by ocean mixing, which is quantified by transient tracers distributions including chlorofluocarbons (CFCs)<sup>22</sup>; we find that ΔAPO<sub>Cant</sub> = -147.75 ± 3.69 per meg. ΔAPO<sub>Cant</sub> is relatively precise because it excludes the effect of changing ocean biology and circulation on natural carbon fluxes that are included in ΔAPO<sub>Climate</sub>. ΔAPO<sub>AtmD</sub> is derived from ocean model simulations with and without aerosols fertilization (phosphate, iron and nitrogen; Fig. S1)<sup>23</sup>. ΔAPO<sub>AtmD</sub> is uncertain, owing in part to uncertainties in iron availability to photosynthetic organisms, but is relatively small compared with the other terms: ΔAPO<sub>AtmD</sub> = 7.00 ± 3.50 per meg. From Eq. (1), we therefore find ΔAPO<sub>Climate</sub> = 27.21 ± 16.85 per meg, corresponding to a least-squares linear trend of +1.11 ± 0.68 per meg per year – larger than the trends expected from 26-year natural variations alone in four earth system models (the Community Earth System Model (CESM) and the Geophysical Fluid Dynamics Laboratory (GFDL), Institut Pierre Simon Laplace (IPSL) and University of Victoria (UVic) models). As shown in Fig. 3, a clear increase in ΔAPO<sub>Climate</sub> emerges over the period January 1991 to the end of December 2016.

**Linking atmospheric potential oxygen to ocean warming.** A starting point for understanding ΔAPO<sub>Climate</sub> is to imagine that O<sub>2</sub> and CO<sub>2</sub> behave like inert gases, such that the air-sea fluxes are dominated by temperature-driven solubility changes. In this case, APO would increase by around 0.8 per meg per 10<sup>22</sup> J of warming, with O<sub>2</sub> and CO<sub>2</sub> solubility changes accounting for an increase of +1.0 per meg per 10<sup>22</sup> J partly offset by the N<sub>2</sub> contribution of –0.2 per meg per 10<sup>22</sup> J (Methods). Support for the dominance of solubility in ΔAPO<sub>Climate</sub> can be found in the natural distribution of O<sub>2</sub> and carbon in the ocean. Ocean potential oxygen (OPO) is a dissolved tracer that mirrors APO<sub>Climate</sub> and tracks changes in air-sea O<sub>2</sub> and CO<sub>2</sub> fluxes<sup>17</sup>. Observed OPO abundance is strongly tied to ocean potential temperature (Fig. 4): warming induces OPO loss and cooling induces OPO gain. The OPO-to-temperature trend of –4.43 nmol J<sup>–1</sup> is within 18% of the trend of –3.65
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Observed link between potential oxygen and ocean heat. OPO concentrations *in-situ* (OPO, yellow) and at saturation based on O₂ and CO₂ solubility (OPOsat, grey) as a function of ocean temperature in the GLODAPv2 database⁴⁹.</figcaption></figure>

nmol J⁻¹ expected from solubility alone (OPOsat-to-temperature). Biological effects (related to changes in ocean circulation and photosynthesis/respiration) on CO₂ and O₂ substantially cancel in OPO (Fig. S2), while thermal impacts reinforce each other, with warming waters releasing both O₂ and CO₂ to the atmosphere and increasing ΔAPOₐtlate.

Further support for the dominance of solubility on ΔAPOₐtlate is found on multidecadal timescales in four Earth system models mentioned above, which yield OPO-to-temperature ratios between −4.69 and −4.36 nmol J⁻¹, bracketing the ratio of −4.43 nmol J⁻¹ found in hydrographic observations (Fig. S3). The models also simulate a very close relationship between ΔAPOₐtlate and the change in global ocean heat content (ΔOHC) that occurs during the simulations (1920–2100), with an atmospheric build-up in APO between 0.82 and 0.98 per meg per 10²² J (Figs. S3 and S4) – close to the ratio expected from temperature-driven solubility changes alone (0.8 per meg per 10²² J). By dividing the simulated APO change into separate biological and thermal components, we show that solubility changes account for more than 80% of ΔAPOₐtlate, while biologically driven changes account for 5% to 20% (Fig. S4). This partitioning found in response to transient warming is very similar to the partitioning found in hydrographic data (where solubility and biology contribute 82% and 18% respectively, to the OPO-to-temperature ratio, Fig. 4).

Small differences between individual model ΔAPOₐtlate-to-ΔOHC relationships (0.82 to 0.98 per meg per 10²² J) reflect systematic differences in biological fluxes. Models with stronger biological effects (IPSL and UVic) yield stronger oceanic loss of OPO and stronger release of APO for a given ocean warming (more negative OPO-to-temperature and higher ΔAPOₐtlate-to-ΔOHC, Fig. S3b). Using this relationship, we find that a ΔAPOₐtlate-to-ΔOHC ratio of 0.86 ± 0.03 per meg per 10²² J is compatible with observed OPO-to-temperature ratio. Combining this constrained ΔAPOₐtlate-to-ΔOHC ratio (0.86 ± 0.03 per meg per 10²² J) with the observation-based trend in ΔAPOₐtlate (1.11 ± 0.68 per meg/y), yields a warming trend of 1.29 ± 0.79 × 10²² Jyr⁻¹ between 1991 and 2016. As shown in Fig. 1, this APO-based estimate of ocean heat uptake is centered near the high end of recent ocean-based estimates, while agreeing with all estimates to within the uncertainties. In our estimate, the largest single source of uncertainty is the scale error from the span calibration of the O₂/N₂ analyzer (Table S6), which can be reduced via within-lab and inter-lab comparisons. The uncertainty in the estimated heat uptake since 1991 should therefore be significantly reduced in future updates.

Discussion

Most recent ocean-based estimates of warming suggest an increase relative to prior estimates¹,⁵,⁹. The independent APO estimate is fully in line with this upwards revision, although uncertainties are large. Our central value of ΔOHC would call for a steric sea level rise of 1.53 mm y⁻¹ (see Methods), also in agreement with satellite-derived constraints on thermal expansion, corrected for the fresh water contribution (1.50 ± 0.40 mm y⁻¹)²⁴,²⁵. Our result suggests that the ocean contributes 0.80 ± 0.49 W/m² to the Earth energy imbalance over the 1991–2016 period (Earth surface of 5.1 × 10¹⁴ m²), with implications for equilibrium climate sensitivity and the climate system lagged response to anthropogenic forcing²⁶,²⁷.

We find that the APO-heat coupling (APOₐtlate) is most robust on decadal and longer time-scales. Strong cancelation of biological O₂ and CO₂ fluxes is not expected on all temporal scales. On seasonal time-scales, air-sea O₂ fluxes driven by marine photosynthesis are ~8 times larger than those of CO₂ due to slow equilibration of CO₂²⁸. More complex coupling is also possible on interannual time-scales²⁹, such as the weaker lagged air-sea CO₂ flux compared to O₂ during El Nino events³⁰.
Atmospheric O₂ and CO₂ measurements have previously been applied to estimate global land and ocean CO₂ sinks, but relied on ocean heat content estimates and model-based oceanic O₂-to-heat ratios to correct for climate-driven O₂ outgassing³¹⁻³³. Here we have reversed this logic, using estimates of other quantities to constrain the ocean heating. Our approach exploits the APO-heat relationship, which is stronger than the O₂-heat relationship (See Methods for estimate of climate-driven ocean O₂ outgassing based on our results). Further work to constrain the separate contributions of O₂ and CO₂ to APO is needed to refine estimates of ocean and land carbon sinks using atmospheric O₂.

The results presented here were previously published in *Nature*, 563, 105–108 (2018)³⁴. The paper was subsequently retracted due to the underestimation of the uncertainties. Certain systematic errors were treated as random errors, and the uncertainty in the land oxidative ratio on the APO budget was neglected. These issues have been corrected in the current paper. The correction did not substantially change the central estimate of ocean warming but led to a roughly fourfold increase in uncertainties.

## Methods

**Observed changes in atmospheric potential oxygen (ΔAPO<sub>OBS</sub>).** A change in atmospheric potential oxygen concentration (in per meg) is defined following¹⁹:

$$ \Delta(\delta\text{APO})[\text{per meg}] = \Delta(\delta\text{O}_2/\text{N}_2) + \frac{\alpha_B}{\text{X}_{\text{O}_2}} \cdot (\text{X}_{\text{CO}_2} - 350) $$

with

$$ (\delta\text{O}_2/\text{N}_2)[\text{per meg}] = \frac{\text{O}_2/\text{N}_2(\text{sample})}{\text{O}_2/\text{N}_2(\text{reference})} - 1 $$

where Δ(δO₂/N₂) is the atmospheric change in δO₂/N₂ ratios (in per meg), XCO₂ is the CO₂ concentration in the air parcel (in ppm, i.e. µmol mol⁻¹) and 350 is an arbitrary reference, αB (=1.05) is the approximate O₂:CO₂ ratio of terrestrial ecosystems²⁰, and XO₂ (=0.2094) is the reference value of atmospheric mole fraction of O₂ necessary to convert XCO₂ from ppm to per meg units. We use an oxidative ratio αB of 1.05 rather than 1.1, following the recommendation of ref. ²⁰, who argue that woody biomass, which has an oxidative ratio near 1.05, likely dominates the long-term carbon sources and sinks on land.

ΔAPO<sub>OBS</sub> is computed from *in situ* atmospheric changes in CO₂ concentrations and O₂/N₂ ratios²⁰ measured at stations of the Scripps Institution of Oceanography network (available online at http://scrippso2.ucsd.edu)¹⁸,³⁵. The global average ΔAPO<sub>OBS</sub> is based on a weighted average of the three stations with the longest record (1991 to 2016), i.e. La Jolla (32.9°N, 117°W), Alert (82.5°N, 62.5°W) and Cape Grim (40.5°S, 144.5°E), with weights of 0.25, 0.25 and 0.5 respectively³⁶. Station's annual means are based on bi-monthly data fit to a four-harmonic seasonal cycle and a stiff long-term trend¹⁸. The uncertainty on ΔAPO<sub>OBS</sub> was computed by generating 10⁶ time-series with a noise scaled to the random and systematic errors of APO data detailed in Table S3. The uncertainty is taken as the 1-sigma interval (±1 standard deviation) from these 10⁶ realizations (Fig. 3). To assess whether the weighted average of the 3 stations accurately reflects the global trend, we also compared other combinations of stations from 1999 onward when we have reliable data from 9 stations. The differences were consistently smaller than 0.05 per meg yr⁻¹ which is negligible compared to the other combined measurement uncertainties.

**Fossil fuel burning and cement production influence on APO (ΔAPO<sub>FF</sub>).** ΔAPO<sub>FF</sub> is estimated using annual CO₂ emissions from oil, coal, gas, flaring and cement production (ΔCO₂(i) in moles)²¹ weighted by their O₂:C combustion ratios Rᵢ¹⁸:

$$ \Delta\text{APO}_{\text{FF}}[\text{per meg}] = \sum_{i} \frac{\alpha_{\text{B}} - R_{i}}{X_{\text{O2}}} \times \frac{\Delta\text{CO}_{2(i)}}{M_{\text{air}}} $$

where Mair is the number of moles of dry air in the atmosphere (convert moles of CO₂ to ppm).

The uncertainty on ΔAPO<sub>FF</sub> includes uncertainties in CO₂ emissions (ΔCO₂(i)³⁷ and in combustion ratios (Rᵢ in Tab S3)¹⁸. Uncertainties in ΔCO₂(i) and Rᵢ are considered not independent in time and were estimated using an autoregressive model with a correlation time-scale of 20 years following the method of ref. ³⁸ (1000 realizations using Monte-Carlo approach). The uncertainty on ΔAPO<sub>FF</sub> was then estimated by combining the 1000 realizations of ΔCO₂(i) and the 1000 realizations of Rᵢ, yielding a set of 10⁶ estimates of ΔAPO<sub>FF</sub>.

**Ocean anthropogenic carbon uptake influence on APO (ΔAPO<sub>Cant</sub>).** We can divide the ocean CO₂ uptake (ΔCO₂) as the sum of three contributions:

$$ \Delta\text{CO}_2 = \Delta\text{Cant}_0 + \Delta\text{Cant}' + \Delta\text{CO}_{2\text{Climate}} (\text{S1}) $$

where ΔCant₀ is the flux driven by the rise in CO₂ assuming steady ocean circulation (ΔCant₀ is negative, corresponding to uptake by the ocean), ΔCO₂Climate is the flux driven by the action of climate on natural carbon in the ocean (ΔCO₂Climate is positive, i.e. warming reduces the uptake of natural carbon), and ΔCant' is the remainder, which accounts for impact of circulation changes on the uptake of carbon driven by rising CO₂ (ΔCant' is positive, i.e. warming reduces the uptake of Cant). ΔAPO<sub>Cant</sub> can be expressed as the weighted sum of the two terms ΔCant₀ and ΔCant'.
$$ \Delta \text{APO}_{\text{Cant}} [\text{per meg}] = \frac{\alpha_{\text{B}}}{X_{\text{O}_2} \times M_{\text{air}}} \times (\Delta \text{Cant}_0 + \Delta \text{Cant}') $$

where $\Delta \text{Cant}_0$ and $\Delta \text{Cant}'$ are in moles. Note that $\Delta \text{CO}_{2\text{Climate}}$ is accounted for in $\Delta \text{APO}_{\text{Climate}}$.

$\Delta \text{Cant}_0$ is taken from the recent ocean inversion scheme with assimilation of observed potential temperature, salinity, radiocarbon, and CFC-11 of ref. 22 updated to 2016. $\Delta \text{Cant}'$ cannot be derived from observations and was estimated to be 0.05 PgC y$^{-1}$, equivalent to a trend of +0.11 per meg$^{-1}$, using model simulations (see details in section Model-based $\Delta \text{Cant}'$ below).

The uncertainty on $\Delta \text{APO}_{\text{Cant}}$ is related to uncertainties in $\Delta \text{Cant}_0$ and $\Delta \text{Cant}'$. We allow for uncertainty in $\Delta \text{Cant}_0$ following ref. 22 using the 10 sensitivity experiments (on ocean vertical and isopycnal diffusivities, data constraint, gas exchange coefficient etc.) available for the ocean inversion and an estimate of the interannual variability in the ocean sink of a 0.2 PgC y$^{-1}$. We also allow an additional 1% uncertainty (<0.03 PgC/y) in $\Delta \text{Cant}_0$ due to imperfectly known atmospheric $\text{CO}_2$ history$^{39}$, taking account of sensitivity to start date (1765 vs 1791), to degree of temporal smoothing, and to using different versions of the record since 1958 (Mauna Loa record versus average of Mauna Loa and South Pole records). This estimate used a variant of the box-diffusion model$^{40}$, and $\text{CO}_2$ data from ref. 41 and the Scripps $\text{CO}_2$ program (https://doi.org/10.6075/J0542KSG). Uncertainties on $\Delta \text{Cant}'$ are assumed to be 100% of the model-based estimate of $\Delta \text{Cant}'$ (see details in section Model-based $\Delta \text{Cant}'$ below).

### Ocean fertilization and atmospheric deposition of anthropogenic aerosol ($\Delta \text{APO}_{\text{AtmD}}$).

Deposition of anthropogenic aerosol from fossil fuel, biomass burning etc. fertilizes the ocean with nutrients, increases surface photosynthesis and sub-surface respiration$^{42-44}$. The effect of aerosol fertilization is partly counterbalanced by biological processes such as a decline in nitrogen fixation, which would be immediate, and an increase in denitrification in the water column, which would be on time-scales of several 100 years$^{45}$. Fixed anthropogenic nitrogen also fertilizes the land biosphere and coastal oceans by river runoffs, but in these cases, efficient denitrification returns fixed nitrogen to the atmosphere and has little impact on the APO budget on the decadal timescales considered here. The impact of anthropogenic aerosol on $\text{O}_2$, $\text{CO}_2$ and APO air-sea fluxes is evaluated with the IPSL ocean model NEMO-PISCES v2$^{46}$ using the difference between simulations with aerosols and a simulation in which the aerosol deposition is fixed to a constant pre-industrial value (equivalent to year 1850, Fig. S1)$^{23}$. We use four simulations with varying aerosols: one includes the combined effect of nitrogen (N), iron (Fe) and phosphorus (P) aerosol deposition, whereas the 3 others only include their individual contribution (N-only, Fe-only or P-only, Fig. S1 and Table S5). Uncertainties at 1-sigma level on $\Delta \text{APO}_{\text{AtmD}}$ are assumed to be ±50%. See Table S4.

Combined, N, Fe and P deposition accounts for an $\text{O}_2$ outgassing of 19.0 Tmol y$^{-1}$ for the 1980/2007 period (16 Tmol y$^{-1}$ for entire 1960/2007 simulation period) and an oceanic $\text{CO}_2$ uptake of 8.3 Tmoly$^{-1}$ for the 1980–2007 period (6.8 Tmol y$^{-1}$ for entire 1960–2007 simulation period, Fig. S1 and Table S5). The overall impact is $\Delta \text{APO}_{\text{AtmD}} = +0.28$ per meg y$^{-1}$ over 27 years of simulation (1980 to 2007), which we extrapolate to our 1991 to 2016 period. Increased $\text{O}_2$ outgassing accounts for an increase in APO of +0.51 per meg y$^{-1}$, and $\text{CO}_2$ uptake accounts for a change in APO of -0.23 per meg y$^{-1}$ ($\text{APO}_{\text{AtmD}(\text{O}_2)}$ and $\text{APO}_{\text{AtmD}(\text{CO}_2)}$ in Table S3).

The overall effect of N, Fe and P is smaller than the sum of the individual effects (Fig. S1), because of the interplay between the aerosol deposition pattern and nutrient co-limitations in the ocean. Phytoplankton growth in the ocean depends on the availability of the most limiting nutrient. While more available N will promote photosynthesis in regions where N is limiting (for example the tropical Atlantic Ocean), the effect is negligible in regions where Fe, P or any other nutrient are limiting (for example the Southern Ocean) (see Fig. 2 in ref. 23).

To our knowledge this is the first estimate of the impact of anthropogenic aerosol deposition on both $\text{O}_2$ and $\text{CO}_2$ air-sea fluxes at the global scale. Note however that ref. 18 used anthropogenic aerosol N inventories and scaling arguments to estimate an ocean $\text{O}_2$ loss due to anthropogenic N-deposition only of ~10 ± 10 Tmol y$^{-1}$, slightly weaker than our model estimate of 15.5 Tmol y$^{-1}$.

### $\Delta \text{APO}_{\text{Climate}}$ trends and uncertainty analysis. We compute the APO response to climate change ($\Delta \text{APO}_{\text{Climate}}$) via

$$ \Delta \text{APO}_{\text{Climate}} = \Delta \text{APO}_{\text{OBS}} - \Delta \text{APO}_{\text{FF}} - \Delta \text{APO}_{\text{Cant}} - \Delta \text{APO}_{\text{AtmD}} $$

We combine the estimates of $\Delta \text{APO}_{\text{FF}}$, $\Delta \text{APO}_{\text{Cant}}$ and $\Delta \text{APO}_{\text{AtmD}}$ plus estimates of the contribution from variations in the oxidative ratio to obtain 10$^6$ time-series of $\Delta \text{APO}_{\text{FF}} + \Delta \text{APO}_{\text{Cant}} + \Delta \text{APO}_{\text{AtmD}}$ and obtain 10$^6$ time-series of $\Delta \text{APO}_{\text{Climate}}$ using the 10$^6$ time-series of $\Delta \text{APO}_{\text{OBS}}$. We computed the $\Delta \text{APO}_{\text{Climate}}$ trend and its uncertainty based on the distribution of the unweighted least square fits to each of the 10$^6$ ensemble realizations of $\Delta \text{APO}_{\text{Climate}}$ generated by combining all sources of uncertainty. We find a $\Delta \text{APO}_{\text{Climate}}$ trend of 1.11 ± 0.68 per meg yr$^{-1}$ for 1991–2016. The individual contributions to the uncertainty are shown in Table S6.

### Hydrography-based estimates of ocean heat uptake ($\Delta \text{OHC}$). We used four global-ocean estimates of $\Delta \text{OHC}$ based on hydrographic measurements in Fig. 1. Surface to 2000 m warming rates are from ref. 6 (PMEL), ref. 5 (MRI, climate.mri-jma.go.jp/pub/ocean/ts/v7.2/), an updated version of ref. 47 (NCEI, www.nodc.noaa.gov/OC5/3M_HEAT Content/basin_avt_data.html), and ref. 8 (CHEN, http://159.226.119.60/cheng/images_files/TOA_OHC_errorbar_1940_2015_2.txt), with the revised deep ocean (depth > 2000 m) constant linear warming rate of 0.10 ± 0.03 × 10$^{22}$ Jyr$^{-1}$ of ref. 7 based on the global ship-based sections program (GO-SHIP; http://www.go-ship.org)$^{48}$.
**Ocean observations of ocean potential oxygen (OPO).** We used *in-situ* ocean observations from GLODAPv2⁴⁹ combined with an anthropogenic carbon estimate²² interpolated at the location of each sample to compute 78,456 values (Glodapv2 quality control = 0, marginal seas and coastal waters were removed) of Oceanic Potential Oxygen (OPO)¹⁷ via

$$OPO = O_2 + \alpha_B \times C_{pi}^*$$

where $O_2^*$ and $C_{pi}^*$ are the ocean conservative tracers related to air-sea fluxes of $O_2$ and pre-industrial carbon⁵⁰, and $\alpha_B$ is the terrestrial oxidative ratio ($\alpha_B = 1.05$). The thermal component (solubility-driven) of OPO (OPO$_{sat}$) is computed as:

$$OPO_{sat} = O_{2sat} + \alpha_B \times C_{pisat}$$

where $O_{2sat}$ is the dissolved $O_2$ concentration at saturation with the observed temperature and salinity⁵¹ and $C_{pisat}$ is the dissolved inorganic carbon concentration expected at the observed temperature and salinity, and assuming equilibrium with a pre-industrial partial pressure of $CO_2$ of 280 ppm and using pre-formed alkalinity⁵².

**Solubility-driven changes in oceanic and atmospheric potential oxygen.** Figure S2 shows a tight and quasi-linear link between observed OPO and potential temperature (−4.43 nmol/J, $r^2 = 0.95$), similar to the link found between $OPO_{sat}$ and potential temperature (−3.65 nmol/J, $r^2 = 0.93$). This suggests that changes in OPO and hence $\Delta APO_{Climate}$ are driven primarily by changes in thermal air-sea fluxes. In these observations, departures of dissolved oxygen and carbon concentrations ($O_2^*$ and $C_{pi}^*$) from their respective saturation curves ($O_{2sat}$ and $C_{pisat}$) due to biological activity tend to balance (Fig. S2). By contrast, thermal effects reinforce each other ($O_{2sat}$ and $C_{pisat}$ both decrease with increasing temperature) and biological effects compensate each other ($O_2^* > O_{2sat}$ and $C_{pi}^* < C_{pisat}$).

Changes in APO expected from changes in gas solubility in the ocean is an increase of 2.95 nmol per J of warming, which includes the outgassing of $O_2$ and $CO_2$ following $OPO_{sat}$ (3.65 nmol/J) and the release of $N_2$ (0.6 nmol/J) (Fig. S2b). A change of 2.95 nmol per J of warming is equivalent to an increase of 0.8 per meg/10²² J of APO in the atmosphere ($=(2.95 \times 10^{-9})/(3.7 \times 10^{19}) \times 10^{22} = 1.0 \times 10^{-6} = 0.8$ per meg per 10²² J, with $3.7 \times 10^{19}$ the number of moles of $O_2$ in the atmosphere). $O_2$ and $CO_2$ solubility alone yield an increase in APO of 1.2 per meg/10²² J, which is partly counterbalanced by the outgassing of $N_2$ that decreases APO by 0.2 per meg/10²² J (via the increase of the $O_2/N_2$ ratio).

**Earth system model experiments.** We used 4 Earth-system models (ESMs): the Geophysical Fluid Dynamics Laboratory Earth System Models with a nominally level vertical coordinate version GFDL-ESM2M (called GFDL here)⁵³,⁵⁴, the Institut Pierre-Simon Laplace Coupled Model 5 version IPSL-CM5A-LR (IPSL here)⁵⁵, the Community Earth System Model large ensemble CESM-LE (CESM here)⁵⁶ and the UVic model version 2.9 (UVic here)⁵⁷. Evaluation of these models and their biogeochemical components can be found in previous studies⁵⁴,⁵⁶,⁵⁸⁻⁶⁰. GFDL, IPSL and UVic participated in the Coupled Model Intercomparison Project Phase 5 (CMIP5)⁶¹.

For GFDL, IPSL and UVic, we used the CMIP5 business as usual “historical-RCP8.5” scenario, the feedback experiment “esmFdbk3” that only includes warming-driven changes associated with anthropogenic emissions (e.g. radiation effects) and the fixed-climate experiment “esmFixClim3” that only includes the direct biogeochemical effects of increasing atmospheric $CO_2$ (e.g. uptake of anthropogenic carbon, acidification etc.). For CESM, we also used the historical and RCP8.5 experiments and the separation between anthropogenic carbon from the natural carbon available in this model (carbon tracer separation approach). The feedback approach used for GFDL, IPSL and UVic removes all direct biogeochemical effects of rising atmospheric $CO_2$ on the air-sea $O_2$ and $CO_2$ exchanges, whereas the natural carbon tracer separation approach used for CESM still includes the biogeochemical impacts of increasing atmospheric $CO_2$ on the carbon cycle (e.g. acidification) even while it excludes the anthropogenic carbon itself. However, we expect this effect to be small and negligible on our results.

We also used the multicentury preindustrial control simulation “piControl” with no increase in atmospheric $CO_2$ to correct for model drift and to estimate the natural internal variability of $\Delta APO_{Climate}$ (Fig. 2). We used model results over the 1920–2100 period, which were available for the four models.

Model OPO was computed as in the observations. Note that for CESM we removed subsurface regions of high denitrification in the eastern equatorial Pacific and Bay of Bengal where oxygen and $O_2^*$ in this model have unrealistic values⁶².

**Model anthropogenic $\Delta Cant'$.** The component $\Delta Cant'$ was derived from Eq. (S1) ($\Delta Cant' = \Delta CO_2 - \Delta Cant_0 - \Delta CO_{2Climate}$) using CMIP5 model simulations. $\Delta CO_2$ was taken from experiment RCP8.5, $\Delta Cant_0$ from experiment esmFixClim3, and $\Delta CO_{2Climate}$ from experiment esmFdbk3. Note that the control simulation was also used to correct model drift. We estimated $\Delta Cant' = 0.05 \pm 0.05$ PgC yr⁻¹ for the 1991 to 2016, based on the results of the three models, which individually yields $\Delta Cant'$ of 0.0 PgC yr⁻¹ (IPSL), 0.11 PgC yr⁻¹ (GFDL) and 0.11 PgC yr⁻¹ (UVic), and assuming an uncertainty of ±100%. This corresponds to a trend of 0.11 ± 0.11 per meg yr⁻¹.

**Model $\Delta APO_{Climate}$ to $\Delta OHC$ ratios and uncertainty.** Model $\Delta APO_{Climate}$ is computed using individual contributions from $O_2$, $CO_2$ and $N_2$ via:
$$ \Delta APO_{\text{Climate}}[\text{per meg}] = \Delta APO_{(O_2)} + \Delta APO_{(CO_2)} + \Delta APO_{(N_2)} $$

$$ \Delta APO_{\text{Climate}}[\text{per meg}] = \frac{1}{X_{O_2} M_{\text{air}}} \times \left( \Delta F_{O_2} + \alpha_B \times \Delta F_{CO_2} - \frac{X_{O_2}}{X_{N_2}} \times \Delta F_{N_2} \right) $$

where $\Delta F_{O_2}$, $\Delta F_{CO_2}$ and $\Delta F_{N_2}$ are the changes in air-sea fluxes of O₂, CO₂ and N₂ respectively (in moles), $M_{\text{air}}$ is the number of moles of dry air in the atmosphere and $X_{N_2}$ and $X_{O_2}$ are the reference atmospheric mixing ratio of N₂ and O₂ respectively⁶³. O₂ and CO₂ fluxes are simulated in the models. N₂ air-sea fluxes, which impact the O₂ atmospheric mixing ratio (because O₂ is ~20% of the atmospheric composition), are quantified from the global ocean temporal changes in N₂ solubility computed from model changes in temperature and salinity⁶⁴.

The link between long-term changes in APO<sub>Climate</sub> and ocean heat content, i.e. ΔAPO<sub>Climate</sub>-to-ΔOHC ratios, were computed for each model using the 180 years of simulations (1920 to 2100). Resulting ΔAPO<sub>Climate</sub>-to-ΔOHC ratios vary between 0.82 and 0.98 per meg per 10²² J of warming (Fig. S3). These ratios include uncertainty in the natural climate variations at interannual and decadal time-scales and uncertainty in the O₂:C oxidative ratio associated with global gains and losses of O₂ and CO₂ by terrestrial ecosystems ($\alpha_B = 1.05 \pm 0.05$). The uncertainty due to interannual variations was evaluated by computing ΔAPO<sub>Climate</sub>-to-ΔOHC ratios using multiple 26 year-long segments from the 180 year simulations. We obtained 616 ΔAPO<sub>Climate</sub>-to-ΔOHC ratios (154 time-series of 26 years per model) and used the standard deviation between these ratios as a measure of the uncertainty.

We combine results from all models along with modeled and observed hydrographic relationships between OPO and potential temperature to establish an observationally-constrained estimate of the ΔAPO<sub>Climate</sub>-to-ΔOHC ratio, as shown in Fig. S3b. We first carry out this procedure with a fixed oxidative ratio $\alpha_B$ of 1.05, which yields ΔAPO<sub>Climate</sub>-to-ΔOHC = 0.86 ± 0.03 per per 10²² J. Repeating the same procedure using oxidative ratios $\alpha_B$ of 1.0 and 1.1 applied to both APO and OPO yield differences in ΔAPO<sub>Climate</sub>-to-ΔOHC smaller than 0.01 per meg per 10²² J, showing that sensitivity to the oxidative ratio is very small. Note that on Fig. S3b, the model error bars include two contributions to the uncertainties on the simulated ΔAPO<sub>Climate</sub>-to-ΔOHC ratios (interannual variations and O₂:C ratio) that combine to yield ±0.01 per meg per 10²² J for the CESM and GFDL models, ±0.02 per meg per 10²² J for the UVic model and ±0.05 per meg per 10²² J for the IPSL model (1-sigma).

**Steric component of sea level rise.** We evaluated the steric component of sea level rise associated with an OHC of $1.3 \times 10^{22}$ Jyr⁻¹ to be 1.53 mm y⁻¹. Following ref. ⁶⁵, this calculation assumes that 45% of the warming occurs below 700 m and that the steric rise is 1 mm per $0.60 \times 10^{22}$ J above 700 m and 1 mm per $1.15 \times 10^{22}$ J below 700 m (i.e. global steric rise of 1 mm per $0.84 \times 10^{22}$ J). Assuming that 48% of the warming occurs below 700 m (ref. ⁶⁵) would yield a global steric rise of 1 mm per $0.86 \times 10^{22}$ J and change our estimate by less than 3%. Our estimate is also consistent with the recent hydrography-based estimate of the WCRP Global Sea Level Budget Group⁶⁶.

**Link to global ocean de-oxygenation.** Our application of O₂ atmospheric measurements to constrain long-term ocean warming can be compared with earlier work considering warming-driven oceanic O₂ outgassing. Multiplying our warming rate of $1.29 \pm 0.79 \times 10^{22}$ Jyr⁻¹ by the O₂-to-heat ratios simulated by the four ESMs ($-3.70 \pm 0.80$ nmol O₂ J⁻¹), yields an ocean loss of $48 \pm 30$ Tmol O₂ y⁻¹. Adding a loss of ~19 Tmol O₂ y⁻¹ due to anthropogenic aerosols (Table S5), yields a global ocean outgassing of $67 \pm 35$ Tmol O₂ y⁻¹, in the range of previous estimates based on atmospheric data⁶⁷ (~40 Tmol O₂ y⁻¹), ocean data above 1000 m⁶⁸,⁶⁹ (55–65 Tmol O₂ y⁻¹) and global ocean data⁷⁰ (96 ± 42 Tmol O₂ y⁻¹). This calculation implies that ocean CO₂ uptake is reduced by warming at a ratio of ~0.70 nmol of CO₂ per Joule (difference between O₂-to-heat ratio of 3.70 nmol J⁻¹ and OPO-to-heat ratio of 4.43 nmol J⁻¹).

## Data availability

Scripps APO data are available at scrippso2.ucsd.edu/apo-data. Model results are available upon reasonable request to R. W. (IPSL anthropogenic aerosol simulations), L. B. (IPSL-CM5A-LR), M. C. L. (CESM-LE), J. P. D. (GFDL-ESM2M) and W. K. (UVic).

## Code availability

The code used to compute APO<sub>Climate</sub> (https://doi.org/10.5281/zenodo.2571986) is available at https://zenodo.org/record/2571986#.XQ1FYZNKhBw. Earth-system model codes are available online for IPSL-CM5A-LR (cmc.ipsl.fr/ipsl-climate-models), GFDL-ESM2M (mdl-mom5.herokuapp.com/web/docs/project/quickstart), UVic(climate.uvic.ca/model) and CESM (cesm.ucar.edu/models/).

Received: 29 August 2019; Accepted: 12 December 2019;

Published online: 27 December 2019

## References

1. Climate change 2013: the physical science basis: Working Group I contribution to the Fifth assessment report of the Intergovernmental Panel on Climate Change. (Cambridge University Press, 2014).

2. Abraham, J. P. et al. A review of global ocean temperature observations: Implications for ocean heat content estimates and climate change. Rev. Geophys. 51, 450–483 (2013).

3. Riser, S. C. et al. Fifteen years of ocean observations with the global Argo array. Nat. Clim. Change 6, 145–153 (2016).

4. Levitus, S., Antonov, J. I., Boyer, T. P. & Stephens, C. Warming of the World. Ocean. Science 287, 2225–2229 (2000).

5. Ishii, M. et al. Accuracy of Global Upper Ocean Heat Content Estimation Expected from Present Observational Data Sets. SOLA 13, 163–167 (2017).
6. Johnson, G. C. et al. Ocean Heat Content [in "State of the Climate in 2016"). Am. Meteorol. Soc. Bull. **98**, S66–S68 (2017).

7. Desbruyères, D. G., Purkey, S. G., McDonagh, E. L., Johnson, G. C. & King, B. A. Deep and abyssal ocean warming from 35 years of repeat hydrography. *Geophys. Res. Lett.* **43**, 2016GL070413 (2016).

8. Cheng, L. et al. Improved estimates of ocean heat content from 1960 to 2015. *Sci. Adv.* **3**, e1601545 (2017).

9. Cheng, L. et al. XBT Science: Assessment of Instrumental Biases and Errors. *Bull. Am. Meteorol. Soc.* **97**, 924–933 (2016).

10. Boyer, T. et al. Sensitivity of Global Upper-Ocean Heat Content Estimates to Mapping Methods, XBT Bias Corrections, and Baseline Climatologies. *J. Clim.* **29**, 4817–4842 (2016).

11. Allan, R. P. et al. Changes in global net radiative imbalance 1985–2012. *Geophys. Res. Lett.* **41**, 5588–5597 (2014).

12. Palmer, M. D. Reconciling Estimates of Ocean Heating and Earth's Radiation Budget. *Curr. Clim. Change Rep.* **3**, 78–86 (2017).

13. Loeb, N. G. et al. Observed changes in top-of-the-atmosphere radiation and upper-ocean heating consistent within uncertainty. *Nat. Geosci.* **5**, 110 (2012).

14. Trenberth, K. E., Fasullo, J. T. & Balmaseda, M. A. Earth's Energy Imbalance. *J. Clim.* **27**, 3129–3144 (2014).

15. Battle, M. et al. Measurements and models of the atmospheric Ar/N₂ ratio. *Geophys. Res. Lett.* **30**, 1786 (2003).

16. Ritz, S. P., Stocker, T. F. & Severinghaus, J. P. Noble gases as proxies of mean ocean temperature: sensitivity studies using a climate model of reduced complexity. *Quat. Sci. Rev.* **30**, 3728–3741 (2011).

17. Resplandy, L. et al. Constraints on oceanic meridional heat transport from combined measurements of oxygen and carbon. *Clim. Dyn.* **47**, 3335 (2016).

18. Keeling, R. F. & Manning, A. C. Studies of Recent Changes in Atmospheric O₂ Content. In Treatise on Geochemistry 385–404 (Elsevier, 2014).

19. Stephens, B. B. et al. Testing global ocean carbon cycle models using measurements of atmospheric O₂ and CO₂ concentration. *Glob. Biogeochem. Cycles* **12**, 213–230 (1998).

20. Randerson, J. T. et al. Is carbon within the global terrestrial biosphere becoming more oxidized? Implications for trends in atmospheric O₂. *Glob. Change Biol.* **12**, 260–271 (2006).

21. Quéré, C. L. et al. Global Carbon Budget 2016. Earth Syst. Sci. Data **8**, 605–649 (2016).

22. DeVries, T. The oceanic anthropogenic CO₂ sink: Storage, air-sea fluxes, and transports over the industrial era. Glob. Biogeochem. Cycles **28**, 631–647 (2014).

23. Wang, R. et al. Influence of anthropogenic aerosol deposition on the relationship between oceanic productivity and warming. *Geophys. Res. Lett.* **42**, 10745–10754 (2015).

24. Church, J. A. et al. Sea Level Change. In Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (eds. Stocker, T. F. et al.) 1137–1216, https://doi.org/10.1017/CBO9781107415324.026 (Cambridge University Press, 2013).

25. Rietbroek, R., Brunnabend, S.-E., Kusche, J., Schröter, J. & Dahle, C. Revisiting the contemporary sea-level budget on global and regional scales. *Proc. Natl. Acad. Sci.* **113**, 1504–1509 (2016).

26. Hansen, J. et al. Earth's Energy Imbalance: Confirmation and Implications. *Science* **308**, 1431–1435 (2005).

27. Forster, P. M. Inference of Climate Sensitivity from Analysis of Earth's Energy Budget. *Annu. Rev. Earth Planet. Sci.* **44**, 85–106 (2016).

28. Keeling, R. F. & Severinghaus, J. P. Atmospheric oxygen measurements and the carbon cycle. In The carbon cycle (Global Change Institute, Proceedings on the Carbon Cycle) 134–140 (Cambridge University Press, New York, 2000).

29. Resplandy, L., Séférian, R. & Bopp, L. Natural variability of CO₂ and O₂ fluxes: What can we learn from centuries-long climate models simulations? *J. Geophys. Res. Oceans* **120**, 384–404 (2015).

30. Eddebbar, Y. A. et al. Impacts of ENSO on air-sea oxygen exchange: Observations and mechanisms. *Glob. Biogeochem. Cycles* **31**, 2017GB005630 (2017).

31. Keeling, R. F. & Garcia, H. E. The change in oceanic O₂ inventory associated with recent global warming. *Proc. Natl. Acad. Sci.* **99**, 7848–7853 (2002).

32. Bopp, L., Le Quéré, C., Heimann, M., Manning, A. C. & Monfray, P. Climate-induced oceanic oxygen fluxes: Implications for the contemporary carbon budget. *Glob. Biogeochem. Cycles* **16**, 6–1–6–13 (2002).

33. Keeling, C. D., Piper, S. C., Whorf, T. P. & Keeling, R. F. Evolution of natural and anthropogenic fluxes of atmospheric CO₂ from 1957 to 2003. *Tellus B* **63**, 1–22 (2011).

34. Resplandy, L. et al. Quantification of ocean heat uptake from changes in atmospheric O₂ and CO₂ composition. *Nature* **563**, 105–108 (2018).

35. Keeling, R. F., Manning, A. C., Paplawsky, W. J. & Cox, A. C. On the long-term stability of reference gases for atmospheric O₂/N₂ and CO₂ measurements. *Tellus B* **59**, 3–14 (2007).

36. Hamme, R. C. & Keeling, R. F. Ocean ventilation as a driver of interannual variability in atmospheric potential oxygen. *Tellus B* **60**, 706–717 (2008).

37. Andres, R. J., Boden, T. A. & Higdon, D. A new evaluation of the uncertainty associated with CDIAC estimates of fossil fuel carbon dioxide emission. *Tellus B* **66** (2014).

38. Ballantyne, A. P. et al. Audit of the global carbon budget: estimate errors and their impact on uptake uncertainty. *Biogeosciences* **12**, 2565–2584 (2015).

39. Bronselaer, B., Winton, M., Russell, J., Sabine, C. L. & Khatiwala, S. Agreement of CMIP5 Simulated and Observed Ocean Anthropogenic CO₂ Uptake. *Geophys. Res. Lett.* **44**, 12,298–12,305 (2017).

40. Oeschger, H., Siegenthaler, U., Schotterer, U. & Gugelmann, A. A box diffusion model to study the carbon dioxide exchange in nature. *Tellus* **27**, 168–192 (1975).

41. MacFarling Meure, C. et al. Law Dome CO₂, CH₄ and N₂ O ice core records extended to 2000 years BP. *Geophys. Res. Lett.* **33** (2006).

42. Wang, D., Gouhier, T. C., Menge, B. A. & Ganguly, A. R. Intensification and spatial homogenization of coastal upwelling under climate change. *Nature* **518**, 390–394 (2015).

43. Ito, T., Nenes, A., Johnson, M. S., Meskhidze, N. & Deutsch, C. Acceleration of oxygen decline in the tropical Pacific over the past decades by aerosol pollutants. *Nat. Geosci.* **9**, 443–447 (2016).

44. Jickells, T. D. et al. A reevaluation of the magnitude and impacts of anthropogenic atmospheric nitrogen inputs on the ocean. *Glob. Biogeochem. Cycles* **31**, 2016GB005586 (2017).

45. Somes, C. J., Landolfi, A., Koeve, W. & Oschlies, A. Limited impact of atmospheric nitrogen deposition on marine productivity due to biogeochemical feedbacks in a global ocean model. *Geophys. Res. Lett.* **43**, 4500–4509 (2016).

46. Aumont, O., Ethé, C., Tagliabue, A., Bopp, L. & Gehlen, M. PISCES-v2: an ocean biogeochemical model for carbon and ecosystem studies. *Geosci. Model Dev.* **8**, 2465–2513 (2015).

47. Levitus, S. et al. World ocean heat content and thermosteric sea level change (0–2000 m), 1955–2010. *Geophys. Res. Lett.* **39**, L10603 (2012).

48. Talley, L. D. et al. Changes in Ocean Heat, Carbon Content, and Ventilation: A Review of the First Decade of GO-SHIP Global Repeat Hydrography. *Annu. Rev. Mar. Sci.* **8**, null (2016).

49. Olsen, A. et al. The Global Ocean Data Analysis Project version 2 (GLODAPv2) – an internally consistent data product for the world ocean. *Earth Syst. Sci. Data* **8**, 297–323 (2016).

50. Sarmiento, J. L. & Gruber, N. Sinks for anthropogenic carbon. *Phys. Today* **55**, 30–36 (2002).

51. Garcia, H. E. & Gordon, L. I. Oxygen solubility in seawater: Better fitting equations. *Limnol. Oceanogr.* **37**, 1307–1312 (1992).
52. Gruber, N., Sarmiento, J. L. & Stocker, T. F. An improved method for detecting anthropogenic CO2 in the oceans. *Glob. Biogeochem. Cycles* **10**, 809–837 (1996).

53. Dunne, J. P. et al. GFDL's ESM2 Global Coupled Climate-Carbon Earth System Models. Part I: Physical Formulation and Baseline Simulation Characteristics. *J. Clim.* **25**, 6646–6665 (2012).

54. Dunne, J. P. et al. GFDL's ESM2 Global Coupled Climate-Carbon Earth System Models. Part II: Carbon System Formulation and Baseline Simulation Characteristics. *J. Clim.* **26**, 2247–2267 (2013).

55. Séférian, R., Iudicone, D., Bopp, L., Roy, T. & Madec, G. Water Mass Analysis of Effect of Climate Change on Air-Sea CO₂ Fluxes: The Southern Ocean. *J. Clim.* **25**, 3894–3908 (2012).

56. Kay, J. E. et al. The Community Earth System Model (CESM) Large Ensemble Project: A Community Resource for Studying Climate Change in the Presence of Internal Climate Variability. *Bull. Am. Meteorol. Soc.* **96**, 1333–1349 (2014).

57. Keller, D. P., Oschlies, A. & Eby, M. A new marine ecosystem model for the University of Victoria Earth System Climate Model. *Geosci Model Dev* **5**, 1195–1220 (2012).

58. Bopp, L. et al. Multiple stressors of ocean ecosystems in the 21st century: projections with CMIP5 models. *Biogeosciences* **10**, 6225–6245 (2013).

59. Keller, D. P., Kriest, I., Koeve, W. & Oschlies, A. Southern Ocean biological impacts on global ocean oxygen. *Geophys. Res. Lett.* **43**, 6469–6477 (2016).

60. Long, M. C., Deutsch, C. & Ito, T. Finding forced trends in oceanic oxygen. *Glob. Biogeochem. Cycles* **30**, 381–397 (2016).

61. Taylor, K. E., Stouffer, R. J. & Meehl, G. A. An Overview of CMIP5 and the Experiment Design. *Bull. Am. Meteorol. Soc.* **93**, 485–498 (2011).

62. Moore, J. K., Lindsay, K., Doney, S. C., Long, M. C. & Misumi, K. Marine Ecosystem Dynamics and Biogeochemical Cycling in the Community Earth System Model [CESM1(BGC)]: Comparison of the 1990s with the 2090s under the RCP4.5 and RCP8.5 Scenarios. *J. Clim.* **26**, 9291–9312 (2013).

63. Rödenbeck, C., Le Quéré, C., Heimann, M. & Keeling, R. F. Interannual variability in oceanic biogeochemical processes inferred by inversion of atmospheric O₂/N₂ and CO₂ data. *Tellus B* **60**, 685–705 (2008).

64. Hamme, R. C. Mechanisms controlling the global oceanic distribution of the inert gases argon, nitrogen and neon. *Geophys. Res. Lett.* **29** (2002).

65. Trenberth, K. E., Fasullo, J. T., von Schuckmann, K. & Cheng, L. Insights into Earth's Energy Imbalance from Multiple Sources. *J. Clim.* **29**, 7495–7505 (2016).

66. WCRP Global Sea Level Budget Group. Global Sea Level Budget 1993 – Present. *Earth Syst. Sci. Data Discuss.* 1–88, https://doi.org/10.5194/essd-2018-53 (2018).

67. Keeling, R. F., Körtzinger, A. & Gruber, N. Ocean Deoxygenation in a Warming World. *Annu. Rev. Mar. Sci.* **2**, 199–229 (2010).

68. Helm, K. P., Bindoff, N. L. & Church, J. A. Observed decreases in oxygen content of the global ocean. *Geophys. Res. Lett.* **38**, L23602 (2011).

69. Ito, T., Minobe, S., Long, M. C. & Deutsch, C. Upper ocean O₂ trends: 1958–2015. *Geophys. Res. Lett.* **44**, 4214–4223 (2017).

70. Schmidtko, S., Stramma, L. & Visbeck, M. Decline in global oceanic oxygen content during the past five decades. *Nature* **542**, 335–339 (2017).

## Acknowledgements

We thank Mike Winton and eight anonymous reviewers for their useful suggestions. L.R. gratefully acknowledges support of the Climate Program Office of the National Oceanic and Atmospheric Administration grant NA13OAR4310219 and the Princeton Environmental Institute. NCAR is sponsored by the National Science Foundation. We also thank the people maintaining the APO time-series at Scripps andoping th the groups devele models CESM, GFDL, IPSL and UViC used in this study.

## Author contributions

L.R. directed the analysis of the several data sets and models used here and shared responsibility for writing the manuscript; R.F.K. shared responsibility for writing the manuscript; R.W. performed the simulations of anthropogenic aerosols; L.B., J.P.D., M.C.L., W.K. and A.O. provided model results. All authors contributed to the final version of the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s41598-019-56490-z.

Correspondence and requests for materials should be addressed to L.R.

Reprints and permissions information is available at www.nature.com/reprints.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2019