Articles | Volume 21, issue 3
Research article
12 Feb 2024
Research article |  | 12 Feb 2024

Underestimation of multi-decadal global O2 loss due to an optimal interpolation method

Takamitsu Ito, Hernan E. Garcia, Zhankun Wang, Shoshiro Minobe, Matthew C. Long, Just Cebrian, James Reagan, Tim Boyer, Christopher Paver, Courtney Bouchard, Yohei Takano, Seth Bushinsky, Ahron Cervania, and Curtis A. Deutsch

The global ocean's oxygen content has declined significantly over the past several decades and is expected to continue decreasing under global warming, with far-reaching impacts on marine ecosystems and biogeochemical cycling. Determining the oxygen trend, its spatial pattern, and uncertainties from observations is fundamental to our understanding of the changing ocean environment. This study uses a suite of CMIP6 Earth system models to evaluate the biases and uncertainties in oxygen distribution and trends due to sampling sparseness. Model outputs are sub-sampled according to the spatial and temporal distribution of the historical shipboard measurements, and the data gaps are filled by a simple optimal interpolation method using Gaussian covariance with a constant e-folding length scale. Sub-sampled results are compared to full model output, revealing the biases in global and basin-wise oxygen content trends. The simple optimal interpolation underestimates the modeled global deoxygenation trends, capturing approximately two-thirds of the full model trends. The North Atlantic and subpolar North Pacific are relatively well sampled, and the simple optimal interpolation is capable of reconstructing more than 80 % of the oxygen trend in the non-eddying CMIP models. In contrast, pronounced biases are found in the equatorial oceans and the Southern Ocean, where the sampling density is relatively low. The application of the simple optimal interpolation method to the historical dataset estimated the global oxygen loss to be 1.5 % over the past 50 years. However, the ratio of the global oxygen trend between the sub-sampled and full model output has increased the estimated loss rate in the range of 1.7 % to 3.1 % over the past 50 years, which partially overlaps with previous studies. The approach taken in this study can provide a framework for the intercomparison of different statistical gap-filling methods to estimate oxygen content trends and their uncertainties due to sampling sparseness.

1 Introduction

Historical observations indicate that the ocean oxygen (O2) inventory has declined in recent decades, a trend that has been termed ocean deoxygenation (Keeling et al., 2010; Levin, 2018). Ocean heat uptake causes the reduction of oxygen solubility and changes in ocean circulation and biogeochemical processes. Ocean warming and increasing stratification can further decrease O2 exchange between upper and deep layers, further reducing the oceanic O2 inventory. The reduction in dissolved oxygen can have far-reaching impacts on marine habitats (Deutsch et al., 2015; Gruber, 2011; Pörtner and Farrell, 2008; Vaquer-Sunyer and Duarte, 2008).

Figure 1Maps of (a) single-year observation for the O2 anomaly in the year 1991 at 200 m depth; (b) the pentadal composite O2 anomaly centered in the year 1991, covering the period from 1989 to 1993, at 200 m depth; and (c) the optimally interpolated pentadal O2 anomaly based on the data in the middle figure.

The distribution of historical O2 measurements is irregular and sparse. The calculation of changes in the global O2 content requires filling the data gaps in time and space, making it difficult to quantify global trends and their uncertainties. Recent estimates of the global oxygen decline are in the range of 0.5 %–3.3 % (IPCC, 2022) relative to climatological means over the period of 1970–2010 (Helm et al., 2011; Schmidtko et al., 2017; Ito et al., 2017). The wide range in the estimates of ocean deoxygenation can stem from different interpolation methods to estimate global O2 content, different data quality control standards, and different data sources. Previous studies estimating the rates of ocean deoxygenation have relied on the World Ocean Database 2018 (WOD18) (Boyer et al., 2018). WOD represents an international collaboration among national data centers, oceanographic research institutions, and investigators to provide a comprehensive dataset of quality-controlled oceanographic variables. Shipboard observations are more prevalent in the Northern Hemisphere oceans in the warm seasons. Oxygen measurements from a single year (e.g., 1991; Fig. 1a) do not adequately cover the global ocean; a pentadal composite (e.g., 1989–1993; Fig. 1b) is performed to increase the coverage at the expense of averaging out the high-frequency variability on timescales shorter than 5 years. Even so, there are large data gaps in the South Pacific and Indian Ocean. In such a case, optimal interpolation (hereafter, OI) has been widely applied to fill data gaps and to yield a gridded data field (Fig. 1c), which produces the best-fit O2 distribution in the least-square sense, given the covariance structure in the dataset (Wunsch, 1996). One shortcoming of OI application is that it can underestimate O2 trends in data-sparse regions. For regions without any nearby measurements, the mapped field approaches the climatology asymptotically (i.e., towards a zero-oxygen anomaly). If there is a widespread O2 decline but only a fraction of ocean volume is sampled, the OI method will underestimate the declining trend of ocean O2 content.

The objective of this study is to use a suite of Earth system model (ESM) simulations as a test bed to evaluate the uncertainties in ocean deoxygenation rates by sub-sampling model output according to the spatial and temporal distributions of the historical shipboard measurements. Earth system models represent our current understanding of physical and biogeochemical processes, expressed in mathematical equations. These processes and their interactions are numerically integrated forward in time, predicting the trajectory of the Earth's climate system. ESMs generate their own natural variability that reflects chaotic behavior of the natural climate system, but its temporal trajectory does not necessarily match that of the real world. Observed O2 changes may be influenced by both external forcing (such as volcanism and anthropogenic greenhouse gases and aerosol emissions) and natural climate variability. These models are imperfect and often include varying degrees of biases due to inadequate process understanding and the lack of computational resources to resolve critical processes at smaller length and/or timescales. Current Earth system models do not fully reproduce the O2 variability and trends (Oschlies et al., 2018, 2017), and observational data are essential for the evaluation of the model output. In turn, the analysis of model output can inform the range of underlying variability and trends.

This study uses seven different ESMs from the Coupled Model Intercomparisperon Project Phase 6 (CMIP6) that provided dissolved oxygen output. These seven models sample the range of O2 variability and trends that can arise from different model architectures, biogeochemical parameterizations, and modes and phases of natural climate variability. Globally gridded O2 fields from ESMs provide fully sampled states and, thus, perfectly known model trends for the simulated variables. The modeled O2 distribution can also be sub-sampled according to the time-evolving pattern of historical ocean observations to evaluate the effect of sampling sparseness. We purposefully remove information from the model output where there was no in situ measurement. This hypothetical observation of model output, with its realistic data gaps, can be used to evaluate the uncertainties in ocean deoxygenation rates due to both data sparseness and statistical gap-filling approaches. The sub-sampled model output can then be subject to a statistical gap-filling method (OI) to evaluate how well the fully sampled states can be reconstructed. It is of great interest to evaluate to what extent the OI method underestimates the true O2 trend in the context of the simulated deoxygenation.

