GOBAI-O₂: temporally and spatially resolved fields of ocean interior dissolved oxygen over nearly 2 decades

For about 2 decades, oceanographers have been installing oxygen sensors on Argo profiling floats to be deployed throughout the world ocean, with the stated objective of better constraining trends and variability in the ocean's inventory of oxygen. Until now, measurements from these Argo-float-mounted oxygen sensors have been mainly used for localized process studies on air–sea oxygen exchange, upper-ocean primary production, biological pump efficiency, and oxygen minimum zone dynamics. Here, we present a new four-dimensional gridded product of ocean interior oxygen, derived via machine learning algorithms trained on dissolved oxygen observations from Argo-float-mounted sensors and discrete measurements from ship-based surveys and applied to temperature and salinity fields constructed from the global Argo array. The data product is called GOBAI-O₂, which stands for Gridded Ocean Biogeochemistry from Artificial Intelligence – Oxygen (Sharp et al., 2022; https://doi.org/10.25921/z72m-yz67); it covers 86 % of the global ocean area on a 1^∘ × 1^∘ (latitude × longitude) grid, spans the years 2004–2022 with a monthly resolution, and extends from the ocean surface to a depth of 2 km on 58 levels. Two types of machine learning algorithms – random forest regressions and feed-forward neural networks – are used in the development of GOBAI-O₂, and the performance of those algorithms is assessed using real observations and simulated observations from Earth system model output. Machine learning represents a relatively new method for gap filling ocean interior biogeochemical observations and should be explored along with statistical and interpolation-based techniques. GOBAI-O₂ is evaluated through comparisons to the oxygen climatology from the World Ocean Atlas, the mapped oxygen product from the Global Ocean Data Analysis Project and to direct observations from large-scale hydrographic research cruises. Finally, potential uses for GOBAI-O₂ are demonstrated by presenting average oxygen fields on isobaric and isopycnal surfaces, average oxygen fields across vertical–meridional sections, climatological seasonal cycles of oxygen averaged over different pressure layers, and globally integrated time series of oxygen. GOBAI-O₂ indicates a declining trend in the oxygen inventory in the upper 2 km of the global ocean of 0.79 ± 0.04 % per decade between 2004 and 2022.

1 Introduction

The inventory of dissolved oxygen in the global ocean has been declining over recent decades and is projected to continue declining through the current century (Keeling et al., 2010; Breitburg et al., 2018; Bindoff et al., 2019; Stramma and Schmidtko, 2019; Limburg et al., 2020), leading to detrimental consequences for aerobic marine organisms (Pörtner and Farrell, 2008; Sampaio et al., 2021) and changes in biogeochemical cycles, potentially triggering important climatological feedbacks (Gruber, 2004; Berman-Frank et al., 2008). Historical deoxygenation has been inferred from analyses of globally distributed observations (Helm et al., 2011; Schmidtko et al., 2017; Ito et al., 2017) and has been reproduced in Earth system model (ESM) reconstructions (Bopp et al., 2013; Frölicher et al., 2009; Kwiatkowski et al., 2020). Global observational studies have generally indicated a greater degree of deoxygenation than model studies over recent decades, indicating that ESMs may misrepresent the sensitivities of the physical and biological processes leading to deoxygenation, which has implications for the reliability of future projections (Oschlies et al., 2017, 2018; Stramma and Schmidtko, 2021). Model studies, however, are based on gridded output that is continuously resolved in space and time, whereas observational studies rely on the interpolation of measurements from discrete bottle samples and/or profiling sensors. These observational datasets have significant spatiotemporal gaps and may not robustly represent global deoxygenation trends.

Discrete measurements of the dissolved oxygen concentration ([O₂]) are typically made using Winkler titrations (Winkler, 1888; Carpenter, 1965; Langdon, 2010), which are also used to calibrate measurements from electrode (or more recently sometimes optical) dissolved oxygen sensors mounted on conductivity–temperature–depth (CTD) profilers. Globally distributed [O₂] observations from discrete measurements and CTD profilers have been provided by hydrographic programs like the World Ocean Circulation Experiment (WOCE); the Climate and Ocean Variability, Predictability and Change (CLIVAR) program; and the Global Ocean Ship-Based Hydrographic Investigations Program (GO-SHIP). Data from these programs are publicly available and are conveniently compiled into databases such as the World Ocean Database (WOD; Boyer et al., 2018) and the Global Ocean Data Analysis Project (GLODAP; Lauvset et al., 2022b). Although unprecedented spatial coverage is provided by global hydrographic programs, the decadal-scale temporal resolution of WOCE, CLIVAR, and GO-SHIP data precludes robust analyses of year-to-year and/or seasonal variability in [O₂].

Since the mid-2000s, approximately 1800 Argo floats equipped with oxygen sensors have been deployed. Argo floats profile the upper ∼2000 m of the water column every ∼10 d. Many oxygen-sensor-equipped Argo floats have been deployed as part of regional arrays such as the Southern Ocean Carbon and Climate Observations and Modeling (SOCCOM) project and the North Atlantic Aerosols and Marine Ecosystems Study (NAAMES). More recently, the push for a global biogeochemical Argo array has spurred the deployment of oxygen-sensor-equipped Argo floats into more sparsely sampled ocean regions (Johnson and Claustre, 2016; Claustre et al., 2020). As more floats have been deployed, improvements have been made to sensor calibration, data adjustments, and quality control. Notably, pre-deployment drift corrections (D'Asaro and McNeil, 2013; Johnson et al., 2015; Bittig and Körtzinger, 2015; Bushinsky et al., 2016; Drucker and Riser, 2016; Nicholson and Feen, 2017), climatology-based calibrations (Takeshita et al., 2013), calibrations via in-air oxygen measurements (Körtzinger et al., 2005; Fiedler et al., 2013; Bittig and Körtzinger, 2015; Johnson et al., 2015; Bushinsky et al., 2016), post-deployment drift corrections (Johnson et al., 2017; Bittig et al., 2018a), and established procedures for delayed-mode quality control (Maurer et al., 2021) have substantially reduced the uncertainty and increased the reproducibility of optode-based [O₂] measurements on Argo floats.

From the time it began, the Argo-Oxygen program (now oxygen is a measured variable under the Biogeochemical Argo program) intended to document ocean deoxygenation, predict and assess anoxic and hypoxic events, and determine seasonal to interannual changes in export production (Gruber et al., 2010). To date, these goals have been achieved primarily on a regional scale. For example, [O₂] measurements from biogeochemical Argo floats have been used to examine ventilation and air–sea exchange of oxygen in the Southern Ocean (Bushinsky et al., 2017) and during deep water formation in the subpolar North Atlantic (Körtzinger et al., 2004; Piron et al., 2016, 2017; Wolf et al., 2018); denitrification and the spatial extent of the oxygen minimum zone (OMZ) in the Bay of Bengal (Sarma and Udaya Bhaskar, 2018; Johnson et al., 2019; Udaya Bhaskar et al., 2021); and carbon production and export in the Pacific Ocean (Bushinsky and Emerson, 2015, 2018; Yang et al., 2017), Southern Ocean (Stukel and Ducklow, 2017; Arteaga et al., 2019), and North Atlantic Ocean (Alkire et al., 2012; Estapa et al., 2019). Recently, in an early global-scale analysis of [O₂] from the Argo array, Johnson and Bif (2021) used the diel cycle of oxygen measured by the ocean-wide array of biogeochemical Argo floats to constrain net primary production in the surface ocean.

With the work presented here, we seek to capitalize on the collective efforts of global hydrographic programs, Biogeochemical Argo (BGC Argo), and Core Argo (Johnson et al., 2022) to create a novel data product: a four-dimensional monthly record of dissolved oxygen in the global ocean. We combine autonomous observations of [O₂] from BGC Argo floats with discrete observations of [O₂] from hydrographic cruises in the GLODAP database to create a dataset with an extensive spatial and temporal resolution. With this dataset, we train machine learning algorithms on ocean interior predictor variables co-located with [O₂] observations; evaluate those algorithms using real and simulated data; and apply the algorithms to gridded ocean interior predictor variables mapped from Core Argo to produce a gridded [O₂] data product at a monthly resolution from 2004 to 2022, on 58 pressure levels in the upper 2 km of the ocean, and on a near-global 1^∘ × 1^∘ (latitude × longitude) grid.

In this paper, we present the four-dimensional gridded [O₂] product, which we call GOBAI-O₂: Gridded Ocean Biogeochemistry from Artificial Intelligence – Oxygen (Sharp et al., 2022; https://doi.org/10.25921/z72m-yz67). Artificial intelligence (AI) is a broad term describing computerized systems that are able to recognize patterns, make decisions, and otherwise perform tasks previously reserved only for humans; GOBAI-O₂ is built using machine learning (ML), which is a subfield of AI that is focused on training, understanding, and applying algorithms that leverage data to artificially learn and reproduce patterns. We introduce GOBAI-O₂ by analyzing spatial patterns, seasonal cycles, and decadal variability. We also describe the process for creating GOBAI-O₂, show the results of evaluation exercises, assess uncertainty in the gridded [O₂] fields, and compare the data product to other gridded datasets and discrete measurements. GOBAI-O₂ represents the first step in leveraging the emerging global array of BGC Argo floats to produce spatially resolved, time-varying snapshots of global ocean biogeochemical distributions in near-real time. It also emphasizes regions that are good candidates for new observational assets; particularly where gaps in observational coverage coincide with high background variability in [O₂]. Critically, GOBAI-O₂ can be used to address the goals of the Argo-Oxygen program set by Gruber et al. (2010) over a decade ago, providing regional and global insight into ocean deoxygenation and hypoxia on timescales ranging from a few months to multiple years.

2 Methods

2.1 Data sources and processing

Hydrographic cruise data were obtained from the GLODAP version 2022 data product (GLODAPv2.2022; Key et al., 2015; Olsen et al., 2016; Lauvset et al., 2022a, b). GLODAPv2.2022 provides quality-controlled data from throughout the entire water column obtained via discrete analyses of more than 1.4 million water samples collected on 1085 research cruises. Discrete Winkler titration data were chosen rather than CTD oxygen profiles due to the issues with the quality of calibration of a subset of CTD oxygen measurements and the relatively coarse vertical resolution of the final GOBAI-O₂ product, which would not benefit from the high vertical resolution of CTD profiles. Data from GLODAP were chosen rather than data from the WOD or any other database due to the high degree of quality control applied to GLODAP data. Dissolved oxygen is the most represented biogeochemical variable in GLODAPv2.2022, with more than 1.2 million data points from 991 research cruises. Data from GLODAPv2.2022 were filtered to retain only samples collected after 1 January 2004, from 0 to 2500 dbar (decibars), and with a quality control flag of 1 (meaning the data were manually inspected) and quality flags of 2 (good) for both salinity and [O₂]. Temperature is not assigned either flag and is assumed to be of sufficient quality if it is reported (Lauvset et al., 2022a). This filtering left 450 032 data points from 21 513 unique profiles from 393 total cruises (red points in Fig. 1).

Figure 1 Discrete profile locations from oxygen-sensor-equipped Argo floats (blue) and GLODAPv2.2022 cruises (red) from 1 January 2004 to 13 June 2023. Data from these profiles were binned and used to train ML algorithms to estimate [O₂] in each of seven regions: the Atlantic Ocean (Atl.), Pacific Ocean (Pac.), Indian Ocean (Ind.), Arctic Ocean (Arc.), Mediterranean Sea (Med.), northern section of the Southern Ocean (N. Sou.), and southern section of the Southern Ocean (S. Sou.). Overlapping areas between regions are shown in grey (Ovrlp.), where [O₂] estimates are made by taking distance-weighted averages of outputs from two regional ML algorithms. The regional boundaries are presented in numerical form in Table B1.

Float data were obtained from synthetic profile (“Sprof”) files (Bittig et al., 2022) stored in the Argo Global Data Assembly Centres (GDACs) via the OneArgo-Mat toolbox (Frenzel et al., 2022) for MATLAB (MathWorks). At the time data were obtained (13 June 2023), the Argo GDACs contained data from about 1800 floats equipped with [O₂] sensors. Float data were filtered to retain only delayed-mode-adjusted data with quality flags of 1 (good), 2 (probably good), or 8 (interpolated/extrapolated) for pressure, temperature, salinity, and [O₂]. This filtering step ensured that float data had been manually reviewed by a data manager and assigned an appropriate quality flag. This filtering left 27 832 192 data points from 138 180 unique profiles from 1022 total floats (blue points in Fig. 1). Of the float profiles, 51.4 % were quality-controlled by comparisons to climatologies (World Ocean Atlas, WOA, or Commonwealth Scientific and Industrial Research Organisation Atlas of Regional Seas, CARS), 30.3 % were quality-controlled using in-air oxygen measurements, 7.0 % were quality-controlled by comparisons to subsurface measurements (WOD, OMZ assumed to have zero oxygen, or CTD profile upon deployment), 5.3 % were quality-controlled using the in situ optode calibration of Drucker and Riser (2016), 3.3 % were quality-controlled via another method, 1.9 % were not categorized, and the remaining 0.9 % were not adjusted.

The discrete temperature, salinity, and [O₂] data obtained from GLODAPv2.2022 and the Argo GDACs are archived online (see Appendix C; https://doi.org/10.5281/zenodo.7747237, Sharp, 2023b). To ensure that the trained machine learning algorithms were not biased toward BGC Argo float data, which (in their native format) have a higher vertical resolution than GLODAP data, each profile was interpolated to, at most, 58 standard pressure levels (the same pressure levels on which the final GOBAI-O₂ data product is provided). Interpolated temperature, salinity, and [O₂] data from each source are also archived online (see Appendix C; https://doi.org/10.5281/zenodo.7747237, Sharp, 2023b). After interpolation, the total number of GLODAP data points used for algorithm training increased to 1 096 324 and the total number of Argo float data points used for algorithm training decreased to 6 635 749. Co-located, interpolated GLODAP and BGC Argo profiles that fell within the same $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ monthly, depth-dependent grid cells were compared for internal consistency. The float [O₂] values were adjusted according to the procedure in Appendix D to remove the small global discrepancy between co-located ship and float measurements in order to ensure internal consistency between the two datasets.

BGC Argo float and GLODAP cruise data were combined into a single dataset after this bias adjustment, which will be referred to as the “combined dataset” from here on. The combined dataset was grouped into seven overlapping regions (Fig. 1, Table B1). This grouping was intended to account implicitly for similar physical–biogeochemical relationships within large ocean regions and to reduce the computational burden of the machine learning (ML) algorithm fits described below. The regions were initially chosen to imitate the biomes presented by Fay and McKinley (2014), and they were then expanded to relatively large regions bound either by land masses or by overlapping boundaries along constant lines of latitude. The number of profiles made within each $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ box by either a discrete ship cast or Argo float (Fig. A1) provides a measure of the temporal resolution of the combined dataset in addition to the spatial distribution shown in Fig. 1.

Gridded temperature and salinity data to which the trained algorithms were applied were obtained from the latest version of the Roemmich and Gilson (2009) (RG09) Argo Climatology (https://sio-argo.ucsd.edu/RG_Climatology.html, last access: 12 January 2023). The RG09 climatology is an upper-ocean (0–2000 dbar) gridded temperature and salinity product constructed exclusively from Argo observations. Long-term (2004–2018) mean fields of temperature and salinity are provided on 58 pressure levels, along with monthly anomaly fields on each of those pressure levels from 2004 to the present day. The most recent major update of the RG09 climatology was made in 2019, and new monthly anomaly fields are provided in near-real time between major updates. Monthly gridded temperature and salinity were calculated from the RG09 long-term mean and monthly anomaly fields (Fig. A2), and they were then used for the creation of the gridded [O₂] product discussed below.

Output from the National Oceanic and Atmospheric Administration (NOAA) Geophysical Fluid Dynamics Laboratory's Earth System Model Version 4 (GFDL-ESM4; Dunne et al., 2020) was used to assess algorithm performance. Model output was downloaded from the World Climate Research Programme database (https://esgf-node.llnl.gov/projects/cmip6/, last access: 8 April 2022), which hosts data from models participating in Phase 6 of the Coupled Model Intercomparison Project (CMIP6). Potential temperature, practical salinity, and [O₂] were downloaded to coincide with available ocean interior observations (Fig. A3). Historical outputs (2004–2014) and projected outputs under Shared Socioeconomic Pathway (SSP) 2-4.5 (2015–2022) were combined to cover the time period over which observations were available. A spatial mask was applied to retain only GFDL-ESM4 grid cells with corresponding temperature and salinity values in the RG09 climatology, because that is the final grid on which GOBAI-O₂ is produced.

2.2 Algorithm training

The combined dataset was used to train ML algorithms for each region to estimate [O₂] from absolute salinity, conservative temperature, potential density anomaly, hydrostatic pressure, bottom depth, and additional spatiotemporal information to allow for geographic, seasonal, and interannual variation (see Table 1). Although biology is not explicitly accounted for in the ML algorithms, Giglio et al. (2018) demonstrate that, with an appropriately distributed dataset, the inclusion of spatiotemporal variables in algorithm training can implicitly accommodate biological processes.

Table 1 Predictor variables used to train random forest regressions and feed-forward neural networks to predict [O₂]. Unitless quantities are denoted using “–”.

Absolute salinity (S_A) was calculated from practical salinity (S_P), hydrostatic pressure (P), latitude, and longitude. Conservative temperature (θ) was calculated from in situ temperature (T), S_A, and P. Potential density anomaly (σ_θ) was calculated from S_A and θ. These calculations were made using the Gibbs-SeaWater (GSW) Oceanographic Toolbox for MATLAB (McDougall and Barker, 2011). As was done by Carter et al. (2021), longitude was transformed into two separate predictors: $\cos (\text{Longitude}-\mathrm{20}{}^{\circ }\text{E})$ and $\cos (\text{Longitude}-\mathrm{110}{}^{\circ }\text{E})$. Cosine functions were applied to maintain the cyclical nature of longitude as a predictor, and offsets of 20 and 110^∘ E were intended to shift regions where the cosine function has minimum explanatory power over landmasses. Bottom depth was determined by matching each observational location with the corresponding bathymetry from the ETOPO2v2 global relief model (NOAA National Geophysical Data Center, 2006).

Two types of ML algorithms were trained: feed-forward neural networks (FNNs; Beale et al., 2023) and random forest regressions (RFRs; Breiman, 2001). Each algorithm type was trained on the input variables given in Table 1 to produce estimates of [O₂] (Fig. A4). Three separate FNNs were trained for each of the seven basins shown in Fig. 1, with an average of the three taken to obtain one equally weighted FNN result. The FNNs were constructed using the “feedforwardnet” function and trained using the “train” function, both from Version 14.4 of the Deep Learning Toolbox for MATLAB (R2022a). Each FNN was trained using a Levenberg–Marquardt algorithm, with 15 % of the data reserved for testing the network during training steps. Each FNN had two hidden layers, with the following combinations of neurons in the first and second layer, respectively: 20 and 10, 15 and 15, and 10 and 20. One RFR was trained for each of the seven basins shown in Fig. 1. RFRs are ensembles of decision trees, each created with a bootstrapped version of the full dataset chosen randomly with replacement. Each RFR consisted of 600 trees, a minimum leaf size of 10, and 6 of the 11 predictors used for each decision split. These parameters were chosen after some trial and error to strike a balance between computational efficiency and algorithm performance. The MATLAB “treebagger” function was used to train RFRs.

In areas where two regions overlap (see Fig. 1), weighted averages of [O₂] estimates were calculated in overlapping grid cells from each regional algorithm. These averages were weighted by distance from the center latitude line of the overlapping area (e.g., a point at 33^∘ S in the overlapping area between the N. Sou. region, whose northern border extends to 25^∘ S, and the Atl. region, whose southern border extends to 35^∘ S, would be calculated as $\left[{\mathrm{O}}_{\mathrm{2}}\right]=\mathrm{0.8}[{\mathrm{O}}_{\mathrm{2}}{]}_{\text{N.Sou.}}+\mathrm{0.2}[{\mathrm{O}}_{\mathrm{2}}{]}_{\text{Atl.}}$). Overlapping areas were used to mitigate discontinuities at the boundaries between regions in the final gridded product.

The average of FNN and RFR estimates (ENS, for ensemble average) was used as the [O₂] estimate for a given set of input data. This ensemble-averaging procedure was implemented due to insights from previous work showing that averaging the outputs of multiple ML algorithms or linear regression models often outperforms the output from just one approach on its own (Gregor et al., 2017, 2019; Bittig et al., 2018b; Carter et al., 2021; Djeutchouang et al., 2022), likely due to complementary strengths and weaknesses of each approach. For this work, any especially erroneous result from either the FNN or RFR should be mitigated by better results from the other algorithm.

2.3 Algorithm evaluation

We performed two exercises to evaluate the effectiveness of the ML algorithms used to estimate [O₂]. The first exercise involved training separate evaluation algorithms (RFR_Data-Eval and FNN_Data-Eval), as described in Sect. 2.2, using a subset of the observational dataset for training while reserving the remaining subset for assessment. For this exercise, data were split randomly into training (80 %) and assessment (20 %) groups; this split was made according to measurement platform (cruise or float; see Fig. A5) to ensure that inherent correlations among the data points from a single cruise or float did not contribute to the apparent effectiveness of each ML algorithm. Then, [O₂] values from the subset of reserved assessment data were compared to estimates of [O₂] from RFR_Data-Eval, FNN_Data-Eval, and the ensemble average of the two (ENS_Data-Eval). This exercise was intended to evaluate the ability of the ML algorithms to reproduce measured data that were not involved in algorithm training (Sect. 3.1.1).

The second exercise involved training evaluation algorithms (RFR_ESM4-Eval and FNN_ESM4-Eval) using synthetic “profiles” extracted from gridded GFDL-ESM4 output at the times and locations for which observational data were available and then assessing the evaluation algorithms using spatially and temporally continuous monthly GFDL-ESM4 output from 2004 through 2022. For this exercise, synthetic profiles for algorithm training were defined by matching the latitude, longitude, month, and year of each available grid cell from the binned observational dataset with the corresponding GFDL-ESM4 output. This resulted in 75 879 synthetic profiles for algorithm training. The RFR_ESM4-Eval and FNN_ESM4-Eval algorithms were trained as described in Sect. 2.2 with the synthetic training data and then used to produce [O₂] estimates for the complete model output. These [O₂] estimates from RFR_ESM4-Eval, FNN_ESM4-Eval, and an ensemble average of the two (ENS_ESM4-Eval) were compared to [O₂] values from the full GFDL-ESM4 output fields at the grid-cell level. This exercise was intended to evaluate the ability of the ML algorithms to estimate [O₂] in a spatiotemporally resolved Earth system model environment when limited to training data representative of the available collection of ocean oxygen observations (Sect. 3.1.2). The four-dimensional field of [O₂] from ENS_ESM4-Eval that represents a reconstruction of the GFDL-ESM4 environment, which we refer to as GOBAI-O₂-ESM4, can also be used as an analogue for how well GOBAI-O₂ (trained on real observational data; Sect. 2.4) might represent [O₂] variability in the real-world environment. For this reason, the four-dimensional field of differences between GOBAI-O₂-ESM4 and GFDL-ESM4 output were used to inform the evaluation of GOBAI-O₂ uncertainty (Sects. 2.5 and 3.2.4). Additionally, we quantified global means, seasonal-cycle amplitudes, long-term trends, and interannual variabilities in [O₂] across different pressure layers of GOBAI-O₂-ESM4. To evaluate the performance of GOBAI-O₂-ESM4 on a global scale, these metrics are compared to the same metrics for the spatiotemporally resolved GFDL-ESM4 output and subsampled grid cells in GFDL-ESM4 corresponding to observational data coverage (Sect. 3.1.2). Comparisons of global means from GOBAI-O₂-ESM4 to GFDL-ESM4 are also used to approximate uncertainty in oxygen inventories for the assessment of trends (Sect. 3.2.3).

2.4 Creation of GOBAI-O₂

FNNs and RFRs for each of the seven regions shown in Fig. 1 were trained with the full combined dataset, using the predictor variables shown in Table 1, with [O₂] as a target variable. Then, the FNNs and RFRs were applied to S_A, θ, and σ_θ calculated from RG09 temperature and salinity fields, along with spatiotemporal information from RG09 grid cells. Weighted averages were calculated where regions overlapped, and ensemble averages (ENS) were calculated from the FNN and RFR estimates. One exception occurred in the subsurface of the north Arabian Sea (above 21^∘ N, between 900 and 1412.5 dbar), where erroneously high ENS [O₂] values caused by a clearly nonphysical feature in RFR [O₂] values were replaced by FNN [O₂] values only. This produced a monthly gridded [O₂] product in the upper 2 km of the ocean on a global grid from January 2004 to December 2022, i.e., GOBAI-O₂ (Sharp et al., 2022; https://doi.org/10.25921/z72m-yz67; Sect. 3.2.1–3.2.3). GOBAI-O₂ was compared to gridded climatological oxygen fields from the 2018 World Ocean Atlas (WOA18; Garcia et al., 2019; Sect. 3.2.5), the GLODAP mapped data product (Lauvset et al., 2016; Sect. 3.2.5), and discrete measurements of oxygen from select cruises between 2004 and 2022 (Sect. 3.2.6).

2.5 Uncertainty estimation

Similar to previous studies that have estimated uncertainty in observation-based biogeochemical data products (e.g., Landschützer et al., 2014; Gregor and Gruber, 2021; Keppler et al., 2020, 2023), we combine uncertainty from three separate sources – measurement, gridding, and algorithm – to estimate uncertainty in GOBAI-O₂ (Sect. 3.2.4).

Measurement uncertainty (u([O₂])_meas.) is attributable to the [O₂] observations themselves. For this quantity, gridded [O₂] from GOBAI-O₂ is multiplied by 3 %, which is an estimate based on a few factors: the consistency of the GLODAPv2.2022 dataset, the accuracy of BGC Argo observations, the relative proportion of each data source, and the recognized issue of float oxygen sensor response time. The nominal value for the consistency of the GLODAPv2.2022 cruise dataset is stated to be 1 % (Lauvset et al., 2022a). The approximate accuracy of BGC Argo float observations is estimated to be about 3 % of surface [O₂], which is about 4 % to 5 % of average ocean [O₂], for floats quality-controlled by climatological comparisons (Takeshita et al., 2013) and about 3 µmol kg⁻¹ (Johnson et al., 2017; Maurer et al., 2021), which is about 2 % of average ocean [O₂], for floats quality-controlled by in-air measurements. We choose a value for measurement uncertainty closer to the float accuracy estimates because float observations outweigh ship observations by about six to one in the dataset used to construct GOBAI-O₂. The estimates for the accuracy of float [O₂] sensors may be somewhat optimistic, especially when the floats are crossing large vertical oxygen gradients (i.e., steep oxyclines). This is because the response time of oxygen optodes is not instantaneous (Bittig et al., 2014, 2018a; Bittig and Körtzinger, 2017). Without response-time corrections, which are not widely applied to the BGC Argo oxygen dataset, there is the potential for systematic biases where float profiles traverse steep oxyclines, such as in the eastern tropical Pacific Ocean and North Indian Ocean.

Gridding uncertainty (u([O₂])_grid.) is attributable to using a single [O₂] value to represent a four-dimensional box that is coarser in time and space than the resolution of many processes that influence [O₂]. We estimate gridding uncertainty by (1) binning the combined GLODAP and Argo observational dataset to grid cells equal in size to the RG09 grid cells; (2) calculating the standard deviation among the observations in cells with more than 10 observations (Fig. A6); (3) fitting a multivariate polynomial regression relating those standard deviations to pressure, potential density anomaly, and bottom depth; and (4) applying that regression to the RG09 grid to compute an estimated standard deviations (i.e., gridding uncertainty) in each grid cell.

Algorithm uncertainty (u([O₂])_alg.) is attributable to the ML algorithms that estimate [O₂] on the RG09 grid. We estimate algorithm uncertainty using the four-dimensional field of absolute differences between [O₂] from GFDL-ESM4 model output and GOBAI-O₂-ESM4, determined from the GFDL-ESM4 algorithm evaluation exercise described in Sect. 2.3.

The three uncertainty sources were combined in quadrature (assuming independence) to calculate a combined uncertainty estimate for each gridded [O₂] value in GOBAI-O₂ (u([O₂])_tot.):

$$\begin{array}{}\text{(1)}& u\left(\right[{\mathrm{O}}_{\mathrm{2}}]{)}_{\text{tot.}}=\sqrt{u\left(\right[{\mathrm{O}}_{\mathrm{2}}]{)}_{\text{meas.}}^{\mathrm{2}}+u\left(\right[{\mathrm{O}}_{\mathrm{2}}]{)}_{\text{grid.}}^{\mathrm{2}}+u(\left[{\mathrm{O}}_{\mathrm{2}}\right]{)}_{\text{alg.}}^{\mathrm{2}}}.\end{array}$$

3 Results and discussion

3.1 Algorithm evaluation

The evaluation exercises indicated that the ML algorithms trained on the combined GLODAP and Argo observational dataset were effective in their ability to estimate [O₂] and reconstruct seasonal to decadal variability in the global oxygen inventory. Mean offsets ($\mathrm{\Delta }\left[{\mathrm{O}}_{\mathrm{2}}\right]=[{\mathrm{O}}_{\mathrm{2}}{]}_{\text{obs/mod}}-[{\mathrm{O}}_{\mathrm{2}}{]}_{\text{est}}$) and root-mean-square differences (RMSDs) between [O₂] from direct measurements ([O₂]_obs) or GFDL-ESM4 output ([O₂]_mod) and [O₂] estimated from ML algorithms ([O₂]_est) were determined as an assessment of the ability of the algorithms to estimate [O₂] at a grid-cell level (Tables 2, B2, B3, B4; Fig. 2). Mean Δ[O₂] and RMSD determined using [O₂]_est from the ESPER-Mixed model (Carter et al., 2021) – an average of predictions from a neural network and moving window multiple linear regression trained on GLODAPv2.2020 data – were also determined as a point of comparison for the observational-data-based validation test (Tables 2, B4; Fig. A7). In the case of the GFDL-ESM4-based validation test, metrics to summarize means, amplitudes, trends, and variability in integrated mean [O₂] values were determined to demonstrate the ability of the GOBAI-O₂ method to capture seasonal- to decadal-scale variability in oxygen at the global scale (Table 3, Fig. 3). The results of each evaluation exercise are discussed in more detail in the following sections.

Table 2 Regional and global error statistics (mean Δ[O₂] and RMSD) for evaluation exercises using the ensemble average (ENS_Data-Eval) of the FNN_Data-Eval and RFR_Data-Eval algorithms trained on a subset of data from the combined GLODAP and Argo observational dataset and tested with a separate subset of withheld data (left side of the table) or using the ensemble average (ENS_ESM4-Eval) of the FNN_ESM4-Eval and RFR_ESM4-Eval algorithms trained on a subset of output from GFDL-ESM4 (corresponding to locations of available Argo and GLODAP data) and tested using the full field of GFDL-ESM4 output (right side of the table). Error statistics calculated using the ESPER-Mixed model are also shown for comparison to the data-based test. The numbers of data points used in the training and assessment of each algorithm are shown.

Figure 2 Two-dimensional histograms showing offsets of measured versus estimated oxygen ($\mathrm{\Delta }\left[{\mathrm{O}}_{\mathrm{2}}\right]=[{\mathrm{O}}_{\mathrm{2}}{]}_{\text{obs}}-[{\mathrm{O}}_{\mathrm{2}}{]}_{\text{est}}$) for (a, b) withheld observational data and (d, e) modeled versus estimated oxygen ($\mathrm{\Delta }\left[{\mathrm{O}}_{\mathrm{2}}\right]=[{\mathrm{O}}_{\mathrm{2}}{]}_{\text{mod}}-[{\mathrm{O}}_{\mathrm{2}}{]}_{\text{est}}$) for GFDL-ESM4 model output as a function of (a, d) [O₂]_est and (b, e) depth in the water column. Offsets are binned into cells that are 2.5 µmol kg⁻¹ tall in terms of Δ[O₂] and (a, d) 5 µmol kg⁻¹ wide in terms of [O₂]_est or equivalent in width to (b) the interpolated pressure levels of the data or (e) the vertical resolution of GFDL-ESM4 grid cells. The frequency of offsets that fall into a given bin is shown on a logarithmic scale, de-emphasizing the significant clustering around Δ[O₂]=0 in favor of showing the few outliers. Absolute Δ[O₂] values averaged over depth and time for 1^∘ × 1^∘ (latitude × longitude) grid cells in the global ocean for (c) withheld observational data and (f) GFDL-ESM4 model output.

Table 3 Statistics representing the mean values, seasonal-cycle amplitudes, long-term trends, and interannual variabilities in [O₂] from the GFDL-ESM4 model, a reconstruction of [O₂] fields from GFDL-ESM4 using the approach of GOBAI-O₂ (GOBAI-O₂-ESM4), and subsampled grid cells from GFDL-ESM4 when and where real observations are available. Global weighted means (μ) of grid-cell-level values are shown, along with differences (Δ) of the fully resolved GFDL-ESM4 means versus GOBAI-O₂-ESM4 and versus the subsampled GFDL-ESM4 grid cells.

Figure 3 Two-dimensional histograms showing grid-cell-level (a) climatological seasonal amplitudes in monthly mean [O₂] (weighted means according to the size of each pressure layer) from 0 to 200 dbar and (d) trends in annual mean [O₂] from 200 to 1000 dbar between GFDL-ESM4 and GOBAI-O₂-ESM4. Pearson's correlation coefficients between GFDL-ESM4 and GOBAI-O₂-ESM4 for (b) monthly mean [O₂] from 0 to 200 dbar, showing coherence between the surface seasonal cycles, and (e) annual mean [O₂] from 200 to 1000 dbar, showing coherence between the subsurface trends. Absolute difference between GFDL-ESM4 and GOBAI-O₂-ESM4 for (c) climatological seasonal amplitudes in monthly mean [O₂] from 0 to 200 dbar and (f) trends in annual mean [O₂] from 200 to 1000 dbar. In panels (e) and (f), stippling indicates grid cells in which the GFDL-ESM4 trend is not significantly different from zero.

3.1.1 Test with withheld observational data

Estimates of [O₂] using the ENS_Data-Eval algorithms tracked closely with [O₂]_obs and showed no strong systematic biases with [O₂]_est or depth (Fig. 2a, b), although variability in Δ[O₂] was greatest from just below the surface to about 500 dbar. Mean offsets were between −1.8 and 3.2 µmol kg⁻¹ for the seven regions, with a global average of −0.2 µmol kg⁻¹; RMSDs were between 7.1 and 11.0 µmol kg⁻¹ for the seven regions, with a global average of 8.8 µmol kg⁻¹ (Table 2). The slightly negative global average offset suggests somewhat higher estimated than measured [O₂] values, and some of the lowest RMSDs from the ENS_Data-Eval algorithms were found in the Southern Ocean regions (Table 2, Fig. 2c), likely because these regions have a significant number of available training data (Fig. 1). However, this evaluation exercise is influenced by the incomplete subset of data (20 %) used to test the ENS_Data-Eval algorithms. A cross-fold validation (e.g., repeating this exercise with five separate 20 % chunks of data withheld from algorithm training) was prohibitively computationally expensive. Therefore, the associated Δ[O₂] and RMSD values alone are not as instructive as a comparison to the Δ[O₂] and RMSD values obtained from the ESPER-Mixed model (Table 2).

Estimates of [O₂] using ESPER-Mixed (Fig. A7) showed average offsets between −3.9 and −1.4 µmol kg⁻¹ for the seven regions (with a global average of −2.6 µmol kg⁻¹) and RMSDs between 10.0 and 21.6 µmol kg⁻¹ for the seven basins (with a global average of 13.1 µmol kg⁻¹) (Table 2). Again, the negative global average offset suggests higher estimated than measured [O₂] values. Compared with ESPER-Mixed (Carter et al., 2021), the ENS_Data-Eval algorithms performed better, both in terms of Δ[O₂] and RMSD in each individual region and overall. This result is likely a reflection of the fact that the ENS_Data-Eval algorithms were trained with more varied data than the ESPER-Mixed model (Argo and GLODAP compared with just GLODAP) and that the withheld data for which estimates were made also comprised more varied data (both Argo and GLODAP as well). Importantly, when estimates were made for just the GLODAP dataset, the ENS_Data-Eval algorithms still performed better than ESPER-Mixed (Table B4), suggesting that the seasonally resolved float data supply important information to the relationships established during algorithm training.

3.1.2 Test with GFDL-ESM4 output

As introduced in Sect. 2.3, we refer to the four-dimensional field of [O₂]_est values calculated by applying ENS_ESM4-Eval algorithms to GFDL-ESM4 output as GOBAI-O₂-ESM4. [O₂]_est values from GOBAI-O₂-ESM4 tracked closely with [O₂]_mod and showed no significant systematic biases with [O₂] or depth (Fig. 2d, e). Similar to the data-based test, variability in Δ[O₂] was greatest from just below the surface to about 500 m. Average offsets were between −1.3 and 1.0 µmol kg⁻¹ for the seven regions (with a global average of −0.2 µmol kg⁻¹) and RMSDs were between 3.3 and 9.6 µmol kg⁻¹ for the seven basins (with a global average of 6.7 µmol kg⁻¹) (Table 2). The near-zero global average offset suggests that [O₂]_est values from GOBAI-O₂-ESM4 matched well with values from GFDL-ESM4 output. The lowest RMSDs were found in the Southern Ocean and Arctic regions (Table 2, Fig. 2f), again due to the high density of training data in these regions.

In addition to direct comparisons of [O₂] values, GOBAI-O₂-ESM4 effectively captured local decadal-scale and seasonal variability in [O₂] in the GFDL-ESM4 model environment (Figs. 3, A8, A9, A10; Table 3). The average Pearson's correlation coefficient between gridded monthly mean [O₂] integrated from 0 to 200 dbar from GFDL-ESM4 output and GOBAI-O₂-ESM4 was 0.92 ± 0.17 (Fig. 3b), and the seasonal amplitudes differed in magnitude (GFDL-ESM4 minus GOBAI-O₂-ESM4) by 1.8 ± 4.0 µmol kg⁻¹ (Fig. 3c). The average Pearson's correlation coefficient between gridded annual mean [O₂] integrated from 200 to 1000 dbar from GFDL-ESM4 output and GOBAI-O₂-ESM4 was 0.66 ± 0.37 (Fig. 3e), and the trends differed in magnitude (GFDL-ESM4 minus GOBAI-O₂-ESM4) by −0.3 ± 2.1 µmol kg⁻¹ per decade (Fig. 3f).

When considered on the global scale, mean values, seasonal-cycle amplitudes, long-term trends, and interannual variabilities in [O₂] matched well between GFDL-ESM4 output and GOBAI-O₂-ESM4 (Table 3). In almost every case, agreement was far better than it was when simply considering GFDL-ESM4 grid cells for which observations are available over this time period, with no spatiotemporal interpolation. For example, the trend in monthly mean [O₂] integrated from 0 to 2000 dbar was −0.38 µmol kg⁻¹ per decade for GFDL-ESM4 output versus −0.18 µmol kg⁻¹ per decade for GOBAI-O₂-ESM4 (Fig. A11). On the other hand, grid cells where observations are available actually indicated an increase in monthly mean [O₂] integrated from 0 to 2000 dbar of 7.3 µmol kg⁻¹ per decade over this time period when no spatiotemporal interpolation is applied.

Whether the internal variability in GFDL-ESM4 is truly representative of the ocean or is biased in one or more dimensions, the success of GOBAI-O₂-ESM4 in this evaluation exercise demonstrates an ability for the ML algorithms employed here to capture that variability using the current distribution of available [O₂] observations as training data. This bodes well with respect to the ability of GOBAI-O₂, which is trained on actual observational data, to represent decadal-scale and seasonal variability in global ocean oxygen in the real world. However, the GFDL-ESM4 output has undergone substantial spatial and temporal averaging and has no observational uncertainties, and thus the assessed skill can be thought of as an upper limit of the reconstruction skill achievable with the currently available observations.

The results of the exercise with GFDL-ESM4 model output are critical for evaluating the uncertainty in gridded oxygen values in GOBAI-O₂ (Sect. 3.2.4). Further, the spatial distribution of Δ[O₂] (Fig. 2f) and the comparisons of reconstructed to modeled decadal trends and seasonal variability (Figs. 3b, c, e, and f and A8–A10) can help inform our observing efforts (e.g., future cruise planning and BGC Argo float deployments). For example, large Δ[O₂] values in the eastern tropical Pacific and eastern tropical Atlantic, coupled with some negative correlations in annual mean [O₂] and differences in annual trends and seasonal amplitudes, suggest that more observations will be required for GOBAI-O₂ (or likely any observation-based gap-filled [O₂] data product) to fully capture variability in that region.

3.2 GOBAI-O₂ product

3.2.1 Spatial oxygen distribution

The full GOBAI-O₂ product is available at https://doi.org/10.25921/z72m-yz67 (Sharp et al., 2022). Vertical–meridional sections of oxygen (Figs. 4, 5) show that surface oxygen concentrations are generally high, as these waters tend to be near equilibrium with the atmosphere. This is particularly true at high latitudes, where cold, dense waters have a high capacity for dissolved oxygen. Southern Ocean surface waters, however, are generally undersaturated with respect to oxygen (Fig. A12), consistent with observations from previous studies that suggest this undersaturation is the result of O₂-depleted thermocline water upwelling into the mixed layer (Chierici et al., 2004; Reuer et al., 2007; Jonsson et al., 2013) making the Southern Ocean on average an oxygen sink (Gruber et al., 2001; Bushinsky et al., 2017). This phenomenon can also be seen in the equatorial Pacific (Fig. A12). Undersaturation in high-latitude regions that are ice-covered during parts of the year can also be the result of limited air–sea gas exchange when sea ice is present.

Figure 4 Long-term mean [O₂] from GOBAI-O₂ at (a) 300 dbar and from the surface to 2000 dbar in the (b) Pacific, (c) Indian, and (d) Atlantic oceans. Dashed white lines in panel (a) show the locations of the sections in panels (b)–(d). White contour lines in panels (b)–(d) are potential isopycnals (kg m⁻³).

Figure 5 Long-term mean [O₂] from GOBAI-O₂ at (a) ${\mathit{\sigma }}_{\mathit{\theta }}=\mathrm{27}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ and from σ_θ=24 to 27.5 kg m⁻³ in the (b) Pacific, (c) Indian, and (d) Atlantic oceans. Dashed white lines in panel (a) show the locations of the sections in panels (b)–(d). White contour lines in panels (b)–(d) are constant isobars (dbar).

Isobaric maps, isopycnal maps, and vertical–meridional sections with pressure and density vertical coordinates (Figs. 4, 5) also reveal the [O₂] signatures of distinct subsurface water masses. In each basin, well-ventilated Subtropical Mode Water can be identified by relatively high [O₂] at midlatitudes on the 300 dbar surface (Fig. 4a) and along dips in isopycnals plotted against pressure and latitude (Fig. 4b, c, d) or along sloping isobars plotted against density and latitude (Fig. 5b, c, d) within the upper ∼500 dbar. Beneath mode waters in the southern portions of each basin, Antarctic Intermediate Water that originates in the Southern Ocean with a relatively high [O₂] signal is prevalent. Beneath mode waters in the northern portions of the Pacific and Indian basins, respectively, relatively old and oxygen-poor North Pacific Intermediate Water (NPIW) and Red Sea Overflow Water (RSOW) can be observed (Talley et al., 2011). Beneath mode waters in the northern portion of the Atlantic basin, intermediate waters are younger and more highly oxygenated. Near the Equator, subsurface oxygen minima are visible in each basin; this is a result of organic matter export from high production in the surface ocean that fuels strong subsurface respiration and relatively poor ventilation (old waters) in this region. Finally, the signatures of higher-oxygen deep or bottom waters can be observed near the bottom or at high latitudes in each vertical–meridional section.

Oxygen concentrations at 300 dbar (Fig. 4a) are highest in the North Atlantic and Southern oceans – where highly oxygenated, newly formed deep and intermediate waters are formed – and lowest in the northern and equatorial Pacific Ocean and the North Indian Ocean – where the oxygen content of subsurface waters has been greatly reduced by heterotrophic respiration over time. The same can be said for [O₂] on the 27.0 kg m⁻³ σ_θ surface (Fig. 5a). Oxygen concentrations are extremely low in the deep North Pacific Ocean (Figs. 4b, 5b) and North Indian Ocean (Figs. 4c, 5c) due to the accumulated effects of oxygen-depleting respiration over the long lifespans of those water masses (i.e., long time since gas exchange with the atmosphere).

3.2.2 Climatological seasonal oxygen cycles

Seasonal cycles in [O₂] reflect a balance among physical and biological processes (Wang et al., 2022). The climatological hemispheric mean [O₂] integrated over three pressure layers from GOBAI-O₂ (Fig. 6) reveals that the magnitude of the [O₂] seasonal cycle is greatest near the surface and decreases with depth. The amplitude of the [O₂] seasonal cycle in a near-surface layer (0–100 dbar) is about 10.7 µmol kg⁻¹ in the Northern Hemisphere and 8.9 µmol kg⁻¹ in the Southern Hemisphere. Maximum [O₂] in this pressure layer (April/May in the Northern Hemisphere and October/November in the Southern Hemisphere) lags about 2 months behind the temperature minimum, suggesting an interaction between a thermally driven increase in oxygen solubility and biologically driven oxygen production. Minimum [O₂] in the near-surface layer (October in the Northern Hemisphere and March/April in the Southern Hemisphere) is more coincident with the temperature maximum, indicating primary control by a thermally driven decrease in oxygen solubility. The amplitude of the [O₂] seasonal cycle is about 2.4 µmol kg⁻¹ in the Northern Hemisphere and 2.6 µmol kg⁻¹ in the Southern Hemisphere in the intermediate layer (100–600 dbar), whereas it is about 0.2 µmol kg⁻¹ in the Northern Hemisphere and 0.1 µmol kg⁻¹ in the Southern Hemisphere in the deep layer (600–2000 dbar). The timing of maximum [O₂] values is similar between the near-surface layer and the intermediate layer in both hemispheres, indicating the well-mixed nature of the ocean in winter and early spring when [O₂] is high. On the other hand, minimum [O₂] in the intermediate layer lags behind that in the near-surface layer in both hemispheres, possibly reflecting higher stratification in the upper ocean when temperatures are warmer and/or the remineralization of sinking organic matter after summer production. Further analysis of climatological [O₂] cycles from GOBAI-O₂ can provide insight into the physical and biological factors that control surface and subsurface oxygen on regional and global scales.

Figure 6 Climatological seasonal cycles of [O₂] anomalies (monthly [O₂] minus long-term mean [O₂]) integrated globally over three pressure layers: 0–100, 100–600, and 600–2000 dbar. The dotted black line shows the climatological temperature anomaly integrated globally over the 0–100 dbar layer. Shading indicates the standard deviation of the climatological seasonal cycle from 2004 to 2022. The dashed lines show climatological seasonal cycles of [O₂] anomalies from WOA18 over similar depth layers to GOBAI-O₂: 0–100, 100–600, and 600–1500 m.

3.2.3 Interannual oxygen trends and variability

Deoxygenation is evident in GOBAI-O₂ over the past 2 decades, coincident with ocean warming (Fig. 7, Table B5). The spatially weighted rate of deoxygenation in the upper 2 km globally (along with a 95 % confidence interval) is −1.19 ± 0.05 µmol kg⁻¹ per decade (−0.79 ± 0.04 % per decade). The rate of deoxygenation in GOBAI-O₂ varies over depth, with a near-surface-pressure-layer (0–100 dbar) trend in [O₂] of −1.10 ± 0.60 µmol kg⁻¹ per decade (−0.49 ±0.26 % per decade), an intermediate-layer (100–600 dbar) trend in [O₂] of −1.38 ± 0.42 µmol kg⁻¹ per decade (−0.86 ±0.26 % per decade), and a deep-layer (600–2000 dbar) trend in [O₂] of −1.12 ± 0.54 µmol kg⁻¹ per decade (−0.79 ±0.38 % per decade). Interannual variability is greatest in the near-surface layer: when the multiyear trends and seasonal cycles are removed, the standard deviation of annual global mean [O₂] anomalies is 0.52 µmol kg⁻¹ in the near-surface layer compared with 0.26 µmol kg⁻¹ in the intermediate layer and 0.12 µmol kg⁻¹ in the deep layer. Trends and uncertainties were determined by fitting linear least-squares models to spatially weighted monthly mean [O₂] and monthly oxygen inventories integrated over the specified pressure layers, with uncertainties in monthly values determined by comparing GOBAI-O₂-ESM4 to GFDL-ESM4; more information on this is provided in Appendix E.

Figure 7 Annual mean (a) [O₂] anomalies from GOBAI-O₂, (b) temperature anomalies from RG09, and (c) [O₂]_sat. anomalies calculated from RG09 temperature and salinity fields, each integrated globally over three pressure layers: 0–100, 100–600, and 600–2000 dbar. In panel (a), shading represents uncertainty determined as the average difference of mean [O₂] from GOBAI-O₂-ESM4 versus GFDL-ESM4 in each layer. Hovmöller diagrams showing annual mean (d) [O₂] anomalies from GOBAI-O₂, (e) temperature anomalies from RG09, and (f) [O₂]_sat. anomalies calculated from RG09 temperature and salinity fields, each versus pressure in decibars and time from 2004 to 2022. Anomalies in each parameter are calculated as annual mean values minus the long-term mean either (a–c) integrated over a pressure layer or (d–f) on a given pressure level.

Ocean warming has a direct effect on oxygen concentrations by lowering the solubility of O₂ in ocean water (Garcia and Gordon, 1992). Solubility changes explain about 56 % of deoxygenation in the near-surface ocean layer (0–100 dbar), 20 % in the intermediate ocean layer (100–600 dbar), and 15 % in the deep ocean layer (600–2000 dbar) (Fig. 7c, f). The remaining deoxygenation must then be caused by the indirect consequences of ocean warming (such as increased ocean stratification and, hence, decreased subsurface ventilation) or other processes, including changes in oxygen utilization and ocean ventilation variability (Oschlies et al., 2018), the magnitudes of which this analysis does not attempt to deconvolve. The RG09 temperature and salinity fields are constructed such that they relax toward the climatological means during periods of low data density. For this reason, toward the beginning of the time series when fewer observations are available, temperature is biased somewhat high (Fig. 7b, e) and, therefore, O₂ solubility is biased somewhat low (Fig. 7c, f). This artifact may influence GOBAI-O₂ (Fig. 7a, d), as it was constructed using the RG09 temperature and salinity fields; however, its influence is partially mitigated because temporal information included in the training and application of the GOBAI-O₂ algorithms allows for the trend inherent in the underlying oxygen data to be retained.

GOBAI-O₂ trends can be viewed in the context of other recent analyses that have explored long-term changes in ocean oxygen. From the surface to 1000 dbar, the GOBAI-O₂ trend of −0.82 ± 0.11 % per decade from 2004 to 2022 is larger than that assessed by Bindoff et al. (2019) of −0.48 ± 0.35 % per decade from 1970 to 2010 (surface to 1000 m), which takes into account estimates from Helm et al. (2011) (−0.44 ± 0.14 % per decade), Schmidtko et al. (2017) (−0.34 ± 0.35 % per decade), and Ito et al. (2017) (−0.68 ± 0.33 % per decade). In the surface layer (0–100 dbar), the GOBAI-O₂ trend of −0.49 ± 0.26 % per decade can be compared to the Bindoff et al. (2019) assessment of −0.28 ± 0.24 % per decade; in the intermediate layer (100–600 dbar), the GOBAI-O₂ trend of −0.86 ± 0.26 % per decade can be compared to the Bindoff et al. (2019) assessment of −0.52 ± 0.36 % per decade. Considering that these comparisons represent different periods of time such that one should not expect perfect agreement, we find the results encouraging. The somewhat more negative GOBAI-O₂ trends compared with previous estimates suggest a possible acceleration of ocean deoxygenation over the last decade or so, which would be consistent with expectations (Kwiatkowski et al., 2020). Further, agreement between GOBAI-O₂ and other observation-based studies provides additional support for the notion that current ESMs, which exhibit weaker deoxygenation trends (see Sect. 3.1.2), may not fully capture the sensitivities of physical and biological processes leading to deoxygenation (Oschlies et al., 2017, 2018; Stramma and Schmidtko, 2021). This comparison not only places the GOBAI-O₂ trends in a longer-term context but suggests that the enhanced observations and analysis result in a reduced trend uncertainty despite the comparatively shorter 19-year record (±0.04 % per decade) versus the longer but more sparse 40-year records assessed by Bindoff et al. (2019) (±0.14 to ±0.35 % per decade).

The trends presented here represent both natural and potentially anthropogenic variability over the interval between 2004 and 2022 as well as uncertainties in the algorithm predictions (see Sect. 3.2.4). As such, these trends should not be interpreted to be driven exclusively by ocean warming and other associated impacts of anthropogenic climate change; the period of time examined is relatively short, and the domain is not inclusive of the entire global ocean. Accordingly, decadal-scale variability in ocean ventilation, interior circulation, and biological oxygen utilization may exert significant influence over these trends. This is especially true of the regional trends. Finally, a sensitivity test indicated that including versus withholding temporal predictors did not significantly impact the global [O₂] trend. However, a shift from the relative dominance of shipboard observations during the early portion of the GOBAI-O₂ time span to the relative dominance of float observations during the later portion cannot be ignored as a potential contributor to a deoxygenation trend, especially considering the potential for systematic biases in the float [O₂] dataset (see Sect. 2.5). This kind of shift in measurement platforms has precedent for producing spurious trends in oceanographic observational datasets (Rykaczewski and Dunne, 2011).

3.2.4 Uncertainty

GOBAI-O₂ uncertainty fields, which were estimated as described in Sect. 2.5, can be used to assess confidence in multiyear trends and seasonal cycles of [O₂], both on a global and regional scale. Time-averaged uncertainty fields at 150 dbar (Fig. 8) suggest that the largest driver of geographic variability in uncertainty is the algorithm uncertainty. Averaged globally over space and time, u([O₂])_meas. was equal to 4.5 µmol kg⁻¹ (5.6 µmol kg⁻¹ on the 150 dbar level), u([O₂])_grid. was equal to 3.1 µmol kg⁻¹ (5.3 µmol kg⁻¹ on the 150 dbar level), and u([O₂])_alg. was equal to 3.8 µmol kg⁻¹ (6.3 µmol kg⁻¹ on the 150 dbar level). Combined, u([O₂])_tot. (Eq. 1) was equal to 7.6 µmol kg⁻¹ (11.2 µmol kg⁻¹ on the 150 dbar level), which can be compared to the global average RMSD of 8.8 µmol kg⁻¹ determined independently by withholding data from algorithm training (Table 2; Fig. 2a, b, c).

Figure 8 Long-term means of the uncertainty contributors to GOBAI-O₂ at 150 dbar, including (a) measurement uncertainty, (b) gridding uncertainty, (c) algorithm uncertainty, and (d) total uncertainty.

Figure 9 The difference between the climatological mean [O₂] from WOA18 and the long-term mean [O₂] from GOBAI-O₂ ($\mathrm{\Delta }\left[{\mathrm{O}}_{\mathrm{2}}\right]=[{\mathrm{O}}_{\mathrm{2}}{]}_{\text{WOA}}-[{\mathrm{O}}_{\mathrm{2}}{]}_{\text{GOBAI}}$) at (a) 300 m and from the surface to 1500 dbar in the (b) Pacific, (c) Indian, and (d) Atlantic oceans.

Measurement uncertainty provides an estimate of confidence in an [O₂] value assigned to a water sample by direct measurement, gridding uncertainty provides an estimate of confidence that the [O₂] value provided for a four-dimensional grid cell might represent [O₂] at any point in time and space within that grid cell, and algorithm uncertainty provides an estimate of confidence that the predicted [O₂] value for a given grid cell is appropriate as the average value for that grid cell. Algorithm uncertainty, in particular, depends upon the distribution of data available to train the ML algorithms and the ability of the trained algorithms to represent underlying variability in the system. On the isobar shown in Fig. 8. (150 dbar), the underlying variability is relatively high in [O₂] minimum zones (e.g., near the Equator and on the eastern boundaries of ocean basins), hence the elevated algorithm (and total) uncertainties in those regions. Here, algorithm uncertainty was assessed via the exercise with synthetic data from GFDL-ESM4 (see Sects. 2.3 and 3.1.2).

Algorithm uncertainty should generally decrease as the spatiotemporal coverage of available training data increases. Regionally, algorithm uncertainty depends upon the degree to which the underlying variability in the system is captured by the available training observations and the ability of the ML algorithms to reconstruct that variability from concurrent measurements of other seawater properties. A comparison of the Δ[O₂] map in Fig. 2f or the algorithm uncertainty map in Fig. 8c to the data distribution map in Fig. 1 or Fig. A1 suggests that sparse sampling is primarily to blame for limitations related to algorithm uncertainty. Detailed analysis of GFDL-ESM4 water mass characteristics in the California Current system has also revealed that high uncertainties occur where water masses with similar physical characteristics but different oxygen signatures mix, underscoring that the measurement of additional biogeochemical parameters can supplement the [O₂] estimates presented here, which are based on physical properties and spatiotemporal information. Overall, the significant influence of algorithm uncertainty is consistent with uncertainty analyses conducted for gap-filling methods applied to other ocean biogeochemical variables (e.g., Landschützer et al., 2014; Gregor and Gruber, 2021). For this reason, continued expansion of oxygen observations in undersampled regions will be critical to reduce uncertainty in our gap filling, and ultimately our understanding, of global subsurface oxygen distributions and variability. Similarly, the significant influence of measurement uncertainty underscores the importance of continued development of oxygen sensor calibration (Bittig et al., 2018a) and data quality control (Maurer et al., 2021) from the evolving BGC Argo fleet.

Global mean depth profiles of uncertainty contributors (Fig. A14) emphasize the general attenuation of uncertainty away from the surface, with subsurface maxima of algorithm uncertainty at 200 dbar and total uncertainty at 100 dbar. The algorithm uncertainty maximum corresponds to pressures at which vertical gradients in [O₂] are relatively high (see Fig. 4). Here, small variations in the depths of density surfaces can influence [O₂] at a given pressure level; this variability is challenging to capture, even with potential density as a predictor variable in the ML models (see Table 1). The total uncertainty maximum represents this vertical gradient effect balanced against relatively high measurement uncertainty closer to the surface, associated with higher [O₂].

3.2.5 Comparison to other gridded products

The long-term mean field of [O₂] from GOBAI-O₂ was compared to the corresponding mean field of [O₂] from the WOA18 monthly climatology (Fig. 9) and the climatological field of [O₂] from the GLODAPv2.2016 mapped product (Fig. 10). On average, the GOBAI-O₂ oxygen concentration is 1.4 µmol kg⁻¹ lower than GLODAP and 9.6 µmol kg⁻¹ lower than WOA18. This can be partly explained by the fact that GOBAI-O₂ is centered on the year 2012, whereas observations in GLODAPv2.2016 are centered around 2002 (Lauvset et al., 2016), WOA18 takes into account [O₂] observations dating back to 1965, and global deoxygenation has occurred in recent decades (Bindoff et al., 2019). Spatially, the largest differences occur within and especially near the boundaries of oxygen minimum zones (eastern tropical Pacific, eastern Atlantic coastal, and northern Indian zones), along σ≈27.5 kg m⁻³ in the Southern Ocean, and along σ≈26.75 kg m⁻³ in the North Pacific. It is difficult to determine whether these differences are functions of data availability (ship data for WOA18 and GLODAP versus ship and float data for GOBAI-O₂), representative time period, or mapping method (objective interpolation for WOA18 and GLODAP versus machine learning algorithms for GOBAI-O₂). A future intercomparison exercise between mapping methods using an identical starting dataset could be helpful in diagnosing these differences among gridded products.

Figure 10 The difference between the 2002-centered mean [O₂] from the GLODAPv2.2016 mapped product and the long-term mean [O₂] from GOBAI-O₂ ($\mathrm{\Delta }\left[{\mathrm{O}}_{\mathrm{2}}\right]=[{\mathrm{O}}_{\mathrm{2}}{]}_{\text{GLODAP}}-[{\mathrm{O}}_{\mathrm{2}}{]}_{\text{GOBAI}}$) at (a) 300 m and from the surface to 1500 dbar in the (b) Pacific, (c) Indian, and (d) Atlantic oceans.

3.2.6 Comparison to synoptic in situ measurements

GOBAI-O₂ was compared to direct observations from repeat hydrography cruises, including meridional transects across the Atlantic (A16 in 2013 and A20 in 2021), Pacific (P16 in 2005), and Indian (I08 and I09 in 2016) oceans, as well as a zonal transect across the Pacific Ocean (P02 in 2012). This exercise assessed how well monthly [O₂] estimates from GOBAI-O₂ were able to represent high-quality [O₂] measurements at distinct points in time and space. Due to fundamental differences between gridded estimates and point observations, we do not expect every matchup to be perfect. However, we would hope to see general coherence in mean values across large-scale ocean sections and to see a pattern of differences that make sense given our a priori expectations.

For the cruise datasets examined, GOBAI-O₂ estimates matched fairly well with discrete measurements in the mixed layer and below ∼1000 dbar (Fig. 11). At intermediate depths, however, large differences occasionally occur. These large differences tended to cluster around areas with strong vertical gradients in [O₂] (thin contours in Fig. 11 represent increments of 50 µmol in [O₂]). Comparison of Fig. 11 to Fig. A15 provides confidence in our uncertainty evaluation: larger differences between discrete measurements and GOBAI-O₂ occur where u([O₂])_tot. is large. Median biases, mean biases, and RMSDs between direct observations and GOBAI-O₂ are given in Table B6.

Figure 11 Section plots displaying comparisons between discrete observations of [O₂] from repeat hydrography cruises and [O₂] extracted from corresponding grid cells in GOBAI-O₂. The thick lines in each panel represent the mixed-layer depth calculated as the depth at which the potential density anomaly increased to 0.03 kg m⁻³ greater than the potential density anomaly at 10 dbar. The thin lines are contours representing increments of 50 µmol kg⁻¹ in [O₂].

4 Data availability

GOBAI-O₂ is available as a NetCDF file at https://doi.org/10.25921/z72m-yz67 (Sharp et al., 2022); additional information and animations can be found at https://www.pmel.noaa.gov/gobai/ (last access: 28 September 2023). GLODAPv2.2022 is available at https://doi.org/10.25921/1f4w-0t92 (Lauvset et al., 2022b); the GLODAP database is updated annually and is made available at https://www.glodap.info (last access: 28 September 2023). GFDL-ESM4 model output is available at https://doi.org/10.22033/ESGF/CMIP6.1407 (Krasting et al., 2018). Data from the 2018 World Ocean Atlas are available at https://www.ncei.noaa.gov/archive/accession/NCEI-WOA18 (Garcia et al., 2019; last access: 01 Oct. 2021). The OneArgo-Mat toolbox used to download Argo float data is available at https://doi.org/10.5281/zenodo.6588041 (Frenzel et al., 2022); the toolbox acquires data from two global data assembly centers: Coriolis and US GODAE (https://doi.org/10.17882/42182, Argo, 2000). The Roemmich and Gilson (2009) Argo-based temperature and salinity product is available at https://sio-argo.ucsd.edu/RG_Climatology.html (last access: 12 January 2023).

5 Conclusions

GOBAI-O₂ is a major step toward the fulfillment of the primary goal set out by Gruber et al. (2010): “to determine, on a global-scale, seasonal to decadal time-scale variations in dissolved oxygen concentrations throughout the upper ocean”. Quantifying these variations is important for documenting ocean deoxygenation, determining global net primary productivity and carbon export, and facilitating studies of the oceanic uptake of anthropogenic CO₂. In addition, insights into ocean biogeochemical dynamics, when observations are unavailable, often come from ocean models, and GOBAI-O₂ can bring value to modeling studies by providing fields of [O₂] to be used for boundary conditions and model initialization. GOBAI-O₂ can also be useful as a dynamic reference check for new, sensor-based [O₂] measurements that would otherwise be compared to a static monthly climatology like WOA18. Still, users should carefully consider the associated spatial uncertainty fields, especially when conducting regional analyses. Spatial and temporal errors and discontinuities may be significant when GOBAI-O₂ is analyzed over small areas, but they are mitigated when looking at broader scales.

The uncertainty analysis conducted here confirms that GOBAI-O₂ remains limited, largely by sparse sampling and inadequate representation of [O₂] across strong gradients. The most consequential actions to improve GOBAI-O₂ fields over the next decade will be the continued deployment of Argo floats with oxygen optodes, emphasizing the importance of bolstering the biogeochemical Argo array and expanding the international OneArgo network into high latitudes, the deep ocean, and marginal seas (Roemmich et al., 2019, 2021; Schofield et al., 2022); continued work toward ensuring reliable measurements from those optodes; and the continued collection of discrete dissolved oxygen observations, primarily through the international GO-SHIP program, both for use in [O₂] mapping and for calibration/validation of the Argo oxygen data.

Besides these actions, additional steps can be taken to improve GOBAI-O₂ fields. For one, more predictor variables and ML algorithms can be tested. Different processes dominate [O₂] variability in different regions (Keeling et al., 2010; Oschlies et al., 2018; Garcia-Soto et al., 2021), and certain predictor variables will be better suited for capturing these processes. Also, ML algorithms adapt to data sparseness and modes of variability in different ways (Ritter et al., 2017; Gregor et al., 2019), so estimates in a given region that are worse using one algorithm may be better using another. Therefore, regionally tuned predictors and more diverse ensembles of ML algorithms should lead to increased confidence in estimates of ocean interior [O₂]. Another action that could result in improved fidelity of GOBAI-O₂ fields is the use of predictor variable fields with higher spatial and temporal resolution across sharp biogeochemical gradients. Ocean profiles of temperature and salinity tend to be relatively smooth, so a depth resolution on the order of tens of meters in the upper ocean increasing to hundreds of meters at depth is sufficient for gridded products. Biogeochemical parameters like oxygen, on the other hand, tend to be characterized by profiles with sharp gradients and with distinct minima and maxima in the water column (Sarmiento and Gruber, 2006). These minima and maxima can occur very near the surface or hundreds of meters below it. For this reason, comparisons of GOBAI-O₂ to direct measurements of [O₂] can be uniquely problematic in the ∼100–1000 dbar range when sharp gradients are present (Fig. 11). A complicating factor is the lack of response-time corrections applied to float sensor data (Sect. 2.5), which contributes to uncertainty in observation-based [O₂] products like GOBAI-O₂. Biogeochemical gradients over horizontal space and time can also be sharp, especially in highly dynamic coastal zones and in the surface ocean where the residence time of oxygen is often less than a month (Luz and Barkan, 2000). Recent work from Lyman and Johnson (2023) uses Argo observations coupled with machine learning to provide well-resolved ($\mathrm{7}\mathrm{d}\times \mathrm{1}/\mathrm{4}{}^{\circ }$ grid) ocean heat content maps, and continued development toward maps of temperature and salinity could be helpful for overcoming the issue of resolving sharp biogeochemical gradients. Alternatively, [O₂] estimates could be made using temperature and salinity observations at their original resolution and then mapped onto four-dimensional grids that are uniquely suited in their spatial resolution for biogeochemical parameters. A necessary consideration of the latter option would be computing resources: applying complex ML algorithms to temperature and salinity measurements from Argo floats at their original resolution may prove to be impractical. Finally, observations from additional platforms could be incorporated into approaches like this one to map [O₂] in the global ocean. Ocean gliders and moored profilers have long been equipped with oxygen optodes. These platforms collect data at unique spatiotemporal scales and could add predictive information for [O₂] that is not provided by Argo float observations nor discrete shipboard measurements. To facilitate the incorporation of new data streams into the development of gridded data products, accessible databases should be created and maintained (Testor et al., 2019; Grégoire et al., 2021).

The method used to develop GOBAI-O₂ can be applied in a similar way to other ocean chemical parameters. In addition to dissolved oxygen, the BGC Argo program has deployed floats with sensors for measuring dissolved nitrate, pH, chlorophyll-a, particle backscatter, and downwelling irradiance. Machine learning methods have been used to develop four-dimensional fields of optical properties, i.e., chlorophyll-a and particle backscatter (Sauzède et al., 2015, 2016), and continued refinement of those fields is ongoing (Sauzède et al., 2021). Chemical properties, i.e., nitrate and pH, that exhibit distributions more similar to [O₂] are good candidates for adoption into the GOBAI mapping approach. Together with property estimation algorithms for total alkalinity (Bittig et al., 2018b; Carter et al., 2021), a mapped ocean interior pH product could be used to resolve the entire ocean carbonate system in four dimensions in near-real time.

Ultimately, global changes in the amount of dissolved oxygen in ocean waters will have profound effects on the metabolism of marine organisms (Pörtner and Farrell, 2008; Sampaio et al., 2021) and the cycling of biogeochemically important elements (Gruber, 2004; Berman-Frank et al., 2008). Whereas ocean models agree that the ocean's oxygen inventory has been declining and will continue to decline, disagreement remains as to regional patterns of this deoxygenation. Direct observations are critical for the confirmation or contradiction of model trends. With this work, we have turned to autonomous and discrete observations, with the assistance of machine learning algorithms, to bridge the model–observational gap. We produce and analyze a multiyear gridded product of ocean dissolved oxygen called GOBAI-O₂, independently confirming a phenomenon that has been demonstrated previously: the ocean is losing dissolved oxygen at a rapid rate (0.79 ± 0.04 % per decade in the upper 2 km according to GOBAI-O₂). In addition, we provide this valuable observation-based product for community use. GOBAI-O₂ can be turned to as a reference for [O₂] observations and model boundary conditions, compared with new and existing observational and model-based reconstructions of ocean deoxygenation, and used for critical analyses of seasonal to decadal and regional to global oxygen variability.

Appendix A

Figure A1 The number of profiles (either ship-based or Argo-float-based) from the combined dataset used to train machine learning algorithms to produce GOBAI-O₂ that are contained within each $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ box in the global ocean.

Figure A2 Annual mean in situ (a) temperature and (b) salinity from RG09 (2004–2022) at 20 dbar.

Figure A3 Annual mean in situ (a) temperature, (b) salinity, and (c) dissolved oxygen concentration from GFDL-ESM4 (2004–2021) at 20 m.

Figure A4 A schematic for the random forest regressions (RFRs) and feed-forward neural networks (FNNs). A random subset of the predictors is used for each tree in the RFR, and a randomly chosen predictor is used for each node split. The two hidden layers (H₁ and H₂) in each of the three FNNs have 10 and 20, 15 and 15, and 20 and 10 nodes. Each machine learning algorithm is trained with input data and [O₂] observations and then used to predict [O₂] from new input data.

Figure A5 The spatial distribution of profile data used to (a) train and (b) test the RFR_Data-Eval and FNN_Data-Eval algorithms.

Figure A6 A comparison between the number of observations binned within a four-dimensional grid cell and the standard deviation in [O₂] among those observations. The horizontal black line shows the mean standard deviation (5.21 µmol kg⁻¹).

Figure A7 For withheld Argo and GLODAP data, two-dimensional histograms showing offsets between measured and ESPER-Mixed-calculated oxygen ($\mathrm{\Delta }\left[{\mathrm{O}}_{\mathrm{2}}\right]=[{\mathrm{O}}_{\mathrm{2}}{]}_{\text{meas}}-[{\mathrm{O}}_{\mathrm{2}}{]}_{\text{ESPER}}$) as a function of (a) [O₂]_ESPER and (b) pressure in the water column. Offsets are binned into cells that are 2.5 µmol kg⁻¹ tall in terms of Δ[O₂] and 5 µmol kg⁻¹ wide in terms of (a) [O₂]_ESPER or (b) equivalent in width to the interpolated pressure levels of the data. (c) Absolute Δ[O₂] values averaged over depth and time for 1^∘ × 1^∘ (latitude × longitude) grid cells in the global ocean.

Figure A8 Integrated from 0 to 200 dbar: (a, b) long-term mean [O₂], (d, e) seasonal [O₂] amplitudes, (g, h) trends in [O₂], and (j, k) interannual variability in [O₂] for (a, d, g, j) GFDL-ESM4 and (b, e, h, k) GOBAI-O₂-ESM4, along with (c, f, i, l) the difference between the two.

Figure A9 Integrated from 200 to 1000 dbar: (a, b) long-term mean [O₂], (d, e) seasonal [O₂] amplitudes, (g, h) trends in [O₂], and (j, k) interannual variability in [O₂] for (a, d, g, j) GFDL-ESM4 and (b, e, h, k) GOBAI-O₂-ESM4, along with (c, f, i, l) the difference between the two.

Figure A10 Integrated from 0 to 2000 dbar: (a, b) long-term mean [O₂], (d, e) seasonal [O₂] amplitudes, (g, h) trends in [O₂], and (j, k) interannual variability in [O₂] for (a, d, g, j) GFDL-ESM4 and (b, e, h, k) GOBAI-O₂-ESM4, along with (c, f, i, l) the difference between the two.

Figure A11 Monthly area-weighted mean [O₂] integrated globally from 0 to 2000 dbar from GFDL-ESM4 (blue) and GOBAI-O₂-ESM4 (orange).

Figure A12 Long-term mean oxygen saturation percentage on the uppermost pressure level in GOBAI-O₂.

Figure A13 Monthly mean de-seasonalized (a) [O₂] anomalies from GOBAI-O₂, (b) temperature anomalies from RG09, and (c) [O₂]_sat. anomalies calculated from RG09 temperature and salinity fields, each integrated globally over three pressure layers: 0–100, 100–600, and 600–2000 dbar. In panel (a), shading represents uncertainty determined as the average difference of mean [O₂] from GOBAI-O₂-ESM4 versus GFDL-ESM4 in each layer. Hovmöller diagrams showing monthly mean de-seasonalized (d) [O₂] anomalies from GOBAI-O₂, (e) temperature anomalies from RG09, and (f) [O₂]_sat. anomalies calculated from RG09 temperature and salinity fields, each over depth in decibars from 2004 to 2022. Anomalies in each parameter are calculated as monthly mean values with a seasonal cycle removed and minus the long-term mean either (a–c) integrated over a pressure layer or (d–f) on a given pressure level.

Figure A14 Global mean depth profiles of uncertainty contributors to GOBAI-O₂, including (a) measurement uncertainty, (b) gridding uncertainty, (c) algorithm uncertainty, and (d) total uncertainty. The shaded region represents variability in space and is calculated as the standard deviation on each depth level of the mean uncertainties over time for each grid cell.

Figure A15 Section plots displaying total uncertainty estimates from GOBAI-O₂ that correspond to discrete measurements of [O₂] from repeat hydrography cruises, to be compared to Δ[O₂] values in Fig. 11.

Appendix B

Table B1 Boundaries for the seven large ocean regions used to fit machine learning algorithms.

* Two sets of boundaries are given for the Pacific to accommodate crossing the International Date Line.

Table B2 Error statistics (mean Δ[O₂] and RMSD) for tests using RFR and FNN algorithms trained on a subset of Argo and GLODAP data and tested with a separate subset of withheld data. Also shown are error statistics corresponding to the ensemble average (ENS) of the estimates from both algorithms.

Table B3 Error statistics (mean Δ[O₂] and RMSD) for tests using RFR and FNN algorithms trained on a subset of output from GFDL-ESM4 (corresponding to locations of available Argo and GLODAP data) and tested using a separate subset of withheld output from GFDL-ESM4. Also shown are error statistics corresponding to the ensemble average (ENS) of the estimates from both algorithms.

Table B4 Error statistics (mean Δ[O₂] and RMSD) for tests using RFR and FNN algorithms trained on a subset of Argo and GLODAP data and tested with all available GLODAP data. Also shown are error statistics corresponding to the ensemble average (ENS) of the estimates from both algorithms and corresponding to the ESPER-Mixed model (Carter et al., 2021).

Table B5 Estimated decadal trends and uncertainties in [O₂] (µmol kg⁻¹ per decade) and oxygen inventory (% per decade) in different pressure layers of GOBAI-O₂. Uncertainties are determined according to the procedure in Appendix E, both using the autocorrelation of residuals to the linear least-squares model (Autocov.) and by incorporating estimated uncertainty in global mean GOBAI-O₂ fields (u(ESM4)). The value used to represent uncertainty in each trend (larger value) is in bold.

Table B6 Summary error statistics between direct observations from repeat hydrography cruises and GOBAI-O₂ and WOA18.

Appendix C

As stated above, the original and vertically interpolated observational datasets from the BGC Argo and GLODAP databases that are used to develop GOBAI-O₂ are archived online: https://doi.org/10.5281/zenodo.7747237 (Sharp, 2023b). The algorithms trained on vertically interpolated observational data that were applied to predictor variables to produce GOBAI-O₂ are also archived online: https://doi.org/10.5281/zenodo.7747308 (Sharp, 2023a).

Appendix D

A negative median bias (−1.18 µmol kg⁻¹) in float [O₂] measurements compared with co-located ship [O₂] measurements (below 300 dbar, to avoid the impact of high-frequency variability near the surface) was adjusted by fitting the differences (Δ[O₂]) to a linear least-squares model as a function of float [O₂] and then adding that [O₂]-dependent adjustment back on to the float [O₂] measurements. The Δ[O₂] values as a function of float [O₂] before (Fig. D1a) and after (Fig. D1b) this adjustment are shown. This resulted in a reduced median Δ[O₂] of 0.33 µmol kg⁻¹.

Figure D1 Unadjusted (a) and adjusted (b) matchups between BGC Argo [O₂] measurements (y axis) and GLODAP [O₂] measurements (x axis). The adjustment procedure does not mitigate the scatter between the matchups, but it does reduce the median error.

Appendix E

Trends and associated uncertainties in GOBAI-O₂ were determined via the four-step procedure outlined below.

1. Spatially weighted monthly mean [O₂] values for the entire GOBAI-O₂ domain or within specified pressure layers were calculated from gridded [O₂], using relative grid-cell volumes as weights.

2. A linear least-squares model with a trend, intercept, and annual (12-month) and semiannual (6-month) harmonics was fit to monthly mean [O₂] values. The monthly trend from the least-squares model was multiplied by 120 to obtain a decadal trend of weighted mean [O₂].

3. Uncertainty in the decadal trends was assessed in two different ways, and the largest of the two uncertainty estimates taken for each analyzed pressure layer, indicating that either (a) uncertainty in the linear least-squares model or (b) uncertainty in the GOBAI-O₂ fields was driving uncertainty in the trend. In all cases, the uncertainty estimate from the second method was chosen, indicating that decadal trend uncertainty was controlled primarily by uncertainty in the GOBAI-O₂ fields. The two different methods used to assess uncertainty are outlined below.

   • a.

     Uncertainty was assessed using the autocovariance of residuals from the linear least-squares model via the following procedure:

     • i.

       the standard error on the trend was calculated from the covariance matrix of the linear least-squares model;

     • ii.

       the autocovariance of the residuals from the least-squares model was examined to compute the e-folding timescale, and the effective degrees of freedom were obtained by dividing the number of monthly mean [O₂] values by the e-folding timescale and subtracting the number of least-squares model parameters; and

     • iii.

       the standard error on the trend was scaled by the effective degrees of freedom, multiplied by 2 to obtain a 95 % confidence interval, and multiplied by 120 to obtain an uncertainty on the decadal trend of weighted mean [O₂].

   • b.

     Uncertainty was assessed by incorporating estimated uncertainty in global mean GOBAI-O₂ fields via the following procedure:

     • i.

       uncertainties in monthly mean [O₂] values were determined as the standard deviations of monthly differences between GOBAI-O₂-ESM4 and GFDL-ESM4 (Sect. 3.1.2);

     • ii.

       these uncertainties were used to compute a weight matrix for the linear least-squares fit, and the effective degrees of freedom were obtained as previously described; and

     • iii.

       the standard error on the trend was scaled by the effective degrees of freedom, multiplied by 2 to obtain a 95 % confidence interval, and multiplied by 120 to obtain an uncertainty on the decadal trend of weighted mean [O₂].

4. The process was repeated for oxygen inventories for the entire GOBAI-O₂ domain or within each specified pressure layer; inventories were determined from gridded [O₂], volumes of each grid cell, and densities of each grid cell.