The Atmospheric Potential Oxygen forward Model Intercomparison Project (APO-MIP1): evaluating simulated atmospheric transport of air-sea gas exchange tracers and APO flux products

Atmospheric Potential Oxygen (APO, defined as O₂ + 1.1 × CO₂) is primarily a tracer of ocean biogeochemistry and fossil fuel burning. APO exhibits strong seasonal variability at mid-to-high latitudes, driven mainly by seasonal air-sea O₂ exchange. We present results from the first version of the Atmospheric Potential Oxygen forward Model Intercomparison Project (APO-MIP1), which forward transports three air-sea APO flux products in eight atmospheric transport models or model variants, aiming to evaluate atmospheric transport and flux representations by comparing simulations against surface station, airborne, and shipboard observations of APO. We find significant spread and bias in APO simulations at eastern Pacific surface stations, indicating inconsistencies in representing vertical and coastal atmospheric mixing. A framework using airborne APO observations demonstrates that most atmospheric transport models (ATMs) participating in APO-MIP1 overestimate tracer diffusive mixing across moist isentropes (i.e., diabatic mixing) in mid-latitudes. This framework also enables us to isolate ATM-related biases in simulated APO distributions using independent mixing constraints derived from moist static energy budgets from reanalysis, thereby allowing us to assess large-scale features in air-sea APO flux products. Furthermore, shipboard observations show that ATMs are unable to reproduce seasonal APO gradients over Drake Passage and near Palmer Station, Antarctica, which could arise from uncertainties in APO fluxes or model transport. The transport simulations and flux products from APO-MIP1 provide valuable resources for developing new APO flux inversions and evaluating ocean biogeochemical processes.

1 Introduction

Atmospheric potential oxygen (APO), defined as the weighted sum of O₂ and CO₂ concentration (APO ≈ O₂ + 1.1 CO₂), is an important tracer of fossil fuel burning and ocean biogeochemical processes (Stephens et al., 1998). APO is intended to be unaffected by terrestrial photosynthesis and respiration due to the cancellation of O₂ and CO₂ exchange at an approximate O₂ : C ratio of −1.1 (Severinghaus, 1995). APO exhibits a large seasonal cycle driven mainly by air-sea O₂ exchange due to upper ocean biological activities, deep water ventilation, and thermally induced O₂ solubility changes. Seasonal APO variability is also slightly affected by the air-sea exchange of CO₂ and N₂ (Manning and Keeling, 2006). APO is decreasing in the atmosphere due to fossil fuel combustion, which acts as an O₂ sink and CO₂ source with a more negative O₂ : CO₂ ratio (global mean $\sim -\mathrm{1.4}$) compared to the assumed −1.1 ratio from terrestrial processes. Although fossil fuel combustion contributes to an annual interhemispheric gradient that has lower APO in the Northern Hemisphere, it has only a minor effect on the seasonal cycle globally (Keeling and Manning, 2014).

APO measurements provide critical constraints on seasonal air-sea O₂ fluxes, which have been used to estimate air-sea gas exchange rates and ocean net community production (NCP), and to benchmark marine NCP in Earth system models (Naegler et al., 2007; Nevison et al., 2018, 2012, 2015, 2016). APO has been used for improved partitioning of ocean and land carbon sinks (Friedlingstein et al., 2025; Manning and Keeling, 2006), to constrain ocean heat uptake and meridional heat transport (Resplandy et al., 2016, 2019), and to quantify fossil fuel emissions (Pickers et al., 2022; Rödenbeck et al., 2023). APO measurements are available at surface stations (e.g., Adcock et al., 2023; Battle et al., 2006; Goto et al., 2017; Keeling and Manning, 2014; Manning and Keeling, 2006; Nguyen et al., 2022; Tohjima et al., 2019), on ship transects (e.g., Ishidoya et al., 2016; Pickers et al., 2017; Stephens et al., 2003; Thompson et al., 2007; Tohjima et al., 2012, 2015, 2024), and from aircraft (e.g., Bent, 2014; Ishidoya et al., 2012; Jin et al., 2023; Langenfelds, 2002; Morgan et al., 2021; Stephens et al., 2018, 2021g).

Global-scale air-sea APO fluxes have been estimated from APO measurements and an ATM within a Bayesian inversion framework (Rödenbeck et al., 2008). ATMs are also used to forward transport APO fluxes simulated from ocean biogeochemistry models (Carroll et al., 2020; Yeager et al., 2022) and surface ocean dissolved oxygen (DO) measurements (Garcia and Keeling, 2001; Najjar and Keeling, 2000) to compare with atmospheric observations, providing a basis for model and flux product evaluation (Jin et al., 2023; Keeling et al., 1998; Stephens et al., 1998). However, using atmospheric data to evaluate flux products and to derive fluxes through inversion is fundamentally limited by biases in ATMs, particularly in their representation of vertical transport and diabatic mixing (Jin et al., 2024; Naegler et al., 2007; Nevison et al., 2008; Schuh et al., 2019; Schuh and Jacobson, 2023; Stephens et al., 2007). The systematic uncertainties in transport modeling limit inversions of APO, CO₂, and other greenhouse gases, underscoring the need for independent transport bias assessments to advance global carbon budget constraints.

To address uncertainty in ATMs for studying large-scale tracer atmospheric transport and the corresponding surface fluxes, several community model intercomparison (TransCom) projects have been established for various tracers including CO₂ (Baker et al., 2006; Gurney et al., 2003, 2004; Law et al., 2008; Patra et al., 2008), N₂O (Thompson et al., 2014), SF₆ (Denning et al., 1999), SF₆ and CH₄ jointly (Patra et al., 2011), as well as an age of air tracer (Krol et al., 2018). Blaine (2005) coordinated a TransCom O₂ experiment to compare model simulations of the O₂ seasonal cycle across the Scripps O₂ network. While this experiment provided valuable initial insights into ATM performance in simulating atmospheric O₂ from ocean fluxes, substantial advances in ATMs and more data collected also from aircraft and ships since then motivate an updated intercomparison study with more extensive model-data comparisons and analyses. More recently, CO₂ inversion intercomparisons have been coordinated through the OCO-2 MIP (Crowell et al., 2019; Peiro et al., 2022; Byrne et al., 2023) and the Global Carbon Project (e.g., Friedlingstein et al., 2025). These experiments reveal substantial spread in forward tracer (e.g., CO₂) atmospheric distribution and inverted surface fluxes, driven by different ATMs and inversion setups. The spread in forward transport simulations stems from multiple factors, including the choice of wind fields from various reanalysis products or online simulation, regridding fine resolution meteorological data to coarse model grids, the advection scheme that governs large-scale mixing, and parameterized sub-grid processes, such as boundary layer mixing and deep convection. Despite the complexity of different transport pathways, long-lived tracers (e.g., CO₂ and O₂) at mid-latitudes tend to show tracer distributions that are aligned with moist potential temperature (θ_e) surfaces. This is because θ_e surfaces are preferential surfaces for mixing, leading to rapid along-θ_e mixing and slow cross-θ_e mixing (Bailey et al., 2019; Jin et al., 2021; Miyazaki et al., 2008; Parazoo et al., 2011).

It is a critical challenge to accurately quantify the rate-limiting cross-θ_e mixing time-scales, which are largely driven by diabatic processes including moist convection and radiative cooling. Here, we define “diabatic mixing rates” as diffusivities that are inversely related to cross-θ_e mixing time-scales. These mixing rates are important for determining the large-scale tracer distribution in ATMs. Jin et al. (2024) established a framework to calculate cross-θ_e mixing rates from ATMs and moist static energy (MSE) budgets from reanalysis based on a mass-indexed isentropic coordinate called $M_{{\mathit{\theta }}_e}$ (Jin et al., 2021). This framework also allows cross-θ_e tracer gradients from airborne observations to provide independent constraints on diabatic mixing. Jin et al. (2024) tested four ATMs used in CO₂ inversions, showing that these models tend to have too fast mixing in the mid-latitudes of the Southern Hemisphere in the austral summer. The too fast mixing is also confirmed by the fact that models simulate smaller CO₂ gradients compared to airborne observations, which is an independent constraint on the mixing rate. The mixing rate constraint and CO₂ gradient constraint also have implications for biases in the inverse model estimates, indicating a too large summer-time Southern Ocean (SO) CO₂ sink. This framework provides a system for independently evaluating transport simulations and flux estimates.

Previous TransCom experiments focused primarily on tracers that only have significant sources and sinks over the land, and large seasonal flux cycles tied to the northern terrestrial biosphere. In contrast, APO is a tracer of surface ocean exchange with the largest seasonal variability observed over mid-to-high latitude oceans in both hemispheres. APO offers a distinct perspective for studying atmospheric mixing within and above the marine boundary layer, the long-range tracer transport into and out of the remote Southern Hemisphere, and the ability for inverting tracer flux over the SO from atmospheric measurements.

Here we use output from the APO-MIP1 (Stephens et al., 2025), which generated a suite of forward ATM simulations of APO and its components (air-sea O₂, CO₂, and N₂ flux, and fossil fuel CO₂ emission and O₂ uptake) from different source fields. This effort was initially motivated by a need to support the calibration of hemispheric-scale seasonal air-sea APO flux estimates from spatially and temporally sparse observations from airborne campaigns (e.g., Jin et al., 2023), stations, and ships. Here we focus on the other goals of APO-MIP1 which were to use atmospheric APO observations to characterize errors in ATMs and APO flux products.

In Sect. 2, we describe APO measurements from surface stations, aircraft, and ships, and the experimental design of APO-MIP1. In Sect. 3.1, we evaluate simulations against observations, revealing large model spread and errors at eastern Pacific surface stations due to mixing uncertainties, while airborne column-average data show smaller cross-ATMs variability and errors. In Sect. 3.2, we analyze diabatic mixing rates, demonstrating that ATMs generally overestimate mid-latitude mixing in both hemispheres, allowing us to separate transport and flux-related biases. In Sect. 3.3, we examine simulations of shipboard data around Drake Passage and the Antarctic Peninsula, revealing that current ATMs and flux products underestimate meridional gradients in APO seasonal amplitude from 53–65° S. The models also fail to capture the APO contrast between Palmer Station flask samples and nearby in-situ ship data due to limitations in representing local topographic flows with coarse-resolution ATMs. In Sect. 3.4, we discuss the broader implications of our analysis for developing methods to identify processes that introduce transport biases and for improving atmospheric transport modeling.

Figure 1 Geographic distribution of APO observations used in this study: (a) Scripps O₂ Program surface stations (red diamonds) with station codes and inlet elevation in meters above sea level; (b) HIPPO (1 to 5) airborne campaign horizontal flight tracks covering the Pacific Ocean; (c) ORCAS aircraft measurements concentrated in the Drake passage; (d) ATom (1 to 4) airborne campaign horizontal flight tracks covering the Pacific and Atlantic Oceans; and (e) Ship-based measurements from the RV Laurence M. Gould operating in the Drake passage.

2 Materials and methods

2.1 Definition of APO

APO in the unit of per meg (see Keeling et al., 1998) is calculated from atmospheric observations of relative changes in the O₂ $/$ N₂ ratio (per meg) and CO₂ mole fraction (ppm) according to Stephens et al. (1998) as

$$\begin{array}{}\text{(1)}& \mathrm{APO}\approx \mathit{\delta }\left({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\right)+{\displaystyle \frac{\mathrm{1.1}}{X_{{\mathrm{O}}_{\mathrm{2}}}}}\left({\mathrm{CO}}_{\mathrm{2}}-\mathrm{350}\right),\end{array}$$

with

$$\begin{array}{}\text{(2)}& \mathit{\delta }({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}})=\left({\displaystyle \frac{{\left(\frac{{\mathrm{O}}_{\mathrm{2}}}{{\mathrm{N}}_{\mathrm{2}}}\right)}_{\mathrm{sample}}}{{\left(\frac{{\mathrm{O}}_{\mathrm{2}}}{{\mathrm{N}}_{\mathrm{2}}}\right)}_{\mathrm{reference}}}}-\mathrm{1}\right)\cdot {\mathrm{10}}^{\mathrm{6}}.\end{array}$$

The factor 1.1 represents the approximate exchange ratio of O₂ to CO₂ in terrestrial biospheric processes (Severinghaus, 1995). We note that this ratio generally varies from 1.01 to 1.14 in aboveground carbon pools across different temporal and spatial scales (Gallagher et al., 2017; Hockaday et al., 2009; Keeling, 1988; Worrall et al., 2013). This ratio also exhibits diurnal change and varies between respiration and photosynthesis in biosphere-atmosphere O₂ and CO₂ exchanges (Faassen et al., 2023, 2024). With our focus on seasonal variations, we use 1.1 as representative of the O₂ to CO₂ exchange ratio during seasonal growth and decay of terrestrial biota. A sensitivity test in Jin et al. (2023) showed that varying this ratio by ±0.05 only leads to ±5.1 % changes in hemispheric average APO. The impact on APO seasonal cycle amplitude (SCA) is ±1.44 % and ±0.41 % in the Northern and Southern Hemisphere, respectively. $X_{{\mathrm{O}}_{\mathrm{2}}}$ (0.2094) is the reference dry-air mole fraction of O₂ used in the definition of the O₂ scale of the Scripps O₂ Program (Keeling et al., 2020). δ(O₂ $/$ N₂) is expressed in units of per meg, while CO₂ is converted from ppm units to per meg units by subtracting a reference value of 350 ppm and then dividing by $X_{{\mathrm{O}}_{\mathrm{2}}}$. APO observations are typically expressed in per meg units, but they can be converted to ppm equivalent units by multiplying by $X_{{\mathrm{O}}_{\mathrm{2}}}$.

2.2 Atmospheric measurements

The APO-MIP1 (Stephens et al., 2025) required model output sampled to match a collection of surface station, airborne, and shipboard observations, and also accepted optional output at additional locations, at higher time resolution, and for full 3-D fields, as shown in Tables S1 and S2 in the Supplement. Here we evaluate model APO simulations using observation data collected at 10 surface stations, on 10 airborne campaigns from three projects, and one repeated shipboard transect from 50 cruises. We show sampling locations, and horizontal flight and ship tracks in Fig. 1. We use surface station APO measurements (2009 to 2018) from 10 sampling sites mainly in the Pacific from the Scripps O₂ Program surface flask network (Keeling and Manning, 2014; Manning and Keeling, 2006). The airborne measurements (Stephens et al., 2018) were made on the NSF NCAR GV aircraft during the HIAPER Pole-to-Pole Observation project from 2009 to 2011 (HIPPO, Wofsy, 2011) and the O₂ $/$ N₂ Ratio and CO₂ Airborne Southern Ocean Study in 2016 (ORCAS, Stephens et al., 2018), and from the NASA DC-8 aircraft during the Atmospheric Tomography Mission from 2016–2018 (ATom, Thompson et al., 2022). Shipboard measurements were made on transects crossing the Drake Passage by the NSF ARSV Laurence M. Gould from 2012–2017 (Stephens, 2025). Details of surface station, airborne, and shipboard APO measurements are provided in Appendix A.

As the primary focus of this study is the APO seasonal cycle and its latitudinal distribution, we remove interannual trends from the observational data. For surface station and airborne measurements, we remove the long-term trend by subtracting a deseasonalized cubic spline fit (smoothing parameter of 0.8) derived from the global mean APO time series using Scripps O₂ Program data following Hamme and Keeling (2008). For the ship data, we apply a similar detrending procedure but use only South Pole Observatory (SPO) data to derive the long-term trend.

2.3 Components of APO in the atmosphere and prescribed surface fluxes

APO exhibits seasonal variations primarily driven by air-sea exchange ($F_{\mathrm{APO}}^{\mathrm{ocn}}$), which comprises three components: air-sea exchange of O₂ ($F_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ocn}}$), CO₂ ($F_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ocn}}$), and N₂ ($F_{{\mathrm{N}}_{\mathrm{2}}}^{\mathrm{ocn}}$). Additionally, APO is influenced by fossil fuel emission of CO₂ ($F_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ff}}$) and consumption of O₂ ($F_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ff}}$), which together combine to form a sink for APO due to fossil fuel burning ($F_{\mathrm{APO}}^{\mathrm{ff}}$). Fluxes are defined as positive to the atmosphere.

In this study, we primarily simulate APO by performing forward transport of these individual flux components in ATMs, except one inverse model flux product that provides net $F_{\mathrm{APO}}^{\mathrm{ocn}}$ directly. We combined these components to calculate the net atmospheric APO anomalies in units of per meg as

$$\begin{array}{}\text{(3)}& \mathit{\delta }\mathrm{APO}=\mathit{\delta }{\mathrm{APO}}^{\mathrm{ocn}}+\mathit{\delta }{\mathrm{APO}}^{\mathrm{ff}},\end{array}$$

with

$$\begin{array}{}\text{(4)}& \mathit{\delta }{\mathrm{APO}}^{\mathrm{ocn}}={\displaystyle \frac{\mathrm{1}}{X_{{\mathrm{O}}_{\mathrm{2}}}}}\cdot \mathrm{\Delta }{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ocn}}-{\displaystyle \frac{\mathrm{1}}{X_{{\mathrm{N}}_{\mathrm{2}}}}}\cdot \mathrm{\Delta }{{\mathrm{N}}_{\mathrm{2}}}^{\mathrm{ocn}}+{\displaystyle \frac{\mathrm{1.1}}{X_{{\mathrm{O}}_{\mathrm{2}}}}}\cdot \mathrm{\Delta }{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ocn}},\end{array}$$

and

$$\begin{array}{}\text{(5)}& \mathit{\delta }{\mathrm{APO}}^{\mathrm{ff}}={\displaystyle \frac{\mathrm{1}}{X_{{\mathrm{O}}_{\mathrm{2}}}}}\cdot \mathrm{\Delta }{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ff}}+{\displaystyle \frac{\mathrm{1.1}}{X_{{\mathrm{O}}_{\mathrm{2}}}}}\cdot \mathrm{\Delta }{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ff}}.\end{array}$$

where $\mathrm{\Delta }{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ocn}}$, $\mathrm{\Delta }{{\mathrm{N}}_{\mathrm{2}}}^{\mathrm{ocn}}$, $\mathrm{\Delta }{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ocn}}$, $\mathrm{\Delta }{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ff}}$, and $\mathrm{\Delta }{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ff}}$ represents the atmospheric fields in units of deviations in ppm of each flux component ($F_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ocn}}$, $F_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ocn}}$, $F_{{\mathrm{N}}_{\mathrm{2}}}^{\mathrm{ocn}}$, $F_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ff}}$, and $F_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ff}}$) that is forward transport in the ATMs (Stephens et al., 1998). The δ sign denotes tracers in units of per meg.

We utilize three distinct ocean APO flux products: (1) the Jena product, which directly provides $F_{\mathrm{APO}}^{\mathrm{ocn}}$ from an atmospheric APO inversion framework that assimilates surface station measurements (Rödenbeck et al., 2008); (2) the CESM product, an Earth System Model simulation with prognostic ocean biogeochemistry (Yeager et al., 2022; Long et al., 2021a) that generates separate flux components ($F_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ocn}}$ and $F_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ocn}}$); and (3) the DISS product, which provides separate observation-based flux components incorporates surface ocean dissolved oxygen measurements (Garcia and Keeling, 2001; Resplandy et al., 2016) and pCO₂ data (Jersild et al., 2017; Landschützer et al., 2016). $F_{{\mathrm{N}}_{\mathrm{2}}}^{\mathrm{ocn}}$ for CESM and DISS is estimated by scaling ocean heat fluxes from CESM and ERA-5, respectively, using the relationship of Keeling et al. (1993). For fossil fuel contributions, we employ the OCO2MIP product for CO₂ emissions (Basu and Nassar, 2021) and the GridFED database for coupled O₂ and CO₂ fluxes from fossil fuel combustion (Jones et al., 2021). Details of each product are provided in Appendix B. All flux fields were linearly interpolated from their original temporal and spatial resolution to 1° longitude × 1° latitude with daily temporal resolution from 1986 to 2020. When flux data were unavailable in the earlier portion of this time period (Jena and OCO2MIP), we set the corresponding fluxes to zero. Participating modelers were requested to simulate at least from 2009 to 2018, following three years of spin up from 2006 to 2008, and optionally longer (Table 1). In addition to Jena, which is simulated directly, we construct the two ΔAPO^ocn products using Eq. (4) and two ΔAPO^ff products using Eq. (5), as described in Appendix B. Figure 2 illustrates the seasonal and latitudinal flux patterns of these three ocean APO flux products and the fossil fuel APO flux from GridFed, which serves as our primary fossil fuel flux dataset in this study.

Figure 2 Comparison of APO flux patterns from the three air-sea flux products (Jena, CESM, and DISS) and fossil fuel emissions (GridFed), averaged from 2009 to 2018. (a) Hovmöller diagrams showing the spatiotemporal distribution of APO fluxes (mol m⁻² per month) as a function of latitude and month. (b) Seasonal cycles of zonally integrated fluxes for three latitude bands: Northern Extratropics (≥20° N, orange), Tropics (20° S–20° N, lavender), and Southern Extratropics (<20° S, green). (c) Latitudinal profiles of flux seasonal cycle amplitude (SCA, red) and annual mean flux (blue). For the annual mean profiles (blue lines in panel c), only the latitudinal gradients should be interpreted, as the global means may contain biases in the ocean flux products, which are not the focus of this paper.

Table 1 Participating ATMs and model parameters.

2.4 Atmospheric tracer transport models

We simulate each component of APO in the atmosphere using the flux fields described in Sect. 2.3, and eight ATMs (see Table 1). All tracer atmospheric fields are modeled as tracer deviations against an arbitrary background with concentrations in ppm dry air mole fraction (as for CO₂). These tracer mole fractions are later converted to deviations in units of per meg after subtracting the model-specific arbitrary reference according to Eq. (4). We describe key model parameters and setups below.

2.4.1 CAM-SD

The Community Atmosphere Model (CAM) version 6.0 is the atmospheric component of CESM2 (Danabasoglu et al., 2020). The version used here is run online with specified dynamics (SD), wherein the model is constrained with MERRA-2 reanalysis, and uncoupled from the other climate system components. Temperature and horizontal winds (u and v) are nudged to MERRA-2, 8 times per day, with a normalized strength coefficient of 0.25. Shallow convection is parameterized following the Cloud-Layers Unified by Binormals framework (CLUBB, Golaz et al., 2002), and deep convection is parameterized following Zhang and McFarlane (1995). CAM has not been used for tracer inversions, but has been evaluated extensively for its dynamical properties (e.g., Bailey et al., 2019; Kay et al., 2012).

2.4.2 CAMS_LMDZ

CAMS_LMDZ refers here to the offline transport model from the Atmospheric General Circulation Model of Laboratoire de Météorologie Dynamique, called LMDz. LMDz is the atmospheric component of the Earth System Model of Institut Pierre-Simon-Laplace (IPSL). It is also used to drive the offline model CAMS_LMDz, in which case its horizontal winds are nudged to those of the ERA5 reanalysis wind fields (Hersbach et al., 2020). From the computer code of LMDz, CAMS-LMDz only keeps the transport subroutines for advection (Hourdin and Armengaud, 1999), deep convection (Emanuel, 1991), thermals (Rio and Hourdin, 2008), and boundary-layer turbulence (Hourdin et al., 2006). All other processes are replaced by an archive of relevant meteorological variables (air mass fluxes, exchange coefficients, temperature, etc.) built with the full LMDz model at the target spatial resolution, thereby allowing relatively small computing time and resources for the offline model. LMDz ensures the physical consistency of the archive of meteorological variables. The meteorological variables are stored as 3-hourly averages. CAMS_LMDZ has been regularly participating in OCO-2 MIP (Byrne et al., 2023) and TransCom intercomparison studies.

2.4.3 CTE_TM5

TM5 is a tracer transport model used for simulating atmospheric trace gas chemistry and transport (Krol et al., 2005). We refer to it as CTE_TM5 because the model was run with the CarbonTracker-Europe (CTE) shell, but this does not alter the TM5 physics and chemistry. TM5 advection is computed using the slopes advection scheme (Russell and Lerner, 1981) and in this work it is driven by ERA-5 reanalysis wind fields, making it an offline model. The convection is computed from the convective entrainment and detrainment fluxes from the ERA-5 reanalysis. Free tropospheric diffusion is computed using the formulation by Louis (1979). Diffusion in the boundary layer is computed using the parametrization by Holtslag and Boville (1993), where the diurnal variability in the boundary layer height is computed using Vogelezang and Holtslag (1996). TM5 is widely used in inversions and regularly participates in MIPs, for different tracers at different model resolutions and driven with different wind reanalysis products (for example, Byrne et al., 2023; Friedlingstein et al., 2025; Gaubert et al., 2019; Krol et al., 2018).

2.4.4 TM3

TM3 (Heimann and Körner, 2003) is an offline atmospheric tracer transport model, in the present runs driven by meteorological fields from the NCEP reanalysis (Kalnay et al., 1996). It was run here on a spatial resolution of 5° longitude, about 3.8° latitude, and 19 vertical layers. The advection uses the slopes scheme (Russell and Lerner, 1981), which is the same as in TM5. Boundary layer mixing is parameterized according to Louis (1979). Vertical mixing due to sub-gridscale cumulus clouds is calculated using the mass flux scheme of Tiedtke (1989). TM3 is the ATM used in Jena APO inversion (Rödenbeck et al., 2008), which is one of the flux products used in this study.

2.4.5 MIROC4-ACTM

MIROC4-ACTM is a new generation Model for Interdisciplinary Research on Climate (MIROC, version 4.0; Watanabe et al., 2008) atmospheric general circulation model (AGCM)-based chemistry-transport model (ACTM; Patra et al., 2018). This AGCM is evolved from the Center for Climate System Research, University of Tokyo (CCSR)/National Institute for Environmental Studies (NIES)/Frontier Research Center for Global Change, JAMSTEC (FRCGC) AGCM version 5.7b (Numaguti et al., 1997). The MIROC4 AGCM propagates only explicitly resolved gravity waves into the stratosphere through the implementation of a hybrid vertical coordinate system compared to its predecessor AGCM5.7b. The MIROC4 AGCM online-simulated horizontal winds and temperature are nudged to the Japanese 55-year Reanalysis (JRA-55) at 6-hourly time intervals (Kobayashi et al., 2015). MIROC4-ACTM produces “age-of-air” up to about 5 years in the tropical upper stratosphere (∼1 hPa) and about 6 years in the polar middle stratosphere (∼10 hPa), in agreement with observational estimates. The convective transport and inter-hemispheric transport of tracers in the model are validated using ²²²Radon and sulphur hexafluoride (SF₆), respectively (Patra et al., 2018).

2.4.6 NICAM-TM_gl5 and NICAM-TM_gl6

NICAM-TM is an atmospheric transport model based on the Nonhydrostatic Icosahedral Atmospheric Model (NICAM) (Niwa et al., 2011; Satoh et al., 2014). In this study, we used the offline mode of NICAM-TM, which uses air mass fluxes, vertical diffusion coefficients and other meteorological variables; those data are calculated in advance by an online calculation of NICAM, in which horizontal winds are nudged toward the JRA-55 data. In NICAM, the air mass fluxes are calculated consistently with the continuity equation while conserving tracer masses, which do not require any numerical mass fixing (Niwa et al., 2011). For APO-MIP1, two horizontal resolutions were used: “glevel-5” (gl5) and “glevel-6” (gl6), whose mean grid intervals are 223 and 112 km, respectively. The number of the vertical model layers is 40 and the top of the model domain is at approximately 45 km. The vertical diffusion coefficients are calculated with the MYNN (Mellor and Yamada, 1974; Nakanishi and Niino, 2004) Level 2 scheme (Noda et al., 2010). The cumulus parameterization scheme used in NICAM-TM is Chikira and Sugiyama (2010). Model performance for atmospheric constituent transport can be found in Niwa et al. (2011, 2012).

2.4.7 NIES

NIES-TM-FLEXPART is a coupled transport model combining Eulerian (NIES-TM) and Lagrangian (FLEXPART) models. It is a transport modeling component of the variational flux inverse modeling system NIES-TM-FLEXPART-Variational (NTFVAR, Maksyutov et al., 2021). The NIES Transport Model (NIES-TM) is an offline model, originally developed in the 1990s (Maksyutov et al., 2008). In this study, the NIES-TM v.21 is used, which improves SF₆ transport and tropopause height over the former v.08.1 (Belikov et al., 2013), as evaluated in Krol et al. (2018), due to (a) using ERA5 hourly wind data, including vertical wind on model coordinates, on 137 model levels and a 0.625° grid for preparation of the 4-hourly average mass fluxes on 42 hybrid-pressure levels, (b) transporting first-order moments (Russell and Lerner, 1981; Van Leer, 1977) for advection, (c) applying penetrative convection rate and turbulent diffusivity supplied by the ERA5 reanalysis (Hersbach et al., 2020). The version v.21 is the same as used in the OCO-2 MIP (Byrne et al., 2023). NIES-TM is coupled with the Lagrangian model FLEXPART (Stohl et al., 2005) to provide refinement to the near field transport during the last 3 d prior to the observation event as presented by Belikov et al. (2016). FLEXPART model v.8.0 is driven by 6-hourly JRA-55 winds, interpolated to 40 hybrid pressure levels and 1.25°×1.25° resolution. The surface flux footprints are produced by FLEXPART at 1°×1° resolution and daily time step.

2.5 Outputs from transport models

For each ATM, we required simulations for all species sampled to match with the observation locations and times in a subset of the full ObsPack CO₂ files GLOBALVIEWplus v7.0 ObsPack (Schuldt et al., 2021), excluding the model spin-up period. This subset corresponds to existing APO observations that are analyzed in this study from Scripps O₂ Program surface stations, NSF NCAR airborne observations, and NSF NCAR and AIST/JMA shipboard programs. The full list of these records is in Table S1. We note that, while the HIPPO, ORCAS, ATom, and Gould ObsPack files contain CO₂ observations from different instruments, their 10-sec sampling times align with the NSF NCAR APO measurements, except during calibration periods for either instrument.

We also received optional output, which includes the full set of ObsPack files, 3-D atmospheric fields, meteorological variables, additional ship data, and output at additional fixed sites (Table S2). Further details are provided in the APO-MIP1 protocol available at Stephens et al. (2025). We obtained output matching the full set of ObsPack files from four ATMs, which will be useful for future network design. We obtained daily mean 3-D gridded concentration fields from six ATMs. These fields support the calculation of diabatic mixing rates, which we use to evaluate ATMs and the flux products, following the method of Jin et al. (2024). Details are in Sect. 3.2. We also received hourly (from two versions of NICAM) or 3-hourly (from NIES) output for an extensive list of sites with past or ongoing APO measurements, and co-located samples for ship sampling programs of NIES VOS, AIST RV Mirai, and UEA Cap San Lorenzo (Hamburg Süd) from three models. These data are not analyzed in this study, but are made available at Stephens et al. (2025).

3 Results and discussion

3.1 APO model-observation comparisons at surface stations and along aircraft flight tracks

3.1.1 APO seasonal and latitudinal variations at surface stations

We show observations and model simulations of APO seasonal cycles at 10 surface stations of the Scripps O₂ program network in Fig. 3. We present annual mean values, seasonal cycle amplitudes (SCA), and phase from both observations and model simulations at these surface stations in Fig. 4, with model-observation differences shown as colors. Observations show clear meridional gradients in APO annual means (Fig. 5a), with higher values in the Southern Hemisphere than Northern Hemisphere, and a southern tropical “bulge” evident at SMO and in the airborne data centered on 15° S (Battle et al., 2006; Gruber et al., 2001; Stephens et al., 1998). The APO SCA shows higher values in the high latitudes of both hemispheres, with larger amplitudes in the Southern Hemisphere compared to the Northern Hemisphere, yet reaches its maximum at the northern mid-latitude station Cold Bay (CBA) (Fig. 4b). The seasonal phase exhibits an approximately 6-month difference between hemispheres, while remaining relatively uniform within each hemisphere (Fig. 4c).

Figure 3 Comparison of simulated and observed APO seasonal cycles at 10 surface stations (Fig. 1a), organized from southern high-latitudes (left) to northern high-latitudes (right). In each panel, the black line represents observations, while colored lines show simulations from different transport models. Each row of panels corresponds to the three different flux products (Jena, CESM, and DISS). In each panel, the y-axis shows APO anomalies in per meg units, and the x-axis shows months from January to December. We note that, for LJO and CBA simulations using the Jena fluxes, a different y-axis range (three times larger) is used compared to the other panels. Observations and model simulations at each station are first detrended using a multiple-station weighted average trend. We calculate monthly mean seasonal APO from 2009 to 2018 for both observations and model simulations.

Figure 4 Evaluation of APO (a) annual mean relative to a multi-station global mean, (b) seasonal cycle amplitude, and (c) seasonal minimum day across surface stations using different flux-transport model combinations. For each panel, results are organized by flux products (JENA, CESM, DISS) in columns and transport models in rows, with observations on the top. The metrics are printed in black, with background colors indicating biases relative to observations. Positive bias is shown in red, and negative bias is shown in blue. Stations from left to right are organized by latitudes from south to north.

Figure 5 APO annual means (a, d), SCA (b, e), and seasonal minimum day (c, f) derived from airborne observations. In (a)–(c), we show latitude-pressure distributions, with data binned into 10° latitude by 100 mbar pressure boxes. In (d)–(f), we show 1000–400 mbar column-averaged (black) and 900 mbar interpolated (blue) values, and also surface station observations (2009 to 2018). Annual mean is derived from a two-harmonic fit with constant offset, where the global multi-station trend has been subtracted to detrend the airborne observations and center the values around zero globally. SCA is calculated as the peak-to-trough amplitude of the two-harmonic fit, and seasonal minimum day is calculated as the day of seasonal trough of the two-harmonic fit.

The higher annual mean APO in the Southern Hemisphere and the southern tropical “bulge” is a result of southward O₂ and CO₂ transport by the oceans, further amplified by net APO uptake in the Northern Hemisphere from fossil fuel burning (Keeling and Manning, 2014; Stephens et al., 1998). The larger APO SCA in mid- to high-latitudes reflects more pronounced seasonal flux cycles resulting from larger marine net primary production (NPP) and sea surface temperature changes in these regions. The thermal and biological effects on APO SCA are further enhanced at eastern Pacific coastal sites (e.g., LJO and CBA), where the shallow marine boundary layer traps high-APO air masses during summer. The 180 d phase difference between the two hemispheres is a result of different seasonal heating and cooling, as well as the biological cycle.

3.1.2 Biases in APO-MIP1 simulations at surface stations

APO-MIP1 simulations of APO annual means and seasonal cycles at surface stations broadly agree with observations (Figs. 3 and 4). Simulations driven by CESM fluxes show the best agreement with observed APO features. For annual mean spatial patterns (Fig. 4a), CESM- and DISS-driven simulations show comparable performance in representing the southern tropical “bulge” and north-south gradient in annual means, while significantly outperforming simulations using the Jena flux model in northern stations. The main limitation of simulations using CESM fluxes is an overestimation of annual mean APO values across Pacific sites in the Southern Hemisphere, and an underestimation at LJO. Simulations using DISS fluxes also underestimate the annual mean APO at LJO.

APO SCA is well represented in simulations driven by CESM flux, but the SCA at LJO is significantly underestimated in all ATMs except CAM-SD. The underestimation is caused by an overly weak summer-time APO peak (Fig. 3), which also leads to the small annual mean presented above. Simulations using DISS flux generally underestimate SCA, especially in the high latitudes. Simulations using Jena flux, however, generally overestimate the SCA in the mid- to high-latitudes. We find largest SCA biases and cross-ATMs spread at LJO and CBA when using the Jena flux. The biases and model spread are closely related to underrepresentation in ATMs, and will be discussed in the next section. We note that the model biases and spread observed at surface stations are smaller than those reported in the previous TransCom-O₂ experiment (Blaine, 2005), indicating improved atmospheric transport modeling.

Phase simulations using CESM flux are consistent with observations at most stations, except at two northern low-latitude stations, KUM and MLO, where we find too late seasonal minimum day by up to two weeks. Simulations using DISS flux show even larger biases, with earlier seasonal minimum days at all southern and northern low-latitude stations.

3.1.3 Impact of ATM mixing biases

We find APO-MIP1 simulations have large model spread and biases at two northern mid-latitudes stations, LJO and CBA (Fig. 3), especially simulations using Jena fluxes. We note that the interdependence of transport models and fluxes in inversions can be seen for the Jena flux product simulations at LJO (Figs. 3 and 4). As expected, we see good agreement with observations for the Jena flux product transported by the same model used in the Jena APO inversion (Jena_TM3). However, all other ATMs overestimate summertime APO, and consequently SCA, for the Jena flux product at LJO, CBA, and BRW. All other ATMs also simulate too negative wintertime APO at LJO. These biases suggest a stronger regional APO source in the Jena flux product that could have resulted from too rapid dilution of surface flux signals at LJO in both summer and winter.

Surface station simulations using CESM flux (Figs. 3 and 4) also reveal elevated model spread and observation deviations at LJO and CBA. At LJO, all ATMs underestimate summertime APO, and consequently SCA, implying too weak upwind outgassing fluxes. The relative magnitude of simulated summer-time peaks for CESM at LJO and CBA maintains a consistent pattern across different flux products, with CAM-SD consistently showing the highest values and Jena_TM3 the lowest, regardless of the flux product used, suggesting consistent biases in the ATMs.

This substantial cross-ATMs variability highlights the challenges in accurately representing complex atmospheric vertical transport processes in regions where strong temperature inversions and stratocumulus clouds significantly influence vertical mixing (Naegler et al., 2007; Nevison et al., 2008). The Jena flux product, derived from an inversion that assimilates these station data, relies on the TM3 tracer transport model (Rödenbeck et al., 2008). Previous studies indicate that TM3 consistently overestimates vertical mixing over the Eastern Pacific, leading to larger inverted seasonal fluxes to match station observations (Jin et al., 2023; Naegler et al., 2007). Our analysis suggests that in comparison to Jena_TM3, vertical mixing is weaker in the two versions of NICAM, CAM-SD, MIROC4-ACTM, and CTE_TM5, which show larger summer-time APO anomalies at LJO and CBA. This pattern is consistent across the three flux products considered.

The larger model spread at northern coastal sites (e.g., LJO, CBA, and BRW) also highlights the limitations of current coarse-resolution ATMs in representing horizontal coastal flows and land-sea breezes. At LJO, samples are collected only during steady west wind (from the ocean) conditions (Keeling et al., 1998). However, ATMs failed to capture the actual small-scale atmospheric conditions associated with on-shore winds during episodic storm systems, which leads to significant underestimation of oceanic influence (Keeling et al., 1998). APO, as a tracer of air-sea gas exchange, is particularly sensitive to the dilution effects in coarse-resolution models.

3.1.4 APO seasonal and latitudinal variations along flight tracks and biases in APO-MIP1

We present zonal averages of APO annual means, SCA, and seasonal minimum days derived from airborne data, grouped into 10° latitude and 100 mbar bands in Fig. 5a–c (full seasonal cycles in Fig. S1 in the Supplement). We further calculate these three metrics as column-average (black) and at 900 mbar (blue) in Fig. 5d–f, where we also compare them with surface station data (shown as red points). The airborne data show patterns similar to those seen at surface stations but provide detailed vertical structures. The vertical profiles consistently show larger SCA at low altitudes, indicating that the main drivers of SCA are near the surface, while annual means and seasonal phases remain uniform across altitudes. Airborne column averages show increasing SCA and decreasing annual means from low to high latitudes, with similar SCA and annual mean values north of 50° N (Fig. 5d and e), whereas station observations show peaks in the mid-latitudes (LJO and CBA) due to high-APO air masses being trapped below the summer marine boundary layer. This trapping effect is also evident in airborne data interpolated to 900 mbar.

Figure 6 Comparison of column-average (1000–400 mbar) APO features across latitude from aircraft observations and model simulations using three different flux products (Jena, CESM, and DISS). The figure is organized into three sets of panels showing (a) annual mean APO relative to a multi-station global mean, (b) SCA, and (c) seasonal minimum day. For each feature, we show latitudinal distributions of observations (black lines) and model simulations (colored lines). We note that the global mean value has been subtracted from the annual mean values (a) at each latitude to highlight spatial patterns. We show the column-average (400–1000 mbar) seasonal cycles of observed and simulated APO for each 10° latitude band in Fig. S1.

We also calculate APO annual means, SCA, and phases using aircraft simulations from APO-MIP1 (full seasonal cycles in Fig. S1) and compare simulated and observed column averages (1000–400 mbar average) in Fig. 6, with biases in column averages and vertical profiles shown in Figs. S2–S5, respectively. Airborne observation-model comparisons complement those using surface station data. We find similar model biases to those seen in surface data, for example, larger SCA at northern high latitudes with the Jena flux product and smaller SCA at high latitudes with the DISS flux product. The airborne data also reveal three key biases that are not resolved at surface stations. Observations suggest a consistent near-zero annual mean APO in the Southern Hemisphere (south of 30° S), with a spike between 40 and 50° S. However, all three flux products show gradually decreasing annual mean APO south of 30° S, with CESM and DISS flux products showing a smaller spike in magnitude between 40 and 50° S. Simulations using CESM and DISS flux products show larger annual mean values in the northern mid-latitudes (40–60° N). Additionally, simulations using the Jena flux product in the low northern latitudes show a seasonal minimum day similar to the Southern Hemisphere phase. This bias is caused by low-latitude flux features in the Jena inversion that largely replicate the Southern Hemisphere cycle, likely due to limited observational constraints in this region (Jin et al., 2023).

Our analysis demonstrates that global airborne measurements provide distinct advantages over station data for evaluating large-scale flux patterns due to the reduced sensitivity of column averages to boundary-layer ATM transport uncertainties. While surface stations show substantial cross-model spread in simulated APO (Figs. 3 and 4), column-averaged airborne simulations (Fig. 6) reveal remarkable consistency across ATMs when driven by the same flux product. This consistency suggests that column-averaged measurements effectively integrate over local transport features that often dominate surface observations. Here we establish CESM as the most realistic flux product among the three products. The better agreement between observations and CESM-driven simulations provides a more reliable baseline for isolating and quantifying transport-related discrepancies in individual ATMs.

3.2 Evaluation of diabatic mixing rates diagnosed from transport models

In this section, we evaluate the mixing timescale across mid-latitude moist isentropes of each ATM using the framework developed in Jin et al. (2024). This framework was applied to identify biases in four ATMs in the mid-latitude Southern Hemisphere using two independent constraints: (1) diagnosed diabatic mixing rates, and (2) cross-isentrope CO₂ gradients. Here we extend the framework to use APO gradients, to include two more reanalysis products, and the analysis in the Northern Hemisphere. We evaluate six of the eight ATMs participating in APO-MIP1 that provide 3-D atmospheric fields (CAM-SD, CTE_TM5, Jena_TM3, NICAM-TM_gl5, NICAM-TM_gl6, and MIROC4-ACTM), which are required to diagnose diabatic mixing rates. Diabatic mixing rates and APO gradients are diagnosed based on the mass-indexed isentropic coordinate $M_{{\mathit{\theta }}_e}$, which was first introduced by Jin et al. (2021). For each pair of transport models and flux products, we resolve cross-$M_{{\mathit{\theta }}_e}$ diabatic mixing rates and cross-$M_{{\mathit{\theta }}_e}$ APO gradients in the mid-latitudes of both hemispheres. We use observation-based diabatic mixing constraints diagnosed from four meteorological reanalyses, and observed APO gradient constraints calculated from three airborne campaigns. The detailed methodology for calculating $M_{{\mathit{\theta }}_e}$ surfaces, diabatic mixing rates, and cross-$M_{{\mathit{\theta }}_e}$ APO gradients is provided in Appendix C.

We show the climatological monthly mean diabatic mixing rates of two $M_{{\mathit{\theta }}_e}$ surfaces in the Southern Hemisphere in Fig. 7, as well as schematics of the geographic distribution of the two $M_{{\mathit{\theta }}_e}$ surfaces. For each ATM, mixing rates in Fig. 7 are calculated from APO and averaged over three realizations diagnosed from using three flux products. The reanalysis mixing rates are calculated from moist static energy (MSE) budget and shown as average and 1σ spread over the four reanalysis products. The six ATMs and the reanalyses show diabatic mixing rates with clear seasonal cycles, suggesting more rapid mixing across isentropes in the austral winter than summer. ATMs generally overestimate diabatic mixing rates, especially in the summer and winter, when there are large cross-$M_{{\mathit{\theta }}_e}$ APO gradients that lead to well-defined mixing rates. Among the six ATMs, CTE_TM5 and Jena_TM3 show too rapid mixing that is biased high in all seasons. The other four ATMs align better with reanalysis, but still show significant overestimation for most of the year. MIROC4-ACTM shows the best performance. These findings align with Jin et al. (2024), which previously identified that the southern hemisphere summer-time mixing rates are overestimated in ATMs used for CO₂ inversions, with consistent results for the three ATMs (MIROC4-ACTM, Jena_TM3, and CTE_TM5) being used in both studies.

Figure 7 Climatological monthly diabatic mixing rates across the (a) 30 and (b) 45 (10¹⁶ kg) M_θe surfaces in the Southern Hemisphere. ATM-diagnosed mixing rates are derived from six ATMs in APO-MIP1 that provide 3-D APO fields. Error bars represent the 1σ spread across the 30 and 45×10¹⁶ kg M_θe of three flux products used here. Black lines represent MSE-diagnosed mixing rates as the average of four reanalysis MSE budgets, while the gray shaded regions represent the 1σ spread. (c) Schematic showing latitude-pressure distribution of troposphere zonal annual average $M_{{\mathit{\theta }}_e}$, and (d) annual average near-surface $M_{{\mathit{\theta }}_e}$ contours of the 30 and 45 (10¹⁶ kg) surfaces, computed from MERRA-2 reanalysis for the year 2009. These two $M_{{\mathit{\theta }}_e}$ surfaces have very small seasonal meridional variability.

We find that biases in diagnosed diabatic mixing rates correlate with biases in cross-$M_{{\mathit{\theta }}_e}$ APO gradients in each season, with stronger diabatic mixing leading to smaller APO gradients (Fig. 8). Figure 8 shows the ATM-diagnosed diabatic mixing rates and simulated APO gradients (points) across six transport models and three flux products at two $M_{{\mathit{\theta }}_e}$ surfaces (30 and 45×10¹⁶ kg $M_{{\mathit{\theta }}_e}$) for three selected 2-month periods in the Southern Hemisphere. The points suggest clear linear relationships between diagnosed mixing rates and simulated APO gradients for each flux product (shown as fit lines for each flux product). The linear relationships persist across all seasons and $M_{{\mathit{\theta }}_e}$ surfaces, though with varying slopes depending on the underlying fluxes (Fig. 8). ATMs generally underestimate cross-$M_{{\mathit{\theta }}_e}$ absolute APO gradients (i.e., a closer to zero gradient) at both $M_{{\mathit{\theta }}_e}$ surfaces, corresponding to the overestimation of diabatic mixing rates in these models. For each flux product, biases in cross-$M_{{\mathit{\theta }}_e}$ APO gradients are always larger in fast mixing ATMs (e.g., Jena_TM3 and CTE_TM5) compared to slow mixing ATMs (e.g., two versions of NICAM-TM, MIROC4-ACTM, and CAM-SD), with MIROC4-ACTM showing the best agreement. For each transport model, the simulated gradient shows clear spread across different flux products. The largest spread occurs in austral winter and spring (Fig. 8c and d), when simulations with the DISS fluxes show much larger gradients compared to CESM or Jena fluxes. We note that the direct comparison of simulated and observed gradients for individual models is complicated by the interplay of ATM biases and flux product biases.

Figure 8 Using MSE-based diabatic mixing rates and airborne observations of cross-isentrope APO gradients to evaluate ATMs and flux models. Each panel compares model-diagnosed diabatic mixing rates (x-axis) and cross-$M_{{\mathit{\theta }}_e}$ APO gradients (y-axis) at the 30×10¹⁶ kg $M_{{\mathit{\theta }}_e}$ surface (a, c, e, ∼44° S surface outcrop) and at 45×10¹⁶ kg M_θe (b, d, f, ∼39° S surface outcrop). Results are shown for three seasonal periods: January–February (a, b), June–August (c, d), and October–November (e, f) based on available airborne campaigns. Points represent individual model simulations, with colors indicating flux products (Jena, CESM, DISS) and symbols denoting different ATMs. Vertical gray bands show the 1σ range of MSE-based mixing rates derived from four reanalysis products. Horizontal gray bands indicate the 1σ range of observed APO gradients after spatial and temporal bias correction. Colored lines show linear fits of mixing rates and APO gradients for each flux product across different transport models.

To evaluate flux products independently of transport model biases, we leverage both diabatic mixing rates and APO gradients. For each flux product, the intersection between the mixing rate-gradient linear fit and the MSE-diagnosed mixing rate indicates the expected APO gradient with realistic mixing characteristics. Therefore, we can evaluate large-scale flux features in the flux products by comparing this expected gradient to the observed gradient. Our analysis in Fig. 8 suggests that CESM is the most realistic flux product in the mid-latitude Southern Hemisphere in all seasons. The expected CESM gradients (intersections of thin blue line and vertical gray band) fall within the observation uncertainty range in all seasons and surfaces except austral summer at the 30×10¹⁶ kg $M_{{\mathit{\theta }}_e}$ surface (Fig. 8a), which suggests a slight underestimation of uptake in the CESM product. The expected gradients of the Jena flux product also generally fall within the observation uncertainty range, but shows an even larger underestimation in Fig. 8a. The expected gradients of the DISS flux product have large biases in the mid-latitude Southern Hemisphere. The expected gradient is significantly larger in the austral winter (Fig. 8c and d), and significantly smaller at the 30×10¹⁶ kg $M_{{\mathit{\theta }}_e}$ surface in austral summer (Fig. 8a) and austral spring (Fig. 8e), suggesting seasonal biases in the flux pattern.

Biases in expected gradients relative to observed gradients result from errors in the magnitude and spatial distribution of air-sea APO flux, specifically the difference in flux magnitudes between regions north and south of the target $M_{{\mathit{\theta }}_e}$ surface. For instance, a positive expected gradient bias during austral summer at the 30×10¹⁶ kg $M_{{\mathit{\theta }}_e}$ surface (Fig. 8a) in the DISS product could stem from underestimated outgassing in high southern latitudes, excessive outgassing in lower latitudes, or both. In addition, a flux product could produce realistic expected gradients despite underestimating absolute fluxes both north and south of the $M_{{\mathit{\theta }}_e}$ surface if the difference remains correct. Resolving these inherent ambiguities requires additional observational constraints from surface stations, ships, and aircraft, which we addressed in Sect. 3.1.

While the focus of Jin et al. (2024) was on the mid-latitude Southern Hemisphere, we extend our analysis of the mid-latitude diabatic mixing rates to the Northern Hemisphere at the 45×10¹⁶ kg $M_{{\mathit{\theta }}_e}$ surface (Fig. 9). ATMs also generally overestimate diabatic mixing rates in the Northern Hemisphere, except during summer (JJA). Whereas MSE-diagnosed mixing rates peak in northern summer, ATM-diagnosed mixing rates have their seasonal minimum at this time. We note that APO gradients in ATMs are close to zero during JJA, leading to poorly defined diabatic mixing rates. We carry out the same transport model and flux product analyses in the Northern Hemisphere in January to March (Fig. 10a) and August to October (Fig. 10b). MIROC4-ACTM still demonstrates the closest agreement with reanalysis data in both seasons, and CTE_TM5 shows the largest mixing rate bias. We note that TM3 and TM5 are based on similar parameterization schemes, but TM3 outperforms TM5. In both seasons, the expected gradients inferred from CESM flux align with the airborne observations, while Jena and DISS overestimate and underestimate expected gradients, respectively.

Figure 9 (a) Similar to Fig. 7, but showing climatological monthly diabatic mixing rates across the 45 (10¹⁶ kg) $M_{{\mathit{\theta }}_e}$ surface in the Northern Hemisphere. We note that JJA diabatic mixing rates in ATMs are poorly constrained due to close-to-zero cross-$M_{{\mathit{\theta }}_e}$ APO gradients. (b) Latitude-pressure distribution of zonal average 45×10¹⁶ kg $M_{{\mathit{\theta }}_e}$ surfaces during boreal summer (JJA) and winter (DJF). The two $M_{{\mathit{\theta }}_e}$ surfaces end at the tropopause, which is higher in the summer in the mid-latitudes. (c) Corresponding Earth surface outcrops of the JJA and DJF 45×10¹⁶ kg $M_{{\mathit{\theta }}_e}$ surfaces. Unlike in the Southern Hemisphere where seasonal meridional variations in $M_{{\mathit{\theta }}_e}$ surfaces are small, the Northern Hemisphere shows pronounced seasonal shifts due to different land/ocean heating and cooling cycles.

Figure 10 Similar to Fig. 8, but showing diabatic mixing rates and cross-$M_{{\mathit{\theta }}_e}$ APO gradients in the Northern Hemisphere late winter/early spring (a) and late-summer/early fall (b) of the 45×10¹⁶ kg $M_{{\mathit{\theta }}_e}$ surface. We choose January to March and August to October due to sufficient aircraft sampling and maximum cross-$M_{{\mathit{\theta }}_e}$ APO gradients in these months.

Our attempt to diagnose mixing rates in ATMs in the Northern Hemisphere mid-latitudes using ocean tracers alone is partly limited by the predominantly land surface. We find both summer and winter peaks in seasonal diabatic mixing rates in the northern mid-latitudes, driven by strong convection. Over land, convection peaks in summer due to strong surface heating that creates unstable atmospheric conditions. Over the ocean, however, convection peaks in winter due to larger air-sea temperature differences. Our ATM-diagnosed mixing rates in the Northern Hemisphere may not capture the summer peak because atmospheric mixing processes over land may not be adequately reflected in transport of air-sea APO flux signals, which occurs initially over the ocean. This limitation is particularly significant in the Northern Hemisphere, where zonal mixing is slower (2–4 weeks) due to topographic blocking and stationary wave patterns. We plan to diagnose the land and ocean contrast in atmospheric diabatic mixing in the next APO-MIP1 by also forward transporting land tracers (e.g., CO₂ sources/sinks from the land biosphere). Our method is more robust in the Southern Hemisphere mid-latitudes due to faster zonal mixing (1–2 weeks) and the predominantly ocean surface. We also note that the distinct thermal capacities of land and ocean in the Northern Hemisphere create more complex surface $M_{{\mathit{\theta }}_e}$ outcrops with larger latitudinal shifts across seasons (Jin et al., 2021), as shown in Fig. 9c. We, however, account for these shifts in our analysis.

Our analysis reveals that the ATM-diagnosed diabatic mixing rate primarily reflects an intrinsic characteristic of the transport model, at least in the Southern Hemisphere, showing little sensitivity to the underlying flux pattern, tracers, and land-ocean differences, particularly in models with smaller mixing rates (i.e., two versions of NICAM-TM, MIROC4-ACTM, and CAM-SD). These four models demonstrate consistent mixing rates across different flux products (Figs. 8 and 10). This consistency is further supported by our analysis of diagnosed mixing rates for individual APO components ($\mathrm{\Delta }{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ocn}}$, $\mathrm{\Delta }{{\mathrm{N}}_{\mathrm{2}}}^{\mathrm{ocn}}$, $\mathrm{\Delta }{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ocn}}$, $\mathrm{\Delta }{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ff}}$, and $\mathrm{\Delta }{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ff}}$) transported by ATMs with smaller mixing rates, which yields similar mixing rates despite these tracers having distinct signs, seasonal patterns, and magnitudes (Fig. S6). However, ATMs with faster mixing rate (e.g., Jena_TM3 and CTE_TM5) show large variability both across flux products (Figs. 9 and 10) and across tracers (Fig. S6). Notably, these two models exhibit approximately 50 % slower diagnosed mixing rates for the fossil fuel CO₂ tracer ($\mathrm{\Delta }{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ff}}$) compared to the other ocean flux tracers in the austral summer at the 30×10¹⁶ kg $M_{{\mathit{\theta }}_e}$ surface. We note that the fossil fuel CO₂ tracer has its main source in the Northern Hemisphere, and its mixing at the mid-latitude Southern Hemisphere preferentially occurs in the upper troposphere. In contrast, the air-sea flux tracers have significant sources/sinks over the Southern Ocean with rapid cross-isentrope mixing preferentially in the lower troposphere. This behavior suggests that these models simulate distinctly different mixing patterns between the planetary boundary layer (0–2 km) and the free troposphere. Specifically, these models appear to have excessive vertical mixing in the boundary layer while maintaining more realistic transport in the free troposphere. Our method, however, assumes a constant cross-$M_{{\mathit{\theta }}_e}$ diabatic mixing rate over the entire $M_{{\mathit{\theta }}_e}$ surface. The excessive boundary layer mixing causes the diagnosed mixing rates in these models to be overly sensitive to the specific vertical distribution of air-sea APO flux components.

Our evaluation of ATMs using simulations from APO-MIP1 advances the original framework of Jin et al. (2024) in two key aspects. First, we expand the experimental design by increasing the number of participating ATMs to six and employing three different flux fields with each ATM, generating 18 model realizations. This comprehensive matrix of simulations enables a more systematic evaluation of both transport and flux-related biases. We demonstrate how atmospheric tracer observations can be leveraged to independently evaluate and distinguish between biases in surface fluxes and atmospheric transport models. Second, we enhance the robustness of our MSE-diagnosed mixing rate calculations by incorporating two additional reanalysis products and computing mixing rates at the native high resolution of each reanalysis, rather than averaging to a coarser grid before the calculation. One limitation in our method is that we only use $M_{{\mathit{\theta }}_e}$ calculated from MERRA-2 for each of the transport models rather than using $M_{{\mathit{\theta }}_e}$ calculated from the individual transport model, which in principle can be done by interpolating the temperature and humidity from parent reanalysis to the ATM grid. This limitation would lead to slight inconsistency between the actual $M_{{\mathit{\theta }}_e}$ in the model and the value we assigned to it. However, the differences between $M_{{\mathit{\theta }}_e}$ calculated from different reanalyses remain small and our method ensures consistency in geography of each $M_{{\mathit{\theta }}_e}$ surface (Jin et al., 2021).

3.3 Shipboard model-observation comparison over the Drake Passage

The APO-MIP1 simulations could not reproduce latitudinal variations in APO seasonal cycle amplitude observed from shipboard measurements from 53 to 65° S over the Drake Passage and adjacent to Tierra del Fuego and the Antarctic Peninsula. Observations reveal a strong meridional SCA gradient (−2.1 per meg per degree, with deg positive northward), with SCA increasing sharply towards higher southern latitudes (Fig. 11). Model simulations substantially underestimate this latitudinal gradient (Fig. 11), showing weaker slopes averaged across ATMs of −1.2 (Jena), −0.5 (CESM), and 0.8 (DISS) per meg per degree. Notably, these gradients remain generally consistent across different ATMs for each flux product (±0.26, ±0.13, and ±0.29 per meg per degree, respectively), suggesting this may predominantly be a result of zonal-scale latitudinal biases in flux seasonality. Underrepresentation of enhanced summertime productivity along the coast of the Antarctic Peninsula in flux products could also play a role. However, the Gould typically only transits waters with elevated chlorophyll south of approximately 62° S while the gradient biases appear further north. Furthermore, seasonally, the SCA biases are caused more by underestimation of the winter/spring drawdown in APO at high latitudes, rather than the smaller underestimation of summertime APO enhancement (Figs. S10 and S11). For CESM, this bias could originate from incomplete process representation in the ocean biogeochemistry model and the underestimation of winter mixed-layer depths in the Pacific sector of the Southern Ocean, which has historically been a problem for Earth System Models (Sallée et al., 2013). The Jena flux product provides the closest match to the observed SCA gradient. However, several limitations remain, which likely stem from the coarse spatial resolution, limited atmospheric observational constraints over the Southern Ocean, and underrepresentation of mixing patterns around the PSA station (see details below and in the Supplement). The DISS flux product is biased due to its underlying assumptions and sparse observational constraints, as discussed in Jin et al. (2023).

Figure 11 Latitudinal distribution of APO SCA across the Drake Passage region (53–65° S) derived from ship observations and model simulations. We calculate SCA by grouping observations and model simulations into 1° latitude bands, shown as points. Model results are color-coded by ATM and organized by flux products in separate panels. The full seasonal cycles of observed and simulated APO of these latitude bands are shown in Fig. S10. We also show SCA observed and simulated for the PSA flask record as open crossed circles (∼64.5° S, shifted 0.7° S for visibility), and for ship data while the Gould is docked at or close to the PSA pier (left-most points, calculated by selecting data from 64.82 to 64.72° S and 64.1 to 64.0° W). The right-most three bands (53 to 55° S) are typically downwind of Tierra del Fuego (Figs. S7–S9). Both observational and model data for each latitude band or at PSA were detrended using corresponding cubic smooth spline fits from SPO. SCA was calculated using two-harmonic fits. The rightmost panel shows the SCA latitudinal gradients (per meg °⁻¹) from 53 to 65° S, with red shading indicating model biases relative to observations. The gradient is calculated as linear fits of SCA from 53 to 65° S for each ATM and flux product pair, and the observations. We exclude CAM-SD in this analysis because the ship data simulation is only available from 2012 to 2015 (i.e., missing 2016 to 2017 data).

Across ATMs, we find systematic differences of up to ±20 % in simulated mean SCA for the entire ship transects over the Drake Passage, independent of the input flux field, with CTE_TM5 consistently producing the smallest SCA and NICAM-TM_gl5 showing the largest. These differences across ATMs are likely caused by differences in marine boundary-layer ventilation in the models. Near-surface mixing over the Southern Ocean is challenging to model, owing to complex boundary-layer structure, strong wind shear, frequent storm systems, SST variations, and poorly represented clouds (Hyder et al., 2018; Knight et al., 2024; Lang et al., 2018; Truong et al., 2020). The coarse-resolution models used here may struggle to capture such phenomena, and the resulting variations in the concentration or dilution of flux signals near the surface drives differences in mean APO SCA. The systematic spread also likely reflects biases in the representation of large-scale diabatic mixing over the high southern latitudes. Models with strong diabatic mixing rates, such as TM5, tend to dilute the meridional gradient of seasonal amplitude through excessive mixing with lower-latitude air masses that have smaller SCAs, resulting in reduced amplitudes at high southern latitudes.

We find that observed SCA at PSA (64.5° S) from SIO flask measurements (∼70 per meg, averaged from 2012 to 2017) is significantly smaller than nearby ship data from 64 to 65° S (∼80 per meg). However, model simulations suggest similar values for both locations. The shipboard measurements are closely tied to the SIO O₂ calibration scale, and any remaining scale differences would be unlikely to affect the seasonal APO SCA. Rather, the observed SCA difference occurs because SIO flask samples collected at PSA predominantly sample descending air masses from the east that have passed over Anvers Island and the Antarctic Peninsula, with peaks above 2000 m (characterized by small APO SCA), whereas the ship samples marine boundary layer air including that over highly productive ocean regions (large APO SCA). As shown in Figs. S7–S9, the SIO flasks are collected from the Terra Lab, on the east side of the station, with a wind selection criteria of 5–205°. Even while docked at Palmer (left-most points in Fig. 11), the Gould measurements show elevated SCA compared to PSA flask samples, because the pier is located to the west of the station with samples filtered to exclude air influenced by the station (Figs. S7–S9). None of the ATMs, regardless of the flux product used, could reconstruct this feature, even though the models were sampled at the flask collection times. This difference is consistent with that seen between 900 mbar airborne samples and PSA flasks (Fig. 5e). The systematic bias points to the lack of resolution or physics that would be necessary, in either the reanalysis products or the ATMs, to accurately capture fine-scale circulation patterns, particularly the distinct air mass origins affecting ship versus station measurements. We note that the Jena flux product has been optimized to match seasonal APO cycles at Cape Grim Observatory (41° S) and at PSA (64.5° S), which may be the reason for its better performance on the SCA latitudinal gradient. It may do even better if the shipboard data were used in the inversion or if the effective sampling altitude of the SIO flasks at PSA were better accounted for.

Our analysis underscores the need for improvements in both ocean biogeochemistry models and ATMs. Future ocean process model developments should include improving accuracy of winter mixed-layer depths and higher-resolution ocean models with enhanced process representation to capture the fine-scale productivity patterns in the Southern Ocean. Additionally, current atmospheric transport models require improved resolution and physics to better represent the complex circulation patterns characteristic of coastal regions.

3.4 Implications for APO and CO₂ inversions and ATM development

Our study motivates a community effort to conduct APO inversions. Estimates of spatial and temporal variations in APO fluxes can improve our understanding of ocean biogeochemical processes and heat transport, and support verification of fossil-fuel emission estimates (Pickers et al., 2022; Rödenbeck et al., 2023). Currently, only one global-scale APO inversion product from Jena CarboScope (Rödenbeck et al., 2008) exists. This product shows excessive seasonal flux amplitudes (Fig. 2) in the southern low-latitudes (∼30 to 0° S) and northern mid-latitudes (∼30 to 60° N) relative to the other two flux products, which show better consistency with aircraft observations in their forward transport simulations (Fig. 6). These biases in Jena APO inversion partly result from limitations in the TM3 model, which exhibits excessive vertical mixing, particularly in the eastern North Pacific, too rapid diabatic mixing in the southern mid-latitudes, and underrepresentation of monsoon dynamics primarily due to coarse resolution (Jin et al., 2023). The large spread and biases in ATMs shown in this study highlight the importance of developing APO inversions using different ATMs and methodologies, as this will improve our ability to fully assess methodological uncertainties and potential biases in inverted air-sea APO flux estimates.

We encourage future inversion efforts to also assimilate column-mean data from airborne campaigns, in addition to sparse surface stations, especially for studying climatological seasonal fluxes. Our study finds that forward simulations from ATMs generally show large spread at northeastern Pacific sites, particularly at LJO and CBA (Fig. 2), where simulations are sensitive to model representation of the marine boundary layer and vertical mixing. The Scripps APO observation network consists mainly of stations along a Pacific transect close to the primary oceanic sources and sinks. Given this limited spatial coverage and our findings of significant vertical mixing biases (e.g., at CBA and LJO) and local wind-direction biases (e.g., at LJO and PSA) in ATMs at the station level, APO inversions that rely solely on these surface observations may be subject to large representation errors. Airborne data, however, provide larger surface footprints and column average metrics that are much less sensitive to vertical mixing biases. Our analysis shows that ATMs are generally consistent with each other in simulating large-scale annual and seasonal column-mean features along flight tracks (Fig. 6). Thus, inversions configured to assimilate airborne column-mean observations would be promising. Further improvement could also be achieved by incorporating shipboard observations to expand zonal coverage, such as from the Gould, across the Atlantic (Pickers et al., 2017), and in the Western Pacific (Tohjima et al., 2012). The study of Jin et al. (2023) used a different configuration of the Jena inversion that also assimilated Japanese ship-based observations across the western Pacific (Tohjima et al., 2012) from 40° S to 50° N. Forward transport of APO fluxes in that configuration aligns better with station and airborne data compared to the configuration used in this study, particularly in reducing the SCA bias in the tropics, suggesting better flux representations.

Biases in diabatic mixing diagnosed from ATMs (Sect. 3.2) imply that CO₂ inversions using these ATMs are also likely biased. A previous study showed that summer-time Southern Ocean CO₂ estimates from inversion products are correlated with corresponding simulated summer-time cross-isentrope CO₂ gradients in inversions (Long et al., 2021b). The simulated gradients are shown to be biased too small due to too rapid diabatic mixing bias in ATMs leading to an overestimation of Southern Ocean CO₂ uptake in the summer (Jin et al., 2024). It is likely that biases in ATMs also contribute to the large spread found in OCO-2 MIP and Global Carbon Project (GCP) inversion ensembles (Byrne et al., 2023; Crowell et al., 2019; Friedlingstein et al., 2025; Peiro et al., 2022). We identify several priority areas for understanding biases in ATMs, particularly the inconsistency between diabatic mixing rates diagnosed from the MSE budgets of parent reanalysis and the tracer fields of coarser resolution ATMs identified here. These inconsistencies likely stem from several potential sources: (1) regridding of original reanalyses to the coarser resolution of the ATM grid, (2) for online GCMs using nudging, incomplete matching of the input meteorology, and (3) for offline models, recalculation or parameterization of convective mass fluxes in the coarser ATM. The first potential source of error from regridding could be evaluated by comparing MSE-based diabatic mixing rates from the parent and regridded fields as long as all components of MSE were included in the regridding. The second potential source of error from nudging could be evaluated by comparing MSE-based diabatic mixing rates from the regridded parent model and the nudged online simulation. Finally, the third potential source of error from recalculating or parameterizing vertical mass fluxes could be evaluated by comparing the MSE-based diabatic mixing rates from the regridded parent model and the tracer-based mixing rates from the ATM. It is notable that diabatic mixing rates diagnosed from two online models, MIROC4-ACTM and CAM-SD, which do not require regridding, are generally consistent with observations, with MIROC4-ACTM showing the best performance among all models (Figs. 7–10).

An important consideration is that the real atmosphere mixes MSE and tracers at different spatial and temporal scales. In the Northern Hemisphere, APO fluxes initially mix vertically over oceans, while strong CO₂ fluxes initially mix vertically over land. In contrast, MSE fluxes mix initially over both land and ocean. Due to the large land area in the Northern Hemisphere, the zonal mixing time scale is much longer (∼2–4 weeks) so that diabatic mixing rates diagnosed from APO or CO₂ tracers could differ from each other and from those diagnosed from MSE tracers. In the Southern Hemisphere mid-latitudes, these potential differences are much smaller due to the predominance of ocean and rapid zonal mixing (∼1–2 weeks). In general, the timescales for diabatic mixing are longer than the timescales of zonal mixing, which support our approach of using tracer fluxes over both ocean and land to evaluate zonal-mean diabatic mixing. Future work should also develop metrics for quantifying along-isentrope (adiabatic) transport to complement our understanding of tracer mixing across isentropes. The timescales of adiabatic mixing influences tracer gradients along isentropic surfaces, which in turn affects diabatic mixing differently in the upper versus lower troposphere. It is also necessary to examine the sensitivity of mixing rates to model resolution, particularly vertical levels at the interface between the boundary layer and free troposphere, and boundary layer schemes. These ATM improvements are essential for enhancing both forward simulations and inverse estimates of surface fluxes.

4 Summary and outlook

We conducted the Atmospheric Potential Oxygen forward Model Intercomparison Project (APO-MIP1) to generate forward simulations of APO and its components using different flux products and eight ATMs. This effort provides model APO simulations at surface stations, along aircraft flight paths, and on ships that can be directly compared with observations. Additionally, we provide 3-D APO fields from six of the eight ATMs. We use simulations from APO-MIP1 to evaluate eight ATMs and three flux products by comparing simulations against observations from surface stations, aircraft, and ships.

We find that model simulations of APO seasonal cycles using a given flux product show considerable summer-time spread at northern surface stations, particularly at two eastern Pacific stations, LJO and CBA (Fig. 3). The bias stems from challenges in accurately representing complex atmospheric vertical transport processes, marine boundary layer mixing, and coastal horizontal mixing in these regions. These findings highlight the limitations of current APO inversions that rely on a single ATM (i.e., TM3 used in Jena APO inversion) and sparse surface observations. However, model simulations of column-average APO resolved from sampling aircraft tracks are consistent across different ATMs, emphasizing the importance of airborne measurements for constraining large-scale flux features.

Using airborne observations and a moist-isentropic coordinate framework, we demonstrate that most ATMs overestimate diabatic mixing rates in the mid-latitudes of both hemispheres when compared to mixing rates derived from energy budgets of reanalyses. Among all ATMs used here, Jena_TM3 and CTE_TM5 show the largest biases. These constraints also enable us to separate flux biases from transport-related biases, allowing independent evaluation of flux models, which show that the CESM flux product is the best among the three flux products used in this study. This prognostic model outperforms two observation based products because of sparse atmospheric and surface observations, limitations in ATM used in atmospheric inversion, and because seasonal APO fluxes are driven by physical and biological processes that CESM represents well.

We encourage the broader community to develop new APO inversions, which could provide independent constraints on ocean biogeochemical processes and improve our understanding of the ocean carbon sink. Model simulations from APO-MIP1 can be used in other applications, including the calibration of methods for estimating seasonal air-sea APO fluxes from global atmospheric observations (e.g., Jin et al., 2023), constraining the representation of regional to global marine production in Earth system models (e.g., Nevison et al., 2012, 2015, 2018), and for understanding ESM biases in seasonal air-sea CO₂ exchange related to both thermal and non-thermal forcings. The transport simulations can also support the evaluation of long-term trends in O₂ : CO₂ ratios over the Southern Ocean based on surface station gradients, useful for assessing biogeochemical responses to climate change.

We expect APO-MIP1 to continue evolving as an active collaboration examining atmospheric tracer transport and air-sea O₂ flux estimates. The current implementation excluded the air-sea CO₂ component and long-term flux trends from the Jena flux product, and does not include interannual and long-term flux trends in the DISS flux product, making these simulations unsuitable for interpreting interannual to long-term air-sea O₂ fluxes features. Thus, we only analyze APO seasonal cycles and meridional gradients here. The next phase of APO-MIP1 will address these limitations by incorporating updated inversion flux fields based on a larger set of atmospheric APO observations and including interannual variability. We will expand the scope by including terrestrial O₂ flux fields for O₂-specific analyses and seasonal-only component fluxes to investigate rectifier effects. The seasonal rectifier effect refers to the creation of non-zero annual mean atmospheric concentration gradients at surface stations even with balanced seasonal O₂ fluxes. This occurs when fluxes correlate with seasonal variations in atmospheric mixing. For example, strong summer O₂ outgassing combined with shallow PBL heights concentrates APO near the surface, while higher winter PBL dilutes the O₂ uptake signal, resulting in observed annual mean APO gradients even when the annual mean flux is zero. Additionally, we plan to update air-sea O₂ fluxes derived from surface ocean dissolved oxygen measurements by replacing Garcia and Keeling (2001) with fluxes calculated from recent machine learning interpolation of dissolved oxygen products (Gouretski et al., 2024; Ito et al., 2024; Sharp et al., 2023). We encourage broader participation from diverse modeling groups in the next phase of APO-MIP1.

Appendix A: Surface station, airborne, and shipboard APO measurements

The surface station APO observations from the Scripps O₂ program have been described in (Keeling et al., 1998). Briefly, flask triplicates have been collected at biweekly to monthly frequency during clean background air conditions at a network of sites for over three decades, and returned to Scripps for analysis using interferometric and mass-spectrometric techniques. Here we use monthly data that was averaged from roughly bi-weekly data. The flask measurements are first adjusted to the middle of each month, parallel to the mean seasonal cycle for that station, before averaging. The APO-MIP1 output for these stations was reported matching the ObsPack CO₂ files from the Scripps O₂ Program, to take advantage of the established ObsPack format. These CO₂ measurements correspond to the same flask air on which O₂ is measured. The model output is treated in the same way as the observations to generate monthly means.

Airborne APO measurements from HIPPO, ORCAS, and ATom campaigns were made in situ with the NSF NCAR Airborne Oxygen Instrument (AO2), using a vacuum-ultraviolet absorption technique to measure O₂ and a single-cell infrared gas analyzer to measure CO₂ (Stephens et al., 2021g). AO2 produces measurements every 2.5 s, which are averaged to 10 s frequency for merging with other aircraft data. To correct for flight-specific sampling offsets, the in situ AO2 data were adjusted to agree with flask measurements collected during each flight using the NSF NCAR/Scripps Medusa flask sampler on a flight-by-flight average basis (Jin et al., 2023; Stephens et al., 2021g).

HIPPO and ATom had nearly pole-to-pole coverage, and from near surface (150–300 m) to above the tropopause. HIPPO consisted of five campaigns between 2009 and 2011, and most data were collected above the Pacific. ATom consisted of four campaigns between 2016 and 2018, and each campaign had a Pacific transect and an Atlantic transect. ORCAS was a 6-week campaign with dense temporal sampling over the Drake Passage and ocean areas adjacent to the tip of South America and the Antarctic Peninsula. The APO-MIP1 output for these aircraft measurements was reported matching the ObsPack CO₂ files for each campaign. These data are also at 10 s frequency but correspond to different instruments with different calibration intervals. To match the observed and model time series, we mask observations when model output is not available, and vice versa. We also exclude any stratospheric data, with the stratosphere defined as water vapor concentrations below 50 ppm and either ozone concentrations exceeding 150 ppb, or detrended N₂O levels (normalized to 2009) below 319 ppb (Jin et al., 2021). Water vapor and ozone were measured by the NOAA UAS Chromatograph for Atmospheric Trace Species instrument (Hintsa et al., 2021). N₂O was measured by the Harvard Quantum Cascade Laser System instrument (Santoni et al., 2014). We filter the airborne data to exclude continental or urban boundary-layer air sampled while landing, taking off, or conducting missed approaches at airports (Jin et al., 2021).

Shipboard APO measurements from the ARSV L. M. Gould were made in situ during over 90 transects of Drake Passage on 50 cruises between 2012 and 2017 using a fuel-cell method for O₂ and a two-cell non-dispersive infrared gas analyzer for CO₂. The instrumentation was similar to a previously developed tower system (Stephens et al., 2003), but adapted and optimized for shipboard use. The instrument produces measurements at 1 min frequency. The cruises occurred in all months of the year but are more sparse during austral winter. The Gould operated almost exclusively between Punta Arenas, Chile and Palmer Station, Antarctica, in support of resupplying and transferring personnel to Palmer Station. The cruises span from 53 to 65° S in all months, and extend as far as 70° S during summer months. The APO-MIP1 output for the Gould was reported matching the ObsPack CO₂ file from the NOAA underway pCO₂ system. This system measures atmospheric CO₂ for 15 min every 2 h. To match the observed and model time series, we first calculate hourly means for each and then mask observations when model output is not available, and vice versa.

The resolved APO annual mean and seasonal cycles have negligible measurement uncertainty compared to model spread because we average data over long time series for stations and over large spatial domains for aircraft and ships, effectively reducing the already small short-term instrument imprecision.

Appendix B: APO flux products

B1 Air-sea APO flux products

The first air-sea APO flux product (Jena) is air-sea APO flux from the Jena CarboScope APO Inversion (version ID: apo99X_v2021), which is available directly as $F_{\mathrm{APO}}^{\mathrm{ocn}}$ (update of Rödenbeck et al., 2008). In this inversion, the posterior fluxes (variable name: apoflux_ocean) were optimized to best match observed APO at 9 stations in the Scripps O₂ Program surface network (Manning and Keeling, 2006) and at 2 stations from the National Institute for Environmental Studies (Tohjima et al., 2012). The prior air-sea CO₂ flux was not included in the forward simulations here. We note that the exclusion of prior air-sea CO₂ flux has only minimal impact on the simulated APO seasonal cycle and north-to-south annual gradient but reduces the tropical “bulge” of annual mean by approximately 1 per meg and results in close to zero long-term APO trend. The Jena product is available from 1999 to 2020 originally with spatial resolution of 2° latitude × 2.5° longitude at daily intervals, converted to 1°×1°. The Jena inversion used the TM3 transport model, which is also one of the models participating in APO-MIP1. In the case of TM3 forward transport simulation, the Jena inversion posterior fluxes have been re-run forward through the ATM, and thus this combination of fluxes and transport should agree well at the surface stations used for inversion optimization.

The second air-sea APO flux product (CESM) uses air-sea O₂, CO₂, and N₂ flux components from the Community Earth System Model (CESM2) Forced Ocean-Sea-Ice (FOSI) simulation (Yeager et al., 2022), which is forced by atmospheric fields from JRA55-do reanalysis (Tsujino et al., 2018) and prognostic ocean biogeochemistry using the Marine Biogeochemistry Library (MARBL, Long et al., 2021a). The model directly produces $F_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ocn}}$ and $F_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ocn}}$ , while $F_{{\mathrm{N}}_{\mathrm{2}}}^{\mathrm{ocn}}$ is calculated by scaling the ocean heat flux (Q, W m⁻²) output using the relationship from Keeling and Shertz (1992) following

$$\begin{array}{}\text{(B1)}& F_{{\mathrm{N}}_{\mathrm{2}}}^{\mathrm{ocn}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{1.3}}}\cdot {\displaystyle \frac{\mathrm{d}S}{\mathrm{d}T}}\cdot {\displaystyle \frac{Q}{C_p}},\end{array}$$

where $\mathrm{d}S/\mathrm{d}T$ (mol kg⁻¹ °C⁻¹) is the temperature derivative of solubility using solubility coefficients from Hamme and Emerson (2004). C_p represents the specific heat capacity of seawater, which is assumed to be 3993 J kg⁻¹ °C⁻¹. The factor of $\mathrm{1}/\mathrm{1.3}$ is to adjust the seasonal amplitude due to the temporal lag between tracer flux and heat flux, as proposed by Jin et al. (2007).

These three CESM flux components have a resolution of 1° latitude × 1° longitude grid with the North Pole displaced to Greenland. All fields are available from 1958 to 2020, but we only use fluxes from 1986 to 2020. $F_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ocn}}$ and $F_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ocn}}$ are output from the model at daily resolution, whereas $F_{{\mathrm{N}}_{\mathrm{2}}}^{\mathrm{ocn}}$ is calculated from monthly model heat fluxes then interpolated to daily resolution. This version of CESM was designed to initialize a seasonal-to-multiyear large ensemble (SMYLE) of coupled simulations for evaluating predictability. It is forced by observed meteorology starting in 1958, at which point it branches off of a FOSI configuration using JRA55-do atmospheric fields as surface boundary conditions (Yeager et al., 2022). The FOSI simulation consists of six consecutive cycles of 1958–2018 forcing, with the sixth cycle (used for SMYLE) extended through 2020. Annual mean heat fluxes from this configuration show a small cooling drift over the historical period, and thus the inferred annual mean and long-term trend of O₂ and N₂ flux should not be interpreted as realistic.

The third air-sea APO flux product (DISS) uses bottom-up air-sea O₂ and CO₂ flux estimates derived primarily from dissolved gas measurements. $F_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ocn}}$ consists of a seasonal component calculated from the dissolved O₂ measurement based climatology of Garcia and Keeling (2001), with seasonal amplitude scaled by 0.82 according to Naegler et al. (2006), and an annual mean component from the ocean inversion of Resplandy et al. (2016) for 21 regions using transport from MITgcm-ECCO. Bent (2014) reported that the 0.82 scaling factor significantly improved agreement between GK flux and HIPPO observations, based on simulations using one ATM (a different MIROC4-ACTM configuration). However, our results show that applying this 0.82 scaling factor actually leads to an underestimation of modelled column-mean APO SCA when comparing with the combined HIPPO, ORCAS, and ATom observations at high latitudes in both hemispheres. The seasonal component (1.125°×1.125° × monthly) was linearly regridded to 1°×1° × daily resolution. For the annual mean component, the original regional values (21 regions) were spatially interpolated to 1°×1° resolution while conserving the total sum within each region, then temporally interpolated to daily values. We use $F_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ocn}}$ from the machine learning interpolation of pCO₂ based air-sea CO₂ fluxes (Jersild et al., 2017; Landschützer et al., 2016). The version of this product that we used provides fluxes from 1982 to 2020, with resolution of 1° latitude × 1° longitude × monthly, which we interpolated to daily. We use Eq. (B1) to calculate $F_{{\mathrm{N}}_{\mathrm{2}}}^{\mathrm{ocn}}$ with heat fluxes from ERA5 reanalyses (Hersbach et al., 2020), which is available from 1979 onwards, with resolution of 0.25° latitude × 0.25° longitude × monthly. Sea-surface temperature (SST) estimates required to calculate $\mathrm{d}S/\mathrm{d}T$ (Eq. B1) are from World Ocean Atlas (WOA) v2018 with resolution of 1° latitude × 1° longitude × monthly. SST is available as a 1981 to 2010 climatology but we use it repeatedly for 1986 to 2020.

B2 Fossil fuel APO uptake products

We used two products for $F_{\mathrm{APO}}^{\mathrm{ff}}$. The first product (GridFED) uses fossil CO₂ emission and O₂ uptake fluxes from Jones et al. (2021), downloaded from Jones et al. (2022). This product is available from 1959 to 2020, with resolution of 0.1° latitude × 0.1° longitude × monthly, which we interpolate to daily.

The second product (OCO2MIP) use $F_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ff}}$ as prepared for the OCO-2 Model Intercomparison Project (MIP) version 10, downloaded from Basu and Nassar (2021), with resolution of 1° latitude × 1° longitude × hourly. This $F_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ff}}$ product uses fossil fuel CO₂ emission from ODIAC (Oda et al., 2018) for 2000 to 2019. For 2020, the flux was scaled from 2019 using the ratio of 2020 to 2019 global emissions reported by Liu et al. (2020). $F_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ff}}$ is not available from this product, but we scale the atmospheric field of $\mathrm{\Delta }{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ff}}$ by a factor of −1.4 to estimate $\mathrm{\Delta }{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ff}}$ (Keeling, 1988; Steinbach et al., 2011). We primarily use GridFED, except for CAMS_LMDZ where we use OCO2MIP instead, because $F_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{ff}}$ from GridFED is missing for years after 2015. The differences between these two products are negligible compared to the magnitude of ocean-driven APO variations, for the seasonal metrics considered here.

Appendix C: Calculation of $M_{{\mathit{\theta }}_e}$, cross-$M_{{\mathit{\theta }}_e}$ diabatic mixing rates and APO gradients

The mass-indexed moist isentropic coordinate $M_{{\mathit{\theta }}_e}$ is defined as the total dry air mass under a specific moist isentropic surface (θ_e) in the troposphere of a given hemisphere. Surfaces of constant $M_{{\mathit{\theta }}_e}$ are parallel to surfaces of constant θ_e but the relationship changes with season, as the atmosphere warms and cools. $M_{{\mathit{\theta }}_e}$ surfaces have air mass (10¹⁶ kg) as the unit, and are adjusted to conserve dry air mass below the surface at any instant in time. $M_{{\mathit{\theta }}_e}$ is calculated as a function of θ_e and time following

$$\begin{array}{}\text{(C1)}& M_{{\mathit{\theta }}_e}\left(xt\right)=\sum M_x\left(t\right)\mathrm{|}{\mathit{\theta }}_{e_x}<{\mathit{\theta }}_e,\end{array}$$

where x indicates an individual grid cell of the atmospheric field, M_x(t) is the dry air mass of each grid cell x at time t, and ${\mathit{\theta }}_{e_x}$ is the equivalent potential temperature of the grid cell. For a given θ_e threshold, the corresponding $M_{{\mathit{\theta }}_e}$ value is calculated by integrating the air mass of all grid cells with θ_e value smaller than the threshold. We only integrate air mass in the troposphere, which is defined here as potential vorticity unit (PVU) smaller than 2. At each time step, this calculation yields a unique value of $M_{{\mathit{\theta }}_e}$ for each value of θ_e as well as a 3-D field of atmospheric $M_{{\mathit{\theta }}_e}$. Following the spatial pattern of θ_e, $M_{{\mathit{\theta }}_e}$ values generally increase from low to high altitudes and from poles to equator. We generate daily $M_{{\mathit{\theta }}_e}$ fields using four different reanalysis products (MERRA-2, JRA-55, JRA-3Q, and ERA5) at their native resolution, avoiding potential information loss from grid interpolation (Gelaro et al., 2017; Hersbach et al., 2020; Kobayashi et al., 2015; Kosaka et al., 2024).

The calculation of diabatic mixing rates in ATMs is based on a box model approach, which uses $M_{{\mathit{\theta }}_e}$ as boundaries. A schematic of the box model is available as Fig. 1 of Jin et al. (2024). The box model invokes tracer air mass balance, which recognizes tracer inventory change (M_i, Tmol) of each $M_{{\mathit{\theta }}_e}$ box equal to the sum of surface fluxes (F_i, Tmol d⁻¹) and the diabatic transport between boxes ($T_{i,i+\mathrm{1}}$, Tmol d⁻¹, positive poleward). The transport term is considered as a diffusive system, which is parameterized as the product of diabatic mixing rate across the $M_{{\mathit{\theta }}_e}$ boundary ($D_{i,i+\mathrm{1}}$, (10¹⁶ kg)² d⁻¹) and the tracer concentration (χ_i+1, Tmol tracer per kg air mass) gradient between two boxes. The full mass balance follows

$$\begin{array}{}\text{(C2)}& {\displaystyle \frac{\partial M_i}{\partial t}}=\left\{\begin{array}{ll}{F}_i+T_{i,i+\mathrm{1}}& \mathrm{if}i=\mathrm{1}\\ F_i+T_{i,i+\mathrm{1}}-T_{i-\mathrm{1},i}& \mathrm{if}i>\mathrm{1}\end{array}\right.,\end{array}$$

with

$$\begin{array}{}\text{(C3)}& T_{i,i+\mathrm{1}}=D_{i,i+\mathrm{1}}\cdot {\displaystyle \frac{{\mathit{\chi }}_{i+\mathrm{1}}-{\mathit{\chi }}_i}{\mathrm{\Delta }{M}_{{\mathit{\theta }}_e}}}.\end{array}$$

In these equations, i is the number label of the box and is set to be 1 at the highest latitude, $\mathrm{\Delta }{M}_{{\mathit{\theta }}_e}$ is the distance in $M_{{\mathit{\theta }}_e}$ coordinates between box centers, which for evenly spaced boxes as used here, is the same as the total air mass of each box. In this study, we set the range of each $M_{{\mathit{\theta }}_e}$ box to be 15×10¹⁶ kg air mass, and therefore $\mathrm{\Delta }{M}_{{\mathit{\theta }}_e}$ equals the same value. The diabatic mixing rate (D) can be expressed as

$$\begin{array}{}\text{(C4)}& D_{i,i+\mathrm{1}}\left(t\right)={\displaystyle \frac{\left[\sum _{i^{\prime }=\mathrm{1}}^{i^{\prime }=i}\left(\frac{\mathrm{d}{M}_{i^{\prime }}\left(t\right)}{\mathrm{d}t}-F_{i^{\prime }}\left(t\right)\right)\right]}{\left[{\mathit{\chi }}_{i+\mathrm{1}}\left(t\right)-{\mathit{\chi }}_i\left(t\right)\right]}}\cdot \mathrm{\Delta }{M}_{{\mathit{\theta }}_e}.\end{array}$$

This method effectively reconstructs large-scale tracer transport features (T) in ATMs, as demonstrated in Jin et al. (2024). We note that the diabatic mixing rate is a property of the corresponding $M_{{\mathit{\theta }}_e}$ and is theoretically insensitive to the choice of box sizes. We calculate climatological monthly average (2009 to 2018) diabatic mixing rates for each of the six transport models using the 3-D APO fields from transporting each of the three flux products (Figs. 7 and 9). To assign $M_{{\mathit{\theta }}_e}$ at the model grid locations and times for each ATM, we always use $M_{{\mathit{\theta }}_e}$ from MERRA-2 interpolated to the ATM grid, to ensure spatial consistency. Using other reanalyses only leads to small (<5 %) differences in ATM-diagnosed diabatic mixing rates (Jin et al., 2024).

Independent observational constraints on ATM-diagnosed mixing rates are calculated from moist static energy (MSE) budgets of four meteorological reanalyses (Figs. 7 and 9). MSE is a measure of static energy that is conserved in adiabatic ascent/descent and during latent heat release due to condensation, and naturally aligns with surfaces of θ_e or $M_{{\mathit{\theta }}_e}$. This diagnostic approach offers more robust mixing rate estimates than tracer-based methods in part because MSE maintains consistent, non-zero gradients at each reanalysis time step, unlike chemical tracers. Additionally, MSE-based mixing rates are directly diagnosed from reanalysis on the original grid, avoiding potential artifacts introduced when these fields are interpolated to coarser transport model grids, and any recalculation of vertical mass fluxes and subgrid-scale mixing parameterizations in ATMs.

The MSE-diagnosed mixing rate calculation adapts our tracer box model framework. In this adaptation, we replace tracer inventory (M_i, Tmol) by MSE (S_i, J), replace surface tracer flux (F_i, Tmol d⁻¹) by surface heat flux (Q_i, J d⁻¹), and add an additional term to account for atmospheric radiative energy balance (R_i, J d⁻¹), following

$$\begin{array}{}\text{(C5)}& D_{i,i+\mathrm{1}}\left(t\right)={\displaystyle \frac{\left[\sum _{i^{\prime }=\mathrm{1}}^{i^{\prime }=i}\left(\frac{\mathrm{d}{S}_{i^{\prime }}\left(t\right)}{\mathrm{d}t}-Q_{i^{\prime }}\left(t\right)-R_{i^{\prime }}\left(t\right)\right)\right]}{\left[{\mathit{\chi }}_{i+\mathrm{1}}\left(t\right)-{\mathit{\chi }}_i\left(t\right)\right]}}\cdot \mathrm{\Delta }{M}_{{\mathit{\theta }}_e}.\end{array}$$

We note that the gradient on the denominator in Eq. (C5) represents the MSE density gradient (J kg⁻¹ air mass) across the $M_{{\mathit{\theta }}_e}$ surface. The calculation of these terms requires air temperature, specific humidity, surface heat flux, including surface sensible and latent heat flux, and radiative imbalance from reanalysis. Further details on the process to diagnose mixing rate from both ATMs and reanalyses can be found in Jin et al. (2024).

The cross-$M_{{\mathit{\theta }}_e}$ APO gradient was calculated using data grouped into two adjacent boxes in the $M_{{\mathit{\theta }}_e}$ space, with box centers spanning 15×10¹⁶ kg air mass across the target surface boundary. For each box, we calculate the average APO concentration by trapezoidal integration of detrended APO as a function of $M_{{\mathit{\theta }}_e}$ and dividing by the $M_{{\mathit{\theta }}_e}$ range (Jin et al., 2021). We carry out the calculation for each airborne campaign, using the observations, model flight track output, and 3-D model fields. Flight-track estimated cross-$M_{{\mathit{\theta }}_e}$ APO gradients are not directly comparable to simulated gradients from full 3-D fields, due to spatial and temporal coverage biases in airborne observations. We correct for both biases in the APO airborne observations and model flight track output (detailed in Sect. S1).