The structure of this paper is as follows. The second section describes the analysis method, the data sources, and the Earth system models. The third section describes the results, followed by the interpretation of the results and a conclusion in section four.

Figure 2Number of OSD and CTD oxygen profiles aggregated into 1× 1 longitude–latitude grid cells for 4 decades, from 1965 to 2014. The color scale indicates the number of measurements in log scale.

2 Methods

2.1 Observational data source

We make use of observations from the bottle and conductivity–temperature–depth instruments' (CTD) O2 data in WOD18. Dissolved oxygen is the third most frequently measured chemical tracer in the ocean, following temperature and salinity. There are approximately 2.8 million temperature, 2.4 million salinity, and 0.9 million O2 vertical profiles in the ocean station data (OSD or simply bottle data) reported in WOD18. In addition, CTD data include approximately 1 million temperature and salinity profiles, and 0.2 million O2 profiles. The OSD (i.e., bottle) O2 data are largely located on the margins of the ocean basins and along repeat hydrographic transects (Fig. 2). O2 observations in the OSD profile were typically measured by means of a modified Winkler titration method with a precision of about 1 µmol kg−1 (Carpenter, 1965). Most modern oxygen chemical titration measurements are based on Carpenter's whole-bottle titration method and an amperometric or photometric end detection with a precision of about 0.5–1 µmol kg−1 (or approximately, 0.3 %). The CTD O2 data are based on electrochemical and optical sensors mounted on the CTD rosette samplers, which are periodically calibrated to the Winkler O2 (Grégoire et al., 2021). The coverage of CTD measurements increased after the 1990s, and that of profiling floats has rapidly increased in recent years. However, the overall spatio-temporal coverage of O2 observations from bottle and CTD has decreased since the 1990s. Profiling-float O2 data have increased significantly in the past 10 years; however, the precision is on the order of 1 %–2 % ( 2 µmol kg−1), and the data quality control and calibration are still under development, especially in the upper-ocean oxycline (Bittig et al., 2018; Maurer et al., 2021). Float O2 data have been excluded in this study but will be an important data source, especially after 2010s for the future studies.

2.2 Data pre-processing and optimal interpolation

The pre-processing of the data includes a check for acceptable data quality using the WOD18 quality control (QC) flags. The original WOD18 standard-depth profiles with 102 depth levels are placed into bins which are the 1× 1 longitude–latitude grid cells with 102 vertical depth levels referenced to the standard depths of WOD18. Of the 102 vertical depth levels, 47 levels are in the upper 1000 m.

The target analysis period is after 1965 when the modern oxygen titration method was established by Carpenter, as referenced above. Some of the data from most recent years are not included in the ESMs, as discussed below, so the analysis ends in 2014. The spatially binned quality-controlled data were averaged at a monthly resolution where mean, variance, and sample size are recorded from 1965 to 2014 for the bottle data and from 1987 to 2014 for the CTD O2 data. Next, the monthly mean climatology is determined by calculating the climatological monthly mean combining the bottle and CTD O2 data and then filling data gaps. We are interested in long-term O2 changes which can be calculated as the anomalies from the monthly climatological mean. Departures from the monthly climatology are recorded as O2 anomalies for each bin. The binned data are very sparse at a monthly timescale (Fig. 1). For each year, the monthly anomaly data are averaged into yearly anomalies, neglecting the months with missing data. This step increases spatial data coverage significantly while averaging out high-frequency variability in the data, including changes shorter than the yearly timescale, such as waves and eddies. In addition, a 5-year moving window (pentadal) averaging is applied to the yearly anomaly, neglecting the years with missing data. This further increases the spatial data coverage while averaging out variability on the timescale shorter than 5 years. The resulting, pentadal O2 anomaly data cover the 46-year period from 1967 to 2012.

A relatively simple optimal interpolation (OI) is applied to the pentadal O2 anomaly data for each year to yield the spatially interpolated O2 anomalies following Wunsch (1996). This method provides the least-square estimate of the O2 field on regularly spaced grid cells, minimizing the mean square error of the mapped data for given observations with a covariance function. Stationary and isotropic Gaussian covariance is assumed throughout this study, with the e-folding length scale (Lref) of 1000 km. This particular choice of length scale controls how far an observation can influence the far field together with the assumed noise-to-signal ratio (ε) of 0.2. The Gaussian assumption may be qualitatively reasonable, but the ocean circulation is neither spatially stationary nor uniform. The use of the Gaussian function allows us to avoid calculating and storing the large and complex covariance structure, but it can distort the resulting maps (Fukumori et al., 1991), which is a caveat for this study. A basin mask is used to interpolate data points only within the same ocean basin such as the Atlantic, Pacific, Indian, and Southern oceans. Each 1× 1 grid point is assigned to one of the 53 basins defined in Appendix 1 of Garcia et al. (2019). The binned oxygen vector, X, is expressed as a (N× 1) vector, where N is the number of binned data for a particular basin. The objective map of oxygen climatology, Y, is a (M× 1) vector, where M is the number of grid cells for the basin. The optimal interpolation is applied to each basin as follows:

(1) Y = D E - 1 X ,

where X is the pentadal oxygen anomaly input from the discrete data, and Y is the objective map of (gap-filled) oxygen anomalies. D is a (M×N) data–grid covariance matrix based on the Gaussian function, where Dmn=exp(-Lmn2/Lref2 and Lmn is the distance between the two points. C follows the same definition but for the N×N data–data covariance, and E=C+εI, where ε is the noise-to-signal covariance ratio. For the Southern Ocean, all data points southward of 30 S are used. An example of this process in the year 1991 is shown in Fig. 1. Basin-wise application of optimal interpolation is performed for the O2 anomalies, resulting in yearly (running pentadal) maps for the 46-year period. The O2 anomaly field and its standard error field are recorded.

2.3 Ocean deoxygenation trend

Using the yearly maps of the O2 anomaly field, global and basin-wise O2 contents are calculated as the volume integral over the upper 1000 m, O(t), where t is time since 1967. The magnitude is referenced to the mean value of the first 10 years, where the 10-year (1967–1976) mean O2 contents are subtracted from respective O2 content time series for comparison purposes. Ocean deoxygenation trends are estimated as the slope (a) of the O2 content time series using standard linear regression:


where σtO is the covariance between time and O2, and σtt is the variance in time. (a,b) are the slope and intercept of linear regression. Assuming that the regression errors are normally distributed, the standard error for the slope (ϵa) and intercept (ϵb) can be calculated as follows:


where MSE stands for the mean square error of regression, tn is time at nth data point. The gridded O2 dataset is constructed based on a 5-year running mean. An effective sample size (Neff) is calculated assuming that 5-year data are independent, thus Neff  9 for 46 years of data. These parameters are later used to evaluate the uncertainty and will be used for the comparison between models and observation.

2.4 CMIP6 Earth system models

Two sets of time-varying O2 fields are derived from the ESMs, including the full field and the reconstructions from the sub-sampled model output (Table 1). We selected a subset of Earth system models participating in the Coupled Model Intercomparison Project Phase 6 (CMIP6), and the outputs for their historical simulation are downloaded from the Earth System Grid Federation (, last access: March 2023). The monthly mean O2 output is first re-gridded onto the global 1× 1 longitude–latitude grid for the period of 1965 to 2014. A bilinear interpolation is first performed for the horizontal interpolation, followed by the linear interpolation on the vertical axis to the standard depths of the WOD18. Sub-sampled model output is then generated from the full field where the model output fields were resampled using the same spatial and temporal locations as with the observations. The sub-sampling strategy assumes that a grid box is sampled if, at least, one observation exists within the grid cell at a particular year or month. If so, we retain model data in the sub-sampled dataset. In reality, there could be multiple casts within the same grid and the same year or month, but multiple samples and/or variability within a single cell are not considered. There are slight differences in the land–ocean masks between models, and we use the model topography as they are provided. Similarly to the observational analysis, the sub-sampled monthly O2 climatology is assembled from the sub-sampled data with the optimal interpolation filling the data gaps using Eq. (1). Then O2 anomalies are calculated by subtracting the sub-sampled monthly climatology, and they are first aggregated into annual O2 anomalies, neglecting months without data, followed by the running pentadal averaging. Finally, the basin-wise optimal interpolation is applied to yield the reconstructed O2 anomaly fields using Eq. (2). These procedures are repeated for each of the models in Table 1. For the comparison purposes, the 5-year moving window averaging is applied to the full field.

Table 1List of CMIP6 models used in this study. Variants represent decadal-scale variability ensemble members and are coded according to (r) realization, (i) initialization, (p) physics, and (f) forcing. The first available variant, typically noted as r1i1p1f1, is taken from each model.

Download Print Version | Download XLSX

Figure 3(a) Linear trend of upper-ocean (0–1000 m) column O2 inventory from 1967 to 2012 from the optimally interpolated pentadal O2 anomaly based on the World Ocean Database 2018. (b) Time series of oxygen inventory is plotted for the global domain from 0 to 1000 m depth, including its linear trend and the 95 % confidence interval of the trend line. The confidence interval is calculated using a Monte Carlo method.

3 Results

3.1 Observed and modeled O2 trend maps

The observational trend is first determined based on the optimally interpolated gap-filled WOD18 profiles. The vertically integrated O2 inventory (0–1000 m) trend pattern is shown in Fig. 3a. While regional differences exist, the basin-scale patterns of the observed O2 loss are similar to those in previous studies (Helm et al., 2011; Schmidtko et al., 2017; Ito et al., 2017). In the North Atlantic, overall O2 decline is observed, except for the south of Greenland in the subpolar North Atlantic, where a patch of an increasing trend exists. In the North Pacific, a strong decrease is found in the western subpolar region, spreading from the sea of Okhotsk (Nakanowatari et al., 2007), which may be connected to the reduced ventilation in this region. A weak increase is found in the subtropical North Pacific (Ito et al., 2019), which is related to the multi-decadal natural variability of the North Pacific climate. Oxygen increases are observed in the subtropical Southern Hemisphere oceans and to the south of Greenland. In terms of the global inventory trends, the data suggest a global linear trend of 175 ± 24 Tmol O2 decade−1 or approximately 1.5 % loss over the 50-year period.

Table 2From left to right column: spatial pattern correlation (Pearson's correlation coefficients) between observed and modeled upper-ocean (0–1000 m) column O2 trend and the global trend magnitudes, the pattern correlation between observed modeled O2 trend patterns (full and sub-sampled optimally interpolated model outputs).

Download Print Version | Download XLSX

Figures 5 and 6 show the comparison of the trend pattern between the models listed in Table 1 for the full model field (Fig. 5) and the reconstructed model output (Fig. 6). The modeled O2 trend patterns are moderately correlated to the observations for some of the CMIP6 models, as summarized in Table 2. CanESM5, MPI-ESM1-2-LR, IPSL-CM6A, and MIROC-ES2L exhibit a moderately positive correlation of approximately r=0.3. It is interesting to contrast this result to the hindcast simulation of the earlier generation of ocean biogeochemistry models reported by Stramma et al. (2012). The subset of CMIP6 models in this study is slightly better correlated to observational estimates than the hindcast runs using an earlier generation of models. This is likely due to the improved biogeochemical model structure and parameterization rather than the physical climate forcing. Hindcast simulations are forced by the observed atmospheric variability through the meteorological reanalysis products. In contrast, historical simulations of the CMIP6 models generate natural climate variability that, in general, does not reproduce the phasing of observed variability.

The reconstructed CMIP6 model output is slightly better correlated to the observation than the full model output for the majority (five out of seven) of models, perhaps reflecting the common sampling pattern and gap-filling approach. Comparing the reconstructed and the full field from the same model, the pattern correlation of the O2 trend ranges from 0.37 to 0.69. While it is not perfect, the OI can estimate the general pattern of the full field with moderate correlation for the 1× 1 gridded trend maps. This motivates us to further investigate to what extent OI can estimate the O2 trend for a larger-scale hemispheric and/or global domain.

Figure 4Time series of upper-ocean (0–1000 m) column O2 inventory from CMIP6 models from (solid line) full model output and (dash line) sub-sampled and optimally interpolated model output. Panel (a) is the global, and panels (b) and (c) are the Northern and Southern hemispheres. The range of the vertical axis for the hemispheric inventories is smaller than the global inventory. The inventory anomaly is referenced in relation to the first 10-year averages.


3.2 Global and hemispheric O2 inventory time series

The globally integrated O2 content has a stronger declining trend than in ESMs, and the weak trend bias in models becomes even greater when reconstructed from sub-sampled data (Fig. 4, Table 2). Only one of the models exceeded the observed global trend in the full field (CanESM5, 165 Tmol decade−1; Table 2). When there are no observations nearby, the OI reverts to the background climatology, thus decreasing the amplitude of anomalies. Thus, the estimated O2 content tends to underestimate O2 anomalies in the region of sparse sampling. The magnitude of underestimation depends on the distance from observations, which sets the covariance according to the assumed Gaussian function. Figure 4 further shows that the sub-sampling introduces three decadal-scale peaks in the years 1988, 2000, and 2011 for both observations and some of the models (Fig. 4). These quasi-decadal peaks are not apparent in the full model output. We hypothesize that these quasi-decadal peaks are likely spurious, caused by the sparse sampling pattern. The magnitudes of these apparent spurious peaks are on the order of 100 Tmol O2. To provide a context, they are comparable to the anomalous O2 inventory increase caused by the eruption of Mt. Pinatubo and subsequent ocean cooling and enhanced ocean O2 uptake (Fay et al., 2023).

The global inventory time series is divided into the northern and southern hemispheric components (Fig. 4b and c). Comparing the hemispheric and global inventory time series indicates some notable issues with the sub-sampling. First, northern hemispheric trends in some of the full model outputs (IPSL-CM6A-LR and CanESM5) have similar magnitudes to the observational trends. In some other ESMs, the overall magnitudes of the southern hemispheric trends are similar to the observations (GFDL-ESM4 and CanESM5). Overall, the hemispheric trends demonstrate similar magnitudes between the north and the south for observations and models. The reconstructed model outputs appear to underestimate the magnitude of the trends for all models. The magnitude and the causes of this underestimation are of great interest and will be investigated further in the following sections.

Secondly, the observed quasi-decadal peaks primarily appear in the southern hemispheric inventory (Fig. 4e and f, solid black line), and some of the models reproduce these peaks (GFDL-ESM4, MIROC-ES2L, NorESM-LM2, E3SM1-1) for the reconstructed model output (Fig. 4f). There are no apparent peaks in the full model output, confirming that these features are spurious.

Thirdly, in the Northern Hemisphere (Fig. 4c and d), there is a moderate increase towards the late 1980s, and then it decreases strongly during the 1990s. Two of the Earth system models (GFDL-ESM4 and E3SM1-1) show similar increases in the early period (Fig. 4c); however, they underestimate the decreasing trend after the 1990s. These features are distorted in the reconstructed model output. It is difficult to determine whether the apparent increase in oxygen content is meaningful during the 1980s, but similar features are found in earlier studies focusing on the near-surface waters (Garcia et al., 2005).

Figure 5Modeled linear trend of the column O2 inventory (0–1000 m) from 1967 to 2014 from the full model output.

The models and observations tend to disagree more significantly in the Southern Hemisphere. Modeled inventory trends disagree substantially from one another because of the spurious quasi-decadal noises. While some models exceed the observed magnitude of oxygen decline (GFDL-ESM4, CanESM5), some other models even show increases in the southern hemispheric O2 inventory (NorESM2-LM, IPSL-CM6A-LR).

Figure 6Same as Fig. 5 but for the trend reconstructed from optimal interpolation of sub-sampled model outputs.

3.3 Spatial pattern of O2 trends and basin-wise inventories

To examine O2 trends across ocean basins, we divided the global data into 13 regions according to a basin mask (shown in Fig. S1 in the Supplement). The basin O2 inventories are integrated for each region from the full model output and sub-sampled model (Fig. S2 and S3). Figure 5 shows the spatial patterns of the column O2 trend (0–1000 m) from the full model output. Blue shading shows strong O2 loss, and red indicates O2 increase. The patterns of O2 trends from the reconstructed model output are displayed in Fig. 6. The effect of sub-sampling does not change the spatial pattern; it only affects the trend magnitudes. As expected, the reconstructed model outputs exhibit weaker trend magnitudes. For each region, the inventory time series are displayed in separate figures from Figs. S4 through S16. For the basin-scale deoxygenation trend, the North Atlantic Ocean is the only basin where all models show the same sign of change relative to the observation for the full-field and reconstructed model outputs (Figs. S2 and S3).

Figure 7Spatial data coverage of the pentadal O2 anomaly data from WOD18. The area covered by grid cells with at least one profile is divided by the total area for each basin. Blue lines are Atlantic basins. Magenta lines indicate Pacific basins. Indian basins are in yellow, and cyan is used for the Southern Ocean and the Arctic Ocean. Basin masks are defined in Fig. S1 and are coded by color. Abbreviations for the basin names are as follows: subpolar North Atlantic (SPNA), subtropical North Atlantic (STNA), equatorial Atlantic (EQAT), subtropical South Atlantic (STSA), Mediterranean Sea (MED), subpolar North Pacific (SPNP), subtropical North Pacific (STNP), equatorial Pacific (EQPA), subtropical South Pacific (STSP), equatorial Indian Ocean (EQID), subtropical South Indian Ocean (STSI), Southern Ocean (SO), and Arctic Ocean (AO).


Figure 7 shows the evolution of the spatial data coverage for each basin. To calculate the percent coverage value, the area of grid cells with at least one shipboard profile is divided by the total area of grid cells in each basin. Overall, the North Atlantic and Mediterranean Ocean are the most well observed among the 13 regions. Near-surface waters are better sampled than the deeper layers (400, 700 m). The data coverage evolves over time, depending on the basin. During the 1970s and 80s, there was greater data coverage for near-surface waters (100–200 m), and the near-surface data coverage gradually decreased after the 1980s. However, this pattern is not uniform through the depths. For some regions such as the subpolar North Atlantic, there appears to be no significant decrease in the deeper profile (700 m).

There are several notable features from this comparison. First, models exhibit varying patterns of O2 changes, and the model disagreements are more pronounced in the Southern Hemisphere oceans, even in the full model output. This is consistent with the varying hemispheric-scale trend magnitude, as shown in Fig. 4e and f. Approximately half of the models show increasing or decreasing trends in the subtropical South Pacific. The observed O2 decline is strong in the observations, but its time series is noisy, and the discrepancies between the full-field and the reconstructed model output are large in the Southern Ocean (Fig. S17). This is consistent with the persistently low data coverage in the Southern Ocean (Fig. 7). The Southern Ocean contributes significantly to the spurious quasi-decadal peaks that are visible in the hemispheric and global time series (Fig. 4); thus, the observed trend in the Southern Ocean may include large uncertainty.

Figure 8Basin-wise relationship between fully sampled and sub-sampled O2 trends for seven CMIP6 models. Data points on or near the 1:1 line (solid red) indicate that sub-sampled data adequately reproduced the fully sampled modeled trend. Abbreviations for the basin names are as follows: subpolar North Atlantic (SPNA), subtropical North Atlantic (STNA), equatorial Atlantic (EQAT), subtropical South Atlantic (STSA), Mediterranean Sea (MED), subpolar North Pacific (SPNP), subtropical North Pacific (STNP), equatorial Pacific (EQPA), subtropical South Pacific (STSP), equatorial Indian Ocean (EQID), subtropical South Indian Ocean (STSI), Southern Ocean (SO), and Arctic Ocean (AO).


There are two regions, namely the subpolar and subtropical North Atlantic, where observations and all models agree in terms of the sign of changes. These two regions' inventory time series are displayed in Figs. S4 and S5. In the subpolar North Atlantic, the magnitude of modeled O2 changes brackets the observation, whereas some models (CanESM5, IPSL-CM6A-LR, E3SM1-1) exhibit even stronger O2 loss than observations. In the subtropical North Atlantic, these three models exhibit a similar magnitude of O2 loss as the observations. In the equatorial Atlantic, there is a clear difference between the models and observation. The observation shows a decreasing trend, and not a single model was able to reproduce it. Similarly, the models were not able to reproduce the magnitude of O2 loss in the subpolar North Pacific, with the exception of MPI-ESM-1-2-LR.

3.4 Synthesis

The basin-wise O2 trend is compared between the full-field and reconstructed model output in Fig. 8, assessing the ability of the OI method to reproduce the full-field data. In Fig. 8, the horizontal axis is the full model, and the vertical axis is the reconstructed model output. Each dot indicates simulated O2 trend magnitude for a basin. The solid red line is the 1:1 ratio, indicating where the OI method was able to fully reproduce the trend magnitude. Most of the dots are located between the solid red line and the dashed purple line, indicating that the magnitude of the ocean deoxygenation trend is underestimated due to the OI method applied to sparsely sampled data.

Four regions (subtropical North Atlantic, subpolar North Atlantic, Mediterranean, subpolar North Pacific) performed very well in terms of capturing more than 80 % of the deoxygenation trend in the context of the simulation. These regions are relatively well sampled, and the loss of the trend magnitude due to the OI is minimal. In contrast, the four regions (equatorial Atlantic, equatorial Pacific, equatorial Indian, Southern Ocean) performed very poorly, capturing less than 30 % of the simulated deoxygenation trend. These regions unfortunately are not well represented by the sub-sampled and gap-filled data, showing the limitation of the OI method. The strong negative trend in the Southern Ocean (upper-left panel in Fig. 5) may be highly uncertain, and this is concerning since the Southern Ocean significantly contributes to the global oxygen content. Other basins (subtropical South Atlantic, subtropical North Pacific, subtropical South Pacific, subtropical Indian Ocean, Arctic Ocean) are moderately represented (30 %–80 % of the true trend).

Figure 9Global relationship between fully sampled and sub-sampled model O2 trend. Blue dots indicate the seven CMIP6 models with full global model output. Purple dots indicate the same, except the four poorly represented regions (equatorial Atlantic, equatorial Pacific, equatorial Indian, Southern Ocean) are excluded.


To what extent did the OI method underestimate the global deoxygenation trend? Figure 9 illustrates the relationship between the true global trend (on the x axis) and the estimated trend from the sub-sampled model outputs. Each dot comes from a model from two different ways of aggregating the global trend. The blue dots include all basins regardless of the ability of the OI method to reconstruct the true trend. The purple dots exclude the equatorial basins, as well as the Southern Ocean. The linear regression among the seven models informs us that the sub-sampling and the gap filling with the OI method can capture approximately two-thirds (68 %, purple line) of the true trend, excluding the low-confidence regions. Looking at the distribution among the models, the spread of this ratio is 19 %, as calculated by the standard deviation. If all basins are included, the fraction that is retrieved by the OI method decreases to 58 % (blue line).

One of the implications from the model analysis is that the optimal interpolation method used in this study may result in the significant underestimation of the dissolved oxygen trend in observations. The observation-based global oxygen content trend can be adjusted assuming that the ratio of the deoxygenation trend between the sub-sampled and full model output is approximately two-thirds (68 ± 19 %), as determined by the CMIP6 ESMs. Optimal interpolation of the WOD18 oxygen profiles estimated a 1.5 % O2 decline over the last 50 years, but the true O2 decline may be in the range of 1.7 % to 3.1 %. This partially overlaps with the recent estimates of the global oxygen decline, which is in the range of 0.5 %–3.3 % (IPCC, 2022), but suggests that the low end of that range is very unlikely.

4 Discussion

The premise of this study is that Earth system models can provide useful information about the uncertainties in the global ocean deoxygenation rate due to the sparse sampling and the specific gap-filling method used with observations. The models disagree amongst each other and with the observations because of different and imperfect representations of processes due to model structures, parameterizations, and the presence of natural variability. However, a model can estimate the observational sampling bias by comparing its true model state to one that is reconstructed by sub-sampling the model output according to the pattern of shipboard profiles (bottle and CTD) from the WOD18.

The model-based analyses in this study are generally consistent in showing that sub-sampling with the gap-filling method yields weaker trends than the full model output. The gap-filling method used in this study is a relatively simple implementation of optimal interpolation (OI), which provides the best-fit distribution of O2 anomaly in the least-square sense assuming a Gaussian covariance structure of the data. Our mapping approach is admittedly simple, but this choice has certain benefits, for example, that the results from a simple method are easy to understand and that it is also easy to notice and to correct mistakes. It can be replicated by other groups relatively easily. If the ocean deoxygenation has a wide-spread, large-scale signal, as well as regional hotspots, we anticipate that a simple method should, at least, capture the majority of the large-scale components and some regional features. There are some drawbacks in that it tends to smooth out spatial gradients, and it may not represent regional signals very well in data-poor regions. This OI method essentially predicts a diminishing anomaly when there is no observation nearby with the assumed e-folding length scale. If there is a widespread O2 decrease, the OI can underestimate the trend in a sparsely sampled region. Our result confirmed this tendency for the global deoxygenation trends from the subset of CMIP6 Earth system models. Our analysis, based on seven such models, suggests that approximately two-thirds of the true trend are captured by the reconstructed model output. This conclusion generally applies to all models independently of the model skills in capturing the observed trend for the global and hemispheric inventories (compare left and right column of Fig. 4) and for the basin-wise trends (compare Figs. 5 and 6). Some ocean regions have better coverage than others, and significant regional variations exist for the sampling density and thus the performance of OI. For example, the North Atlantic and subpolar North Pacific are relatively well sampled, and the OI was able to capture more than 80 % of the true trend. In these well-sampled regions, detailed analyses of ocean deoxygenation rates are likely to be fruitful using models and observations using the OI method.

Broadly speaking, the Northern Hemisphere oceans are generally better sampled than the Southern Hemisphere oceans, but the overall trends appear to be equally contributed by both hemispheres (Fig. 4). Basin-wise analysis revealed diverging basin-wise trend patterns among the models (Figs. 5 and 6). There is no consistent pattern in the contributions from different regions to the overall trend regardless of the sampling density. Also, there is no consistent pattern in the sign of multi-decadal O2 trends, except for the North Atlantic Ocean, where all models are in general agreement. This region has the highest sampling density (Fig. 7), and the full-field and reconstructed O2 trends are in good agreement (Fig. 8). Data coverage is not the only factor, but it plays an important role in the performance of the OI method. The North Atlantic is sampled at 20 %–50 % density based on 1× 1 grid cells with decreasing coverage from the surface to deeper depths and from subpolar to subtropical latitudes (Figs. 7 and S17). In the low-sample region, namely, the Southern Ocean, whose data coverage is persistently less than 13 %, the OI method struggled to reconstruct the full-field O2 trends. In this region, historical observations are limited to certain longitudes or latitudes (e.g., Drake Passage) and the repeat hydrographic cruises (Fig. 2), and it was clearly inadequate to represent the full-field data.

It is useful to compare oxygen content trends using multiple gap-filling approaches to assess the uncertainties from different methodologies (IPCC, 2022). The framework developed in this study may be helpful to further deepen such intercomparison studies and to quantify the skill of different gap-filling methods in the context of model output. Such comparison studies may reveal what sampling density is sufficient to reconstruct the real trend. For such an exercise, it is important to select model-derived oxygen fields that include realistic background variability. For the OI method, it is also crucial to have the covariance structure of the O2 field. An important caveat for this study is that the models used here were not eddy resolving, and we also used a Gaussian covariance with a prescribed length scale. Mesoscale ocean eddies are energetic features with characteristic spatial scales of 10–100 km and characteristic timescales of several months. The model outputs did not include this type of internal ocean variability, and the modeled fields did not include the mesoscale noises that are present in the observation. The use of non-eddying models reduces the level of internal ocean variability to be much lower than the observations. Thus, we are not able to address to what extent the trend estimates vary depending on the presence of ocean eddies and smaller-scale variability, which is a caveat in this study. It may be possible to emulate the mesoscale eddy noise (and uncertainties from other factors such as instrumental errors) and to estimate a more realistic covariance structure of O2 fields using outputs from detrended high-resolution simulations, but this is beyond the scope of this paper and is left for future study.

This paper has focused on mapping and the trends of total O2, which consist of two sub-components, oxygen saturation (O2sat) and apparent oxygen utilization (AOU). O2sat strongly depends on temperature, with a minor contribution from salinity, and explains less than half of the observed O2 trend (Ito et al., 2017; Schmidtko et al., 2017). Most uncertainty is associated with the AOU component since temperature is measured with much higher sampling rates, and its mapping uncertainty would be significantly lower than that of O2 or AOU. This paper has also focused on the mapping in depth coordinates. While it is beyond the scope of this paper, interpolating O2 along isopycnal surfaces may potentially reduce the mapping uncertainty. Temperature variation on an isopycnal is much smaller than that of O2 so the O2sat variation would be better constrained along isopycnals. Also, ocean transport in the interior ocean is primarily oriented along isopycnals, and the interpolation on isopycnal surfaces can potentially reduce the spurious errors. However, there could be technical difficulties as the bottle O2 measurements come from discrete bottle samples, and the sampling depths are unlikely to match the location of desired isopycnals, leading to an interpolation error. In the end, one would have to try and evaluate how much uncertainty can be reduced by mapping along density horizons, which would be a promising topic for future study.

While sampling sparseness is likely a major source of uncertainty, there are other sources of uncertainties that remain open for further investigation. Historical O2 profiles may have evolving precision and uncertainty that are difficult to replicate in a model-based study. For example, a Winkler titration performed on a Nansen bottle during the 1960s may have different precision than a more recent Winkler titration done on a Niskin bottle using amperometric or photometric end-point detection methods. Looking ahead, integration of autonomous float O2 data will pose challenges in terms of assessing uncertainties that are changing with the evolution of measurement techniques. Another important area of future investigation would be the uncertainties from natural variability. A recent modeling study (Fay et al., 2023) showed overlapping magnitudes of externally forced and internally generated O2 anomalies in the context of volcanic eruption. Quantification of natural variability is difficult to achieve using observations or a collection of single runs from multiple models. The best approach would be to use multi-model large ensembles with adequate ensemble members of randomized natural climate variability.

Code and data availability

The World Ocean Database (WOD) is available from the NOAA National Center for Environmental Information (, Boyer et al., 2018). CMIP6 model outputs are available from the Earth System Grid Federation (, WCRP, 2023). Optimal interpolation of WOD18 data is available from Ito (2023,


The supplement related to this article is available online at:

Author contributions

All the authors developed the ideas of this study, and TI designed the methods of data analysis and performed most of the analysis. The paper was mainly written by TI, with major comments and revisions by CAD, SM, MCL, HEG, ZW, and YT. All the authors contributed to the article and approved the submitted version.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Special issue statement

This article is part of the special issue “Low-oxygen environments and deoxygenation in open and coastal marine waters”. It is a result of the 53rd International Colloquium on Ocean Dynamics (3rd GO2NE Oxygen Conference), Liège, Belgium, 16–20 May 2022.


The oceanographic data that comprise the WOD have been acquired through many sources and projects, as well as from individual scientists. In addition, many international organizations such as the IODE/GODAR and WDS have facilitated data exchanges, which have provided much data to the WOD. Takamitsu Ito is grateful for funding support from the Department of Energy (grant no. DE-SC0021300), the National Science Foundation (grant no. OCE-2123546), and JSPS International Fellowship for Research in Japan. Shoshiro Minobe is supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI (grant no. JP18H04129). Yohei Takano is supported as part of the Energy Exascale Earth System Model project, funded by the US Department of Energy, Office of Science, Office of Biological and Environmental Research. The submission has been approved and assigned the LA-UR no. LA-UR-23-20603.

Financial support

This research has been supported by the National Science Foundation, the Directorate for Geosciences (grant no. OCE-2123546); the Department of Energy, Biological and Environmental Research (grant no. DE-SC0021300); and the Japan Society for the Promotion of Science (grant no. JP18H04129).

Review statement

This paper was edited by Marilaure Grégoire and reviewed by Lijing Cheng and one anonymous referee.


Bittig, H. C., Körtzinger, A., Neill, C., van Ooijen, E., Plant, J. N., Hahn, J., Johnson, K. S., Yang, B., and Emerson, S. R.: Oxygen Optode Sensors: Principle, Characterization, Calibration, and Application in the Ocean, Frontiers in Marine Science, 4, 429,, 2018. 

Boucher, O., Servonnat, J., Albright, A. L., Aumont, O., Balkanski, Y., Bastrikov, V., Bekki, S., Bonnet, R., Bony, S., Bopp, L., Braconnot, P., Brockmann, P., Cadule, P., Caubel, A., Cheruy, F., Codron, F., Cozic, A., Cugnet, D., D'Andrea, F., Davini, P., de Lavergne, C., Denvil, S., Deshayes, J., Devilliers, M., Ducharne, A., Dufresne, J.-L., Dupont, E., Éthé, C., Fairhead, L., Falletti, L., Flavoni, S., Foujols, M.-A., Gardoll, S., Gastineau, G., Ghattas, J., Grandpeix, J.-Y., Guenet, B., Guez, L. E., Guilyardi, E., Guimberteau, M., Hauglustaine, D., Hourdin, F., Idelkadi, A., Joussaume, S., Kageyama, M., Khodri, M., Krinner, G., Lebas, N., Levavasseur, G., Lévy, C., Li, L., Lott, F., Lurton, T., Luyssaert, S., Madec, G., Madeleine, J.-B., Maignan, F., Marchand, M., Marti, O., Mellul, L., Meurdesoif, Y., Mignot, J., Musat, I., Ottlé, C., Peylin, P., Planton, Y., Polcher, J., Rio, C., Rochetin, N., Rousset, C., Sepulchre, P., Sima, A., Swingedouw, D., Thiéblemont, R., Traore, A. K., Vancoppenolle, M., Vial, J., Vialard, J., Viovy, N., and Vuichard, N.: Presentation and Evaluation of the IPSL-CM6A-LR Climate Model, J. Adv. Model. Earth Sy., 12, e2019MS002010,, 2020. 

Boyer, T. P., Baranova, O. K., Coleman, C., Garcia, H. E., Grodsky, A., Locarnini, R. A., Mishonov, A. V., Paver, C. R., Reagan, J. R., Seidov, D., Smolyar, I. V., Weathers, K., and Zweng, M. M.: World Ocean Database 2018, NOAA Atlas NESDIS 87, edited by: Mishonov, A. V., Silver Spring, MD, 209, 2018 (data available at:, last access: 3 February 2024). 

Burrows, S. M., Maltrud, M., Yang, X., Zhu, Q., Jeffery, N., Shi, X., Ricciuto, D., Wang, S., Bisht, G., Tang, J., Wolfe, J., Harrop, B. E., Singh, B., Brent, L., Baldwin, S., Zhou, T., Cameron-Smith, P., Keen, N., Collier, N., Xu, M., Hunke, E. C., Elliott, S. M., Turner, A. K., Li, H., Wang, H., Golaz, J. C., Bond-Lamberty, B., Hoffman, F. M., Riley, W. J., Thornton, P. E., Calvin, K., and Leung, L. R.: The DOE E3SM v1.1 Biogeochemistry Configuration: Description and Simulated Ecosystem-Climate Responses to Historical Changes in Forcing, J. Adv. Model. Earth Sy., 12, e2019MS001766,, 2020. 

Carpenter, J. H.: THE ACCURACY OF THE WINKLER METHOD FOR DISSOLVED OXYGEN ANALYSIS1, Limnol. Oceanogr., 10, 135–140,, 1965. 

Deutsch, C., Ferrel, A., Seibel, B., Pörtner, H.-O., and Huey, R. B.: Climate change tightens a metabolic constraint on marine habitats, Science, 348, 1132–1135,, 2015. 

Dunne, J. P., Horowitz, L. W., Adcroft, A. J., Ginoux, P., Held, I. M., John, J. G., Krasting, J. P., Malyshev, S., Naik, V., Paulot, F., Shevliakova, E., Stock, C. A., Zadeh, N., Balaji, V., Blanton, C., Dunne, K. A., Dupuis, C., Durachta, J., Dussin, R., Gauthier, P. P. G., Griffies, S. M., Guo, H., Hallberg, R. W., Harrison, M., He, J., Hurlin, W., McHugh, C., Menzel, R., Milly, P. C. D., Nikonov, S., Paynter, D. J., Ploshay, J., Radhakrishnan, A., Rand, K., Reichl, B. G., Robinson, T., Schwarzkopf, D. M., Sentman, L. T., Underwood, S., Vahlenkamp, H., Winton, M., Wittenberg, A. T., Wyman, B., Zeng, Y., and Zhao, M.: The GFDL Earth System Model Version 4.1 (GFDL-ESM 4.1): Overall Coupled Model Description and Simulation Characteristics, J. Adv. Model. Earth Sy., 12, e2019MS002015,, 2020. 

Fay, A. R., McKinley, G. A., Lovenduski, N. S., Eddebbar, Y., Levy, M. N., Long, M. C., Olivarez, H. C., and Rustagi, R. R.: Immediate and Long-Lasting Impacts of the Mt. Pinatubo Eruption on Ocean Oxygen and Carbon Inventories, Global Biogeochem. Cy., 37, e2022GB007513,, 2023. 

Fukumori, I., Martel, F., and Wunsch, C.: The hydrography of the North Atlantic in the early 1980s. An atlas, Prog. Oceanogr., 27, 1–110, 1991. 

Garcia, H. E., Boyer, T. P., Levitus, S., Locarnini, R. A., and Antonov, J.: On the variability of dissolved oxygen and apparent oxygen utilization content for the upper world ocean: 1955 to 1998, Geophys. Res. Lett., 32, L09604,, 2005. 

Garcia, H. E., Weathers, K., Paver, C. R., Smolyar, I., Boyer, T. P., Locarnini, R. A., Zweng, M. M., Mishonov, A. V., Baranova, O. K., Seidov, D., and Reagan, J. R.: World Ocean Atlas 2018, NOAA National Centers for Environmental Information, Silver Springs, MD, 2019. 

Grégoire, M., Garçon, V., Garcia, H., Breitburg, D., Isensee, K., Oschlies, A., Telszewski, M., Barth, A., Bittig, H. C., Carstensen, J., Carval, T., Chai, F., Chavez, F., Conley, D., Coppola, L., Crowe, S., Currie, K., Dai, M., Deflandre, B., Dewitte, B., Diaz, R., Garcia-Robledo, E., Gilbert, D., Giorgetti, A., Glud, R., Gutierrez, D., Hosoda, S., Ishii, M., Jacinto, G., Langdon, C., Lauvset, S. K., Levin, L. A., Limburg, K. E., Mehrtens, H., Montes, I., Naqvi, W., Paulmier, A., Pfeil, B., Pitcher, G., Pouliquen, S., Rabalais, N., Rabouille, C., Recape, V., Roman, M., Rose, K., Rudnick, D., Rummer, J., Schmechtig, C., Schmidtko, S., Seibel, B., Slomp, C., Sumalia, U. R., Tanhua, T., Thierry, V., Uchida, H., Wanninkhof, R., and Yasuhara, M.: A Global Ocean Oxygen Database and Atlas for Assessing and Predicting Deoxygenation and Ocean Health in the Open and Coastal Ocean, Frontiers in Marine Science, 8, 724913,, 2021. 

Gruber, N.: Warming up, turning sour, losing breath: ocean biogeochemistry under global change, Philos. T. R. Soc. A, 369, 1980–1996,, 2011. 

Hajima, T., Watanabe, M., Yamamoto, A., Tatebe, H., Noguchi, M. A., Abe, M., Ohgaito, R., Ito, A., Yamazaki, D., Okajima, H., Ito, A., Takata, K., Ogochi, K., Watanabe, S., and Kawamiya, M.: Development of the MIROC-ES2L Earth system model and the evaluation of biogeochemical processes and feedbacks, Geosci. Model Dev., 13, 2197–2244,, 2020. 

Helm, K. P., Bindoff, N. L., and Church, J. A.: Observed decreases in oxygen content of the global ocean, Geophys. Res. Lett., 38, L23602,, 2011. 

IPCC: Changing Ocean, Marine Ecosystems, and Dependent Communities, in: The Ocean and Cryosphere in a Changing Climate: Special Report of the Intergovernmental Panel on Climate Change, Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, 447–588,, 2022. 

Ito, T.: Optimally interpolated dissolved oxygen based on the World Ocean Database 2018 and CMIP6 models, Zenodo [data set],, 2023. 

Ito, T., Minobe, S., Long, M. C., and Deutsch, C.: Upper ocean O2 trends: 1958–2015, Geophys. Res. Lett., 44, 4214–4223,, 2017. 

Ito, T., Long, M. C., Deutsch, C., Minobe, S., and Sun, D.: Mechanisms of Low-Frequency Oxygen Variability in the North Pacific, Global Biogeochem. Cy., 33, 110–124,, 2019. 

Keeling, R. F., Kortzinger, A., and Gruber, N.: Ocean Deoxygenation in a Warming World, Annu. Rev. Mar. Sci., 2, 199–229,, 2010. 

Levin, L. A.: Manifestation, Drivers, and Emergence of Open Ocean Deoxygenation, Annu. Rev. Mar. Sci., 10, 229–260,, 2018. 

Maurer, T. L., Plant, J. N., and Johnson, K. S.: Delayed-Mode Quality Control of Oxygen, Nitrate, and pH Data on SOCCOM Biogeochemical Profiling Floats, Frontiers in Marine Science, 8, 683207,, 2021. 

Mauritsen, T., Bader, J., Becker, T., Behrens, J., Bittner, M., Brokopf, R., Brovkin, V., Claussen, M., Crueger, T., Esch, M., Fast, I., Fiedler, S., Fläschner, D., Gayler, V., Giorgetta, M., Goll, D. S., Haak, H., Hagemann, S., Hedemann, C., Hohenegger, C., Ilyina, T., Jahns, T., Jimenéz-de-la-Cuesta, D., Jungclaus, J., Kleinen, T., Kloster, S., Kracher, D., Kinne, S., Kleberg, D., Lasslop, G., Kornblueh, L., Marotzke, J., Matei, D., Meraner, K., Mikolajewicz, U., Modali, K., Möbis, B., Müller, W. A., Nabel, J. E. M. S., Nam, C. C. W., Notz, D., Nyawira, S.-S., Paulsen, H., Peters, K., Pincus, R., Pohlmann, H., Pongratz, J., Popp, M., Raddatz, T. J., Rast, S., Redler, R., Reick, C. H., Rohrschneider, T., Schemann, V., Schmidt, H., Schnur, R., Schulzweida, U., Six, K. D., Stein, L., Stemmler, I., Stevens, B., von Storch, J.-S., Tian, F., Voigt, A., Vrese, P., Wieners, K.-H., Wilkenskjeld, S., Winkler, A., and Roeckner, E.: Developments in the MPI-M Earth System Model version 1.2 (MPI-ESM1.2) and Its Response to Increasing CO2, J. Adv. Model. Earth Sy., 11, 998–1038,, 2019.  

Nakanowatari, T., Ohshima, K. I., and Wakatsuchi, M.: Warming and oxygen decrease of intermediate water in the northwestern North Pacific, originating from the Sea of Okhotsk, 1955–2004, Geophys. Res. Lett., 34, L04602,, 2007. 

Oschlies, A., Duteil, O., Getzlaff, J., Koeve, W., Landolfi, A., and Schmidtko, S.: Patterns of deoxygenation: sensitivity to natural and anthropogenic drivers, Philos. T. R. Soc. A, 375, 20160325,, 2017. 

Oschlies, A., Brandt, P., Stramma, L., and Schmidtko, S.: Drivers and mechanisms of ocean deoxygenation, Nat. Geosci., 11, 467–473,, 2018. 

Pörtner, H. O. and Farrell, A. P.: Physiology and Climate Change, Science, 322, 690–692,, 2008. 

Schmidtko, S., Stramma, L., and Visbeck, M.: Decline in global oceanic oxygen content during the past five decades, Nature, 542, 335–339,, 2017. 

Seland, Ø., Bentsen, M., Olivié, D., Toniazzo, T., Gjermundsen, A., Graff, L. S., Debernard, J. B., Gupta, A. K., He, Y.-C., Kirkevåg, A., Schwinger, J., Tjiputra, J., Aas, K. S., Bethke, I., Fan, Y., Griesfeller, J., Grini, A., Guo, C., Ilicak, M., Karset, I. H. H., Landgren, O., Liakka, J., Moseid, K. O., Nummelin, A., Spensberger, C., Tang, H., Zhang, Z., Heinze, C., Iversen, T., and Schulz, M.: Overview of the Norwegian Earth System Model (NorESM2) and key climate response of CMIP6 DECK, historical, and scenario simulations, Geosci. Model Dev., 13, 6165–6200,, 2020. 

Stramma, L., Oschlies, A., and Schmidtko, S.: Mismatch between observed and modeled trends in dissolved upper-ocean oxygen over the last 50 yr, Biogeosciences, 9, 4045–4057,, 2012. 

Swart, N. C., Cole, J. N. S., Kharin, V. V., Lazare, M., Scinocca, J. F., Gillett, N. P., Anstey, J., Arora, V., Christian, J. R., Hanna, S., Jiao, Y., Lee, W. G., Majaess, F., Saenko, O. A., Seiler, C., Seinen, C., Shao, A., Sigmond, M., Solheim, L., von Salzen, K., Yang, D., and Winter, B.: The Canadian Earth System Model version 5 (CanESM5.0.3), Geosci. Model Dev., 12, 4823–4873,, 2019. 

Vaquer-Sunyer, R. and Duarte, C. M.: Thresholds of hypoxia for marine biodiversity, P. Natl. Acad. Sci. USA, 105, 15452–15457,, 2008. 

WCRP: CMIP6 model outputs,, WCRP, last access: March 2023. 

Wunsch, C.: The Ocean Circulation Inverse Problem, Cambridge University Press, Cambridge,, 1996. 

Short summary
This study aims to estimate how much oceanic oxygen has been lost and its uncertainties. One major source of uncertainty comes from the statistical gap-filling methods. Outputs from Earth system models are used to generate synthetic observations where oxygen data are extracted from the model output at the location and time of historical oceanographic cruises. Reconstructed oxygen trend is approximately two-thirds of the true trend.
Final-revised paper