Introduction

The global ocean has experienced a loss of oxygen, or deoxygenation, since the mid-twentieth century (Helm et al., 2011;Ito et al., 2017;Keeling et al., 2010), and will likely continue to deoxygenate as global warming persists (Bopp et al., 2013;Kwiatkowski et al., 2020). Insufficient oxygen availability, or hypoxia, can impair growth and reproduction in marine organisms, perturb predator-prey dynamics, and lead to fish mortality events (McCormick & Levin, 2017;Rabalais et al., 2002;Vaquer-Sunyer & Duarte, 2008). Oxygen availability is already a limiting factor for marine life, making much of the ocean interior inhospitable for macrofauna (Bertrand et al., 2011;Pörtner & Knust, 2007). Future ocean deoxygenation is expected to further restrict marine habitats (Bindoff et al., 2019;Deutsch et al., 2020;Morée et al., 2023;Seibel, 2011), force migrations of marine species (Pinsky et al., 2013), and reduce the resiliency of ecosystem function and services (Breitburg, 2002;Gattuso et al., 2015;Sandifer & Sutton-Grier, 2014).

Oxygen in the ocean is supplied by water masses that follow ventilation pathways connecting the ocean surface and interior (MacGilchrist et al., 2017;Marshall et al., 1993;Sallée et al., 2010;Wyrtki, 1962). Concentrations of oxygen are approximately saturated at the surface and gradually depleted by biological respiration along these pathways (Redfield, 1942). This results in the formation of high-and low-oxygen regions, including oxygen minimum zones (OMZs) at thermocline depths (about 100-1,000 m) in tropical (30°S-30°N) oceans where the renewal of waters from the surface is slow (Karstensen et al., 2008;Luyten et al., 1983b). Global warming reduces oxygen solubility and slows down the main ventilation pathways supplying oxygen to the ocean interior (Bopp et al., 2017;Gnanadesikan et al., 2007;Oschlies et al., 2018;Palter & Trossman, 2018). This oxygen loss is of particular concern in tropical oceans where oxygen levels are already critically low (Stramma et al., 2012).

Numerous studies based on Earth System Model (ESM) projections have shown a consistent signal of global deoxygenation with warming, but have reported uncertainty as to whether tropical OMZs will intensify in the future (Bopp et al., 2013(Bopp et al., , 2017;;Cabré et al., 2015;Cocco et al., 2013;Kwiatkowski et al., 2020;Resplandy, 2018). Recent studies have, however, found consistencies in OMZ projections across ESMs, demonstrating that lowest oxygen waters consistently, and counter-intuitively, gain oxygen with warming, whereas the outer and more oxygenated layers of OMZs lose oxygen with warming (Busecke et al., 2022;Ditkovsky et al., 2023). This pattern was found using an oxygen-coordinate approach-sorting the tropical ocean from lowest to highest oxygen levels-which isolates water masses according to their progression along ventilation pathways (younger waters close to the surface have higher oxygen concentrations). A recent study by Portela et al. (2024) also finds that the use of an oxygen-coordinate can simplify the description of global ventilation pathways. Indeed, ESMs capture the main water masses and ventilation pathways in the ocean, but their locations may not exactly align across the different models or with observed patterns in geographic-space. There are thus limitations in comparing observations, model historical simulations and projections strictly in a geographic-space, but also opportunities offered by physically or biologically informed frameworks, particularly those that use an oxygen-based coordinate (Busecke et al., 2022;Ditkovsky et al., 2023;Portela et al., 2024).

Here, we introduce an "oxygen water mass framework" to constrain and understand projections of ocean oxygen and OMZs in an ensemble of 14 ESMs from the Coupled Model Intercomparison Project phase 6 (CMIP6; Eyring et al., 2016) that include oxygen. This framework separates oxygen changes into two components: transformation of water masses in "ventilation-space," and redistribution of water masses in "geographic-space." Within this framework, we compare two subsets of ESMs: those that best reproduce the geography of oxygen water masses in the thermocline and those that poorly represent it. We find consistent global patterns of transformations in ventilation-space across all ESMs, with deoxygenation concentrated in well-ventilated waters, whereas OMZs remain stable, but find that the high-skilled model group systematically projects lower levels of deoxygenation. The high-skilled model group also shows higher precision in future projections of water mass redistribution, clarifying the sign of oxygen change over much of the tropical ocean thermocline.

Methods and Materials

Oxygen Water Mass Framework: Transformation in Ventilation-Space and Redistribution in Geographic-Space

Future oxygen changes at a fixed geographic location can arise in two ways (see schematic in Figure 1): (a) water masses can be transformed by physical, chemical, or biological processes, such as air-sea exchanges, biological respiration, or turbulent mixing between water masses (Groeskamp et al., 2016;Walin, 1982), and (b) water masses can be redistributed in space and reshaped by processes including (but not limited to) changes in ocean circulation, which results in strong local trends that sum to zero globally (Zika et al., 2021). Here, we define water masses in a "ventilation-space" by the oxygen-percentile relation, O 2 * Öp, tÜ, which tracks the percent of ocean volume, p, with oxygen concentrations less than O 2 * at time t. Oxygen is a useful tracer to define water masses and characterize their relative age along ventilation pathways, that is, the average time since that water mass has been exposed to the surface (e.g., Richards & Redfield, 1955;Talley, 2011).

where O 2 Öx, y, z, tÜ is the geographic field of dissolved oxygen. We then rescale and invert this function to obtain an oxygen-percentile relationship.

where O * 2 Ö p, tÜ is the time-dependent cumulative distribution function of oxygen with respect to percentile, p, and V T is the total ocean volume. We compute O * 2 Öp, tÜ numerically using bins of width 1% of V T . The first bins group the lowest oxygen values and correspond to the least ventilated water masses of the ocean, whereas the last bins group the highest oxygen values and correspond to the best ventilated water masses of the ocean. The cumulative oxygen content below a given percentile, Φ O 2 Öp, tÜ, is computed by cumulatively integrating O 2 * Öp, tÜ over p (Figure 1b). Ocean processes transform water masses by changing O * 2 Öp, tÜ, or equivalently Φ O 2 Öp, tÜ. In addition to boundary fluxes ÖFÜ and mixing processes ÖMÜ which transform water masses defined by conservative tracers (e.g., conservative temperature; Sohail et al., 2021), we must also consider internal biogeochemcial sources and sinks ÖSÜ when tracking the transformation of oxygen water masses (see Appendix A for details).

In ventilation-space, oxygen water masses are represented by O * 2 Öp, tÜ or Φ O 2 Öp, tÜ. However, these water masses can also be represented in geographic-space (latitude, longitude, and depth). In geographic-space, these water masses are represented by the spatial distribution of oxygen-percentiles, τÖx, y, z, tÜ. Φ O 2 Öp, tÜ contains information about the amount of oxygen present in the ocean, but not where it is located. Meanwhile, τÖx, y, z, tÜ contains information about where oxygen is distributed, but not how much there is. Future oxygen changes at a geographic location can therefore be expressed as the sum of a transformation effect which incorporates the changes in the oxygen-percentile relation in ventilation-space (O 2 * Ö p, tÜ or Φ O 2 Ö p, tÜ) for the water masses present at that location, and a redistribution effect which incorporates changes in the geographic-space distribution of oxygenpercentiles ÖτÖx, y, z, tÜÜ at that location (Figure 1; see Appendix B for details).

ESMs and Observations

In this paper, we use simulations from the Coupled Model Intercomparison Project Phase 6 (CMIP6) archive (Eyring et al., 2016;O'Neill et al., 2016) for 14 ESMs which provide dissolved oxygen, potential temperature, and salinity fields (see Table S1 in Supporting Information S1 for model details). Historical mean fields are computed from years 1950-2014 of historical simulations. Future projections are taken from the high-emission SSP5-8.5 and middle-of-the-road SSP2-4.5 scenario simulations. Dissolved oxygen, potential temperature, and salinity fields are available for SSP5-8.5 from 14 ESMs and for SSP2-4.5 from 13 ESMs. Ideal age fields are available for SSP5-8.5 from 10 ESMs and for SSP2-4.5 from 4 ESMs. The model drift for all fields are corrected using linear trends from pre-industrial control simulations (using the xMIP package; Busecke et al., 2023). For each model experiment, we use the average of up to five members when available (Table S1 in Supporting Information S1).

For comparison to observations, we use the oxygen climatology (objectively analyzed one-degree gridded field from all available years) from the World Ocean Atlas 2018 (Garcia et al., 2019). We also compare historical simulations to the interpolated time-resolved dissolved oxygen product presented by Ito (2021), based on the World Ocean Database. All simulations and data products were regridded to a uniform 1°⇥ 1°global horizontal grid and 100 vertical layers using the xesmf and xgcm packages (Abernathey et al., 2020;Zhuang et al., 2024). GSW-python was used to calculate oxygen saturation fields from potential temperature and salinity (Firing et al., 2021). The xhistogram package was used to help facilitate the transformation of ocean fields to oxygenpercentile coordinates (Abernathey et al., 2022).

Key Water Masses

We investigate the full spectrum of oxygen water masses in this paper, but we also highlight three key water masses for simplicity. We define "OMZs" using the 10th percentile (10% of ocean volume with lowest oxygen concentrations), "poorly ventilated waters" using the 50th percentile (50% of ocean volume with lowest oxygen concentrations), and "well-ventilated waters" as the reciprocal of poorly ventilated waters (50% of ocean volume with the highest oxygen concentrations). The 10th and 50th percentiles correspond to oxygen concentration thresholds of 82 and 177 μmol/kg, respectively, in observations from the World Ocean Atlas 2018 climatology (hereafter WOA18, see O 2 * in Figure 1b). These thresholds are similar to previous studies that used fixed oxygen thresholds to define water masses, such as OMZs where O 2 < 60 or 80 μmol/kg and low oxygenated waters where O 2 < 120 μmol/kg (e.g., Bopp et al., 2013;Busecke et al., 2022;Cabré et al., 2015;Cocco et al., 2013). However, unlike with fixed oxygen thresholds, the use of fixed percentile volumes provides a consistent way to track the evolution of water mass properties and content budgets over time.

In ventilation-space, the ensemble of 14 ESMs captures remarkably well the observed oxygen-percentile relation O 2

* ÖpÜ of WOA18 (Figure 1b), indicating that they reproduce the oxygen content of observed water masses despite showing mismatches in their location in geographic space (see details in Section 3.1 and Figure 2). Specifically, the ESM median reproduces the OMZs and the poorly ventilated waters bounded by the 10th and 50th oxygen percentiles within 10 μmol/kg of the WOA18 thresholds (see O 2 * in Figure 1b), as well as the oxygen content of OMZs (3% of ocean's oxygen), poorly ventilated (⇠1/3 of ocean's oxygen) and well-ventilated waters (2/3 of ocean's oxygen, see Φ O 2 in Figure 1b).

ESM Evaluation Using Oxygen Water Masses

We evaluate the skill of each ESM at simulating the historical geographic distribution of thermocline water masses (as defined by oxygen-percentiles) by computing a distance in ventilation-space, D, between simulations and observations for each water column (on a 1°by 1°grid):

Journal of Geophysical Research: Oceans 10.1029/2025JC022333

where f Öx, y, pÜ à PDF z âτÖx, y, zÜä is the probability density function of oxygen percentiles between 100 and 1,000 m that characterizes the water mass composition of the thermocline. Therefore, we can interpret the units of this metric as a water mass fraction distance. We compute this distance with N à 10 bins, each with a width of 10 percentiles. We ultimately rank ESMs by the global average of this distance,

yÜdx dy: Ö5Ü

Where A T is the total ocean area. Note that a perfect representation would yield a value of D à 0, whereas the

. The value of D for each ESM used in the study is listed in Table S2 in Supporting Information S1.

To quantify the precision of ESM ensemble projections in the tropical thermocline, we use two metrics: (a) the area between 30°S and 30°N over which at least 85% of ensemble members agree on the sign of oxygen trend and (b) the interquartile range (IQR) of the ensemble oxygen trends at each horizontal gridpoint averaged between 30°S and 30°N. For both metrics, oxygen trends are integrated over the thermocline between 100 and 1,000 m. An 85% agreement requires that 12 of 14 or 11 of 13 (for 7-ESM ensembles) or 23 of 27 (for 14-ESM ensemble) entries agree across SSP5-8.5 and SSP2-4.5 experiments. Journal of Geophysical Research: Oceans 10.1029/2025JC022333

Trends With Respect to Sea Surface Warming

In this paper, all temporal trends are presented as trends with respect to global sea surface warming. We compute time series of global mean sea surface temperature using the first regridded depth layer of potential temperature in SSP5-8.5 and SSP2-4.5 experiments. Trends with respect to surface warming are computed as the slope of a linear regression with global mean sea surface temperature and given in units of per °C.

Results

Challenges of Modeling Oxygen Water Mass Geography

Figure 2 not only illustrates the challenge of comparing observations and ESM simulations in geographic-space but also underlying consistencies that can be reconciled in ventilation-space. In geographic-space, the distribution of water masses defined by oxygen-percentiles ÖτÖx, y, zÜÜ exhibits similar patterns in observations from WOA18 and two example ESMs from the CMIP6 ensemble used here (GFDL-ESM4 and CanESM5-CanOE, Figure 2). Poorly ventilated waters (50th percentile) span all basins in the tropics at 300 m and extend into the abyssal ocean down to 5,000 m depth, and are separated from the well-ventilated water masses in subtropical and subpolar gyres by sharp horizontal and vertical oxygen gradients (see dashed black contours in Figure 2). Indeed, the general distribution of oxygen water masses is set to first order by ocean circulation patterns including gyres, subtropical cells and equatorial jets in the tropical thermocline (Brandt et al., 2010(Brandt et al., , 2015;;Calil, 2023;Duteil et al., 2014). Thus, the geography of these water masses is closely linked with the representation of these circulation features (see Figure S1 in Supporting Information S1 for circulation patterns in the tropical Pacific for GFDL-ESM4 and CanESM5-CanOE). OMZs (10th percentile) are found in ocean shadow zones where ocean circulation is weak (Luyten et al., 1983b), specifically in the tropical thermocline in each basin, and extending northward at depth in the Pacific Ocean (see solid black contours in Figure 2). By the end of the twenty-first century under highemissions forcing (SSP5-8.5 scenario), these broad features are reshaped and redistributed similarly in both ESMs (see red contours in Figure 2), in accordance with changing ocean circulation (Figure S1 in Supporting Information S1). This includes a systematic westward expansion of the Tropical Pacific OMZ (Busecke et al., 2022) and a slight poleward shift of poorly ventilated waters in the upper ocean as gyre circulations adjust to changing stratification and winds (Yang et al., 2020).

Despite these broad similarities in their historical and future distributions, the exact geographic shape of water mass boundaries is unique to each model and observations (Figure 2). For example, the Tropical Pacific OMZ extends across the basin at 300 m in WOA18, but only extends about half as far in GFDL-ESM4 and is even more confined to the eastern Pacific in 2c,and 2e). Similarly, there are differences in the vertical distribution and extent of the poorly ventilated waters in the South Pacific Ocean (Figures 2b, 2d, and 2f). These differences in oxygen distribution across models directly influence the pattern of their future perturbations. For instance, the region of oxygen loss in response to the Pacific OMZ westward expansion and poleward shift depends on the historical position of the OMZ and poorly ventilated water boundaries (Figures 2c-2f).

Oxygen Water Masses and Oxygen Changes in Historical Simulations

We evaluate the skill of ESMs in simulating geographically accurate ventilation pathways using a water mass distance metric, D, between the distribution of oxygen percentiles in the thermocline (100 and 1,000 m) simulated by each ESM and observations from WOA18 (globally averaged to D, see Methods 2.2.2 and Table S2 in Supporting Information S1 for ESMs skill). We compare the full 14-ESM ensemble (ensemble averaged D of 0.37) to the subset of seven ESMs with the highest skill (ensemble averaged D of 0.33) and the subset of seven ESMs with the lowest skill (ensemble averaged D of 0.40). Figures 3a-3c show the median spatial patterns of water mass distance in the three ESM ensembles. The 14 ESM ensemble median exhibits elevated water mass distances (i.e., poor representations) associated with subtropical gyres in each basin, as well as in the North Pacific and polar Southern Ocean (Figure 3a). The core regions of the tropical oxygen minimum zones in the Pacific and Indian Oceans, where a single water mass occupies the entire thermocline, exhibit the lowest water mass distances. When comparing the 7 ESMs with highest skill to the 7 ESMs with lowest skill, the bias in water mass distance decreases most dramatically in the North Pacific and polar Southern Ocean (Figures 3b and 3c), suggesting that the extent and depth of the North Pacific OMZ and the ventilation in the Southern Ocean play a prominent role in the differences between high and low skill models. Meanwhile, the tropical ocean exhibits Journal of Geophysical Research: Oceans 10.1029/2025JC022333 widespread, but generally modest, improvements associated with the positioning of the subtropical gyres and zonal extent of the tropical Pacific OMZ. Maps of water mass distance for individual ESMs are presented in Figure S2 in Supporting Information S1. Interestingly, related models (GFDL-CM4/ESM4, MPI-ESM1-2-LR/ HR, and NorESM2-LM/MM) tend to exhibit very similar patterns and magnitudes of water mass distance despite varying physical resolutions, suggesting that other model aspects (e.g., vertical coordinates, mixing parameterizations, and biological model) may be more important for setting the mean-state water mass representation. For instance, the model pair with the largest difference in water mass distance, CanESM5 and CanESM5-CanOE, differ only in their biogeochemistry.

We also evaluate the two groups of ESMs with respect to historical changes in oxygen content. Figure 3d compares the highest-skill and lowest-skill ensembles to the observed variability in thermocline (100-1,000 m) oxygen content between 1965 and 2014 (Ito, 2021). Note that internal variability (from seasonal to decadal timescales) is damped in the model ensembles from averaging over members and ESMs. The median thermocline deoxygenation is stronger in the highest-skill ensemble, but both groups of models agree well with observed trends in thermocline oxygen content and fall within the quantified uncertainty in the observations. Figure 3e shows the same comparison for changes in oxygen content integrated over the full ocean depth. Over the full ocean depth, the median deoxygenation is weaker in the highest-skill ensemble than the lowest-skill ensemble. Although the median timeseries for the lowest-skill ensemble shows closer agreement to observed variability in this case, both groups of models fall within the quantified uncertainty of the observations for a similar portion of the timeseries. Integrated over both the thermocline and full ocean depth, the highest-skill ensemble exhibits a consistently tighter ensemble spread (shading calculated by ensemble IQR) than the lowest-skill ensemble over the historical period.

Constraining Oxygen Projections

The geographic shifts in water mass boundaries found across ESMs influence the precision of the 14-ESM ensemble oxygen projections. In particular, the 14-ESM ensemble does not constrain whether oxygen will increase or decrease with warming over most of the tropical thermocline (Figure 4a; stippling where <85% of experiments agree on the sign of change) where oxygen projections are most disparate (Bopp et al., 2013;Busecke et al., 2022;Kwiatkowski et al., 2020). This contrasts with the mid-latitudes, where the 14-ESM ensemble robustly projects deoxygenation over almost all of the ocean thermocline, although the sign of change is not constrained in much of the North Atlantic and some limited areas of the Southern Ocean. Here, we demonstrate that selecting ESMs based on historical geographic distribution of thermocline oxygen water masses-used as a Journal of Geophysical Research: Oceans The broad spatial pattern of thermocline oxygen trends are surprisingly similar in the highest-skill, lowest-skill and 14-ESM ensembles, with deoxygenation projected at mid-and high-latitudes and weak changes or oxygenation projected in the tropics, but the precision of the projection varies widely across these ensembles (Figures 4a-4c). The tropical area with robust agreement on the sign of change increases drastically from 22% in the lowest-skill and 26% in the 14-ESMs ensembles to 58% in the highest-skill ensemble (Figure 4d). The largest gain in robustness is in the tropical Pacific Ocean, where the highest-skill ensemble projects more stable oxygen concentrations, and less widespread regions of oxygenation than the 14-ESM and lowest-skill ensembles (Figures 4a-4c). The improved precision is largely in the western tropical Pacific Ocean, rather than in the core regions of OMZs along the eastern boundary of the basin. The tropical Atlantic exhibits the opposite pattern, with increased precision in the poorly ventilated eastern boundary region. The model spread in oxygen trends in the tropics drops from 5.5 μmol/kg/°C in lowest-skill and 4.4 μmol/kg/°C in the 14-ESM ensembles to 3.0 μmol/kg/ °C in the highest-skill ensemble (Figure 4e). Outside of tropical regions, the highest-skill ensemble also suggests weaker deoxygenation than the lowest-skill ensemble in polar deep water formation regions, including the Weddell and Ross Seas in the Southern Ocean and the Labrador Sea in the North Atlantic Ocean.

To further demonstrate the relationship between model skill and precision of projections, we evaluate all possible 7-ESM subsets of the ensemble (3,432 combinations of 7 out of 14 ESMs, or 14-choose-7) using the two metrics: Journal of Geophysical Research: Oceans 10.1029/2025JC022333

area of sign agreement and ensemble spread. We find a strong relationship between the ensemble average skill with both the area of sign agreement and ensemble spread, with skill explaining 41% and 67% of the variance across ensembles, respectively (Figures 4d and 4e). This strongly suggests that the ability of ESMs to accurately represent ocean circulation and ventilation pathways-assessed here through oxygen patterns-is crucial for enhancing the precision of future projections. In the following sections, we leverage the oxygen water mass framework to contrast the response and dynamics of the highest-skill and lowest-skill ensembles and gain understanding on the fate of tropical oxygen.

Ventilation Changes Alleviate Deoxygenation in Skilled Models

The highest-and lowest-skill ensembles have distinct water mass properties in ventilation-space. The highestskill ensemble captures better the observed oxygen-percentile relation O 2

* ÖpÜ of WOA18 (Figure 5a), with steeply increasing oxygen thresholds at the highest and lowest percentile values, as well as OMZ and poorly ventilated water thresholds of 100 [78,199] μmol/kg and 186 [179,192] μmol/kg, comparable to the thresholds of 82 and 177 μmol/kg found in the WOA18 (model ensemble statistics are ensemble median and IQR in brackets). In contrast, the lowest-skill ensemble exhibits a more linear oxygen-percentile relation, and particularly fails to capture the observed relationship within the OMZs, underestimating the observed OMZ threshold by about half at 41 [38,48] μmol/kg. As a result, the highest-skill ensemble better reproduces the observed cumulative oxygen content of water masses, Φ O 2 Öp, tÜ, with the observed Φ O 2 Öp, tÜ falling within the IQR of the highest-skill ensemble for nearly all percentile values, but falling outside of the lowest-skill ensemble IQR for all values of p less than 70%.

In this ventilation-space, the CMIP6 ESMs largely agree on how water mass deoxygenation will be distributed, with well-ventilated waters deoxygenating much faster than poorly ventilated waters and OMZs in all models; however, the highest-skill ensemble projects systematically weaker deoxygenation than the lowest-skill ensemble (Figure 6a). Over the global ocean (100th percentile), the highest-skill ensemble projects a trend of 2.9 [ 3.7, 2.1] Pmol O 2 /°C across both middle of the road (SSP2-4.5) and high emission (SSP5-8.5) experiments between 2015 and 2100, whereas the lowest-skill ensemble projects a trend of 4.2 [ 4.8, 3.0] Pmol O 2 /°C. For both groups of ESMs, about three quarters of this deoxygenation occurs in well-ventilated waters (81% [67%, 82%] in the highest-skill ensemble and 77% [66%, 81%] in the lowest-skill ensemble). We note that the difference in global deoxygenation rates between the highest and lowest-skill ensembles is derived primarily from deep (below 2,000 m) well-ventilated waters, whereas waters above 1,000 m deoxygenate at similar rates for both ensembles (Figure S3 in Supporting Information S1). Averaged over all OMZs, the oxygen content remains quite stable (IQR includes zero change) in both groups of ESMs, with a weak oxygenation of á0.05 [ 0.12, á0.11] Pmol/°C in the highest-skill ensemble and á0.03 [ 0.06, á0.06] Pmol/°C in the lowest-skill ensemble.

ESMs also generally agree on the processes driving these future changes in ventilation-space. Water mass transformations are driven by changes in the air-sea flux of oxygen ÖFÜ, the respiratory sink of oxygen in the ocean interior (S, regulated by respiration rates and water mass age), and the transport of oxygen via mixing across water masses (M, see Appendix A and Figure 1). To gain insights into the relative contributions of these processes, we decompose changes in cumulative oxygen content ÖΦ O 2 Ü using the canonical separation into (a) changes in oxygen saturation (Φ O 2sat ; Figure 6b) which approximates changes in oxygen due to air-sea flux at the Journal of Geophysical Research: Oceans Integrated over the global ocean (i.e., contributions from mixing are zero; M à M O 2sat à M AOU à 0), changes in oxygen saturation contribute comparably to deoxygenation in the median of the highest-skill and lowest-skill ESM ensembles (Figure 6b). This indicates that the total loss of oxygen by efflux at the air-sea interface, F, with warming is comparable in the highest-and lowest-skill ensembles (Φ O 2sat changes at a rate of 1.5 [ 1.7, 1.3] Pmol O 2 /°C vs. 1.5 [ 2.0, 1.4] Pmol O 2 /°C). In contrast, globally integrated changes in AOU are weaker in the highest-skill ensemble than in the lowest-skill ensemble (Φ AOU changes at a rate of 1.1 [ 2.0, 0.8] Pmol O 2 /°C vs. 2.1 [ 2.7, 1.7] Pmol O 2 /°C), indicating a weaker increase in the total respiratory sink (S) and explaining the weaker global ocean deoxygenation in the highest skill ensemble (Figures 6b and 6d). As noted for ensemble differences in the total deoxygenation trend, differences in AOU trends between the highest-skill and lowest-skill ensembles primarily come from the deep ocean (Figure S3 in Supporting Information S1). However, AOU trends in the model pair that differs only in their representation of biogeochemistry (CanESM5 and CanESM5-CanOE, both in the highest-skill ensemble) diverge primarily (75% of difference) in the upper 2,000 m (Figure S3 in Supporting Information S1), suggesting that biogeochemical representation has much stronger influence in the upper ocean compared to the deep ocean. Indeed, we show below that globally integrated increases in the respiratory sink S are likely controlled by changes in water mass age, rather than changes in local respiration rates. The globally integrated estimate of S assumes that surface waters are subducted at saturation, but deviations from this assumption may also influence the estimate (Frölicher et al., 2020;Koeve & Kähler, 2016).

The balance between changes in oxygen saturation ÖΦ O 2sat Ü and AOU ÖΦ AOU Ü varies across water masses (Figures 6b and 6c). Φ O 2sat trends yield similar rates of deoxygenation at nearly all percentiles in both ensembles. In contrast, Φ AOU trends switch from oxygenation at low oxygen-percentiles to deoxygenation at high oxygen- Journal of Geophysical Research: Oceans 10.1029/2025JC022333

percentiles (switch from positive slopes of Φ AOU trends to negative slopes), which homogenizes the oxygen distribution across water masses by supplying oxygen from well-ventilated waters to poorly ventilated waters and OMZs (Figure 6c). As a result of this balance, three regimes emerge in the key water masses. In the OMZs, the deoxygenation by Φ O 2sat changes and oxygenation by Φ AOU changes compensate to maintain a stable oxygen content (see Figure 6b at 10th percentile) in both groups of ESMs, consistent with previous work which has noted the geographic pattern of compensating changes in O 2sat and AOU in tropical OMZs with warming (Bopp et al., 2017;Busecke et al., 2022;Kwiatkowski et al., 2020). In poorly ventilation waters, the deoxygenation by Φ O 2sat changes exceeds the oxygenation by Φ AOU changes in both ensembles, but the Φ AOU -driven oxygenation is higher in the highest-skill ensemble resulting in lower net deoxygenation (Figures 6a-6c). Finally, in wellventilated waters, both Φ O 2sat and Φ AOU changes drive deoxygenation but with again a weaker Φ AOU -driven deoxygenation simulated in the highest-skill ensemble (Figures 6a-6c).

We can qualitatively infer the contributions of F, S and M to water mass deoxygenation. Since we cannot distinguish the contributions of S and M AOU to Φ AOU changes for the non-globally integrated case, we also examine changes in water mass ideal age (the average time since a water mass has been exposed to the surface, available for 14 out of 23 experiments; see Table S1 in Supporting Information S1). The ocean as a whole gets older with warming in all available experiments (Figure 6d). Since global rates of organic matter export are expected to decline with warming (Bopp et al., 2013;Kwiatkowski et al., 2020), water mass age is likely the primary control driving increases in the global respiratory sink. In contrast, OMZs and poorly ventilated waters get younger. OMZs and poorly ventilated waters do not outcrop at the surface (Figure S4 in Supporting Information S1), meaning that these water masses are not directly subject to changes in air-sea fluxes of oxygen, F, nor the resetting of water mass age at the surface. Thus, trends in Φ O 2sat and ideal age in poorly ventilated waters and OMZs must be the effect of mixing. We can therefore infer that poorly ventilated waters and OMZs experience increased mixing with warmer and younger (i.e., shallower) waters relative to older and cooler (i.e., deeper) waters, as has been hypothesized in previous studies (Bryan et al., 2006;Gnanadesikan et al., 2007Gnanadesikan et al., , 2012)). Wellventilated waters must experience the reciprocal effect from mixing, including a supply of cool and old waters which partially offsets the direct effects of changes in F and likely exacerbates deoxygenation driven by changes in Φ AOU and S. We note that, although changes in local respiration rates may also play a role in setting S, these results demonstrate that changing ventilation dynamics can explain the pattern of Φ AOU trends to first-order.

Examining the cumulative mean age trends in both groups of ESMs, we see that the changes in OMZ ventilation appear to be stronger in the highest-skill ensemble, with OMZs getting between 10 and 20 years younger per degree surface warming, compared to about 5 years younger in the lowest-skill ensemble. This suggests that the OMZs in the highest-skill ensemble experience a greater increase in the contribution of younger upper ocean oxygenated waters relative to older less oxygenated waters than in the lowest skill ensemble. Trends in Φ O 2sat support this larger increase in the contribution of younger upper ocean waters in the highest-skill ensemble. Indeed, there is a subtle, but systematic, stronger thermal gain in the highest-skill ensemble across OMZs and poorly ventilated waters, whereas there is no systematic difference in the global thermal signal. Stronger shifts in the mixing of warmer and younger waters into OMZs likely contributes to the wider spread of OMZ oxygen projections in the highest-skill ensemble compared to the lowest-skill ensemble in these low oxygen regions.

Constraints on Water Mass Redistribution in Skilled Models Improve Precision of Oxygen Projections

The geographic redistribution of water masses can lead to strong local oxygen trends that reinforce the changes associated with transformation within water masses. We show that this phenomenon is particularly important for understanding oxygen trends in the tropical thermocline, a region which is dominated by water masses that exhibit weak transformation trends. Figure 7 shows the pattern of OMZ redistribution in response to warming using trends in water mass contribution (i.e., the fraction of integrated waters occupied by OMZs) by depth and latitude. By construction, the highest-skill ensemble captures better the observed distribution of historical water masses in depth and latitude than the lowest-skill ensemble. Specifically it reproduces better the observed extent of OMZs in the deep tropical ocean (too deep in the lowest-skill ensemble) and the presence of OMZs in the northern extratropics (>20% water mass contribution in the highest-skill ensemble vs. <10% in the lowest-skill ensemble; Figures 7a and 7c). It is worth noting that the lowest-skill ensemble better represents OMZ contributions in the northern tropics (Figure 7c). ESMs tend to project a decline of OMZ volume in the tropical Journal of Geophysical Research: Oceans 10.1029/2025JC022333 thermocline (trends less than zero between 100 and 1,000 m), but this effect is weaker and more tightly constrained in the highest-skill ensemble (Figure 7b). Integrated over the tropical thermocline, the highest-skill ensemble projects a shift in OMZ contribution of 0.3 [ 0.6, á0.0]%, whereas the lowest-skill ensemble projects a shift of 1.0 [ 1.3, á0.3]%. Integrated over all depths, the two groups of ESMs project a similar amount of OMZ volume that gets redistributed from the tropics to northern mid-latitudes (negative trends in tropics and positive trends in northern hemisphere, Figure 7d), but more of this exported volume originates below the thermocline in the highest-skill ensemble than in the lowest-skill ensemble (Figure 7b). The highest-skill ensemble also exhibits tighter constraints on this redistribution of OMZ volume from the tropics than the lowest-skill ensemble, particularly between 15°S and 15°N where the ensemble IQRs differ by nearly a factor of two. These results suggest that more accurately simulated ventilation pathways in the tropical thermocline may be more geographically stable with warming than comparatively unrealistic ventilation pathways.

We quantify the relative influence of water mass transformation and geographic redistribution on the magnitudes and uncertainties of zonal mean thermocline oxygen projections in the Pacific, Atlantic, and Indian Ocean basins (Figure 8). At first glance, the meridional patterns of transformation and redistribution in the thermocline appear similar in the highest-and lowest-skill ensembles. In both ensembles, redistribution generally controls the deoxygenation in the northern extra-tropics and the oxygenation in the tropics, whereas transformation controls the deoxygenation in the southern extra-tropics. Yet, the effects of redistribution dominate the spread across ESM projections at all latitudes and in all basins (violet shading), whereas the effects of transformation are relatively well constrained across ESMs (orange shading in Figure 8). We can therefore attribute the uncertain sign of oxygen trends in the tropical thermocline, reported in previous studies (Bopp et al., 2013;Kwiatkowski et al., 2020), primarily to redistribution effects. This uncertainty associated with water mass redistribution is, however, significantly mitigated by selecting the highest-skill ensemble. Compared to the full 14-ESM ensemble, the IQR of total oxygen projections in the highest-skill ensemble is tighter by an average of 47%, 32%, and 28% at tropical latitudes in the Pacific, Atlantic, and Indian Oceans, respectively. Similarly, the IQRs tighten by 29% and 38% on average in the north Pacific and Atlantic Oceans (30°N-60°N). As a result, the agreement on the sign of redistribution and hence of total oxygen changes are improved, especially in the tropical thermocline. The highest-skill ensemble projects a weak deoxygenation of the tropical Pacific Ocean, and an oxygenation of the tropical Atlantic and Indian Oceans, whereas the lowest-skill ensemble projects a weak oxygenation with a high spread that renders the sign of change ambiguous in all tropical oceans (Figures 8b, 8d, and 8f). Journal of Geophysical Research: Oceans 10.1029/2025JC022333

Discussion and Conclusions

We provide a framework for evaluating ESM projections and a new perspective on the processes regulating ocean oxygen in a changing climate by characterizing water masses based on ventilation patterns. We isolate the imprints of transformation by air-sea fluxes, respiration and mixing, as well as redistribution of these water masses. We find that ensembles of models with different mean-state representations of water masses yield distinct projections. ESMs with better thermocline water mass representations project a global deoxygenation that is 30% lower than in the less skilled models, mostly from weaker deoxygenation of the deep ocean. The global efflux of oxygen from the direct effect of warming (i.e., temperature driven changes in oxygen solubility) is consistent across all ESMs, but the highest-skill group exhibits weaker increases in the overall respiratory sink, likely controlled by physical ventilation. Models with highest skill also show a weaker deoxygenation of poorly ventilated waters, in part due to stronger changes in mixing ratios of source waters which oxygenate the OMZs. This result is in agreement with prior work which suggests the mixing ratio between thermocline and deeper water masses as the primary control on declining OMZ water age in ESM projections (Gnanadesikan et al., 2007(Gnanadesikan et al., , 2012)). The supply of oxygen to OMZs in the tropics is tied to both thermocline circulation in the upper ocean and thermohaline circulation in the deep ocean (Bopp et al., 2017;Gnanadesikan et al., 2007;Khatiwala et al., 2012;Oschlies et al., 2018;Wyrtki, 1962), so a holistic understanding of how global ocean dynamics respond to warming is important for constraining oxygen projections at a regional scale (Cabré et al., 2015;Marshall & Speer, 2012;Naveira Garabato et al., 2017). The results of this paper highlight the global nature of ventilation pathways and the importance of interactions between the thermocline and deep ocean. By grouping models based on their representation of thermocline water masses, we isolate distinct patterns of deoxygenation in deep water formation regions, which likely explain the diverging projections for deep ocean deoxygenation in the two model groups.

The contrast between models in the relative influence of solubility and ventilation in controlling future deoxygenation echoes a progression in climate model projections over the past decades. Early coarse-resolution simulations estimated that AOU changes accounted for up to 75% of projected global deoxygenation (Bopp et al., 2002;Matear & Hirst, 2003), whereas more recent generations have estimated this value closer to 50% (Bopp et al., 2017;Frölicher et al., 2009). The higher-skilled CMIP6 models analyzed in this study also estimate about a 50% contribution from AOU changes, whereas the lower-skilled models simulate a 10% stronger contribution from AOU.

Water mass transformation sets the global levels of deoxygenation (by definition), but water mass redistribution is the main driver of heterogeneities and model ensemble spread at regional scales. Redistribution, and therefore projection precision, are better constrained in the highest-skill ensemble of ESMs than in the lowest-skill ensemble. We find that, in addition to simulating a stronger oxygenation of OMZs, models with highest skill Journal of Geophysical Research: Oceans 10.1029/2025JC022333

also simulate more stable OMZ geography (i.e., less redistribution of OMZ volume). This is consistent with the results of Busecke et al. (2019) who showed using high-resolution models that the shape of the OMZ was far less sensitive to climate variability in models that better represented the circulation pathways supplying oxygen to the OMZ (e.g., equatorial current system). This stronger "stability," or resistance to change under perturbation in highest-skill models ultimately leads to weaker and more constrained oxygen projections over much of the tropical thermocline. Specifically, this stability constrains regions in the central and western tropical Pacific that are influenced by changes in the zonal extent of OMZ waters, rather than regions in the eastern tropical Pacific that are likely to always be dominated by OMZs. This reduces the uncertainties that dominated regional oxygen projections in previous model inter-comparison studies (Bopp et al., 2013;Cabré et al., 2015;Cocco et al., 2013;Kwiatkowski et al., 2020).

Projecting changes in aerobic conditions at fixed geographic locations is important for assessing risks associated with ecosystem responses and ecosystem services (Breitburg, 2002). By tightening constraints on water mass redistribution in the tropical thermocline by more than 30% (and more than 40% in the Pacific), we disambiguate the sign of thermocline oxygen change over about half of the tropical Pacific and large areas of the tropical Atlantic and Indian Oceans. This is critical to inform future management of marine resources, including important tuna fisheries in the western tropical Pacific Ocean (Drakou et al., 2018), and similarly in the western Indian Ocean (Kaplan et al., 2014). Specifically, this work contributes to understanding and constraining the oxygen projections used to assess higher-trophic ecosystem response and inform industrial and conservation efforts (Heneghan et al., 2021;Lefort et al., 2015;Tittensor et al., 2018). As such, conceptualizing the ocean as being composed of water masses defined by oxygen-percentiles is useful for both developing a fundamental understanding of the ocean's response to warming, as well as strengthening the practical application of ESM projections.

The oxygen water mass framework introduced here successfully leverages a ventilation-space to filter spatial biases across simulations and elucidate consistent changes in ocean oxygen dynamics under climate change. However, this approach is still subject to limitations. For example, models consistently simulate the Atlantic Ocean OMZ as suboxic, contrary to observations (Séférian et al., 2020), likely due to anomalously weak tropical circulation pathways in coarse resolution models (Brandt et al., 2015;Calil, 2023;Fu et al., 2022). Thus, we take caution in interpreting the results of this analysis in the tropical Atlantic Ocean. Furthermore, water masses defined by oxygen percentiles group geographically distinct water masses that may be subject to disparate local effects. For example, although OMZs in the tropical Pacific and Atlantic Oceans are governed by eastern boundary upwelling dynamics (Bograd et al., 2023), future projections of the Indian Ocean OMZ are governed by unique local circulation features including the Indonesian Throughflow and Persian Gulf outflow (Ditkovsky et al., 2023;Lachkar et al., 2023). Oxygen-percentiles also group together northern and southern waters, yet the Southern Ocean dynamics may respond differently to surface warming than dynamics in the northern basins where continents play a key role (Doddridge et al., 2019;Sallée, 2018). The ventilation-space framework reduces the ocean to a single dimension by aggregating similarly ventilated waters in order to elucidate patterns that are obscured in geographic space. But consequently, it cannot fully capture the range of regional responses, and should not be a substitute for detailed analyses within a specific region of interest.

The approach of this paper is informed by water mass transformation theory (Groeskamp et al., 2019;Walin, 1982), and in particular, recent developments in tracer-percentile frameworks (Holmes et al., 2019;Sohail et al., 2021Sohail et al., , 2022Sohail et al., , 2023) ) and biogeochemical applications (Iudicone et al., 2011). Although almost all studies in water mass transformation have focused on conservative tracer coordinates, this study pioneers the use of a nonconservative tracer coordinate, which is useful for characterizing ventilation and oxygen changes. Compared to traditionally used water mass coordinates like temperature or density, which largely follow latitude and depth, oxygen or oxygen-percentiles capture a more complex spatial partitioning of the ocean-one based on connectivity with the surface or ventilation. This partitioning is prudent because the ventilation of high-oxygen water masses and OMZs likely have different responses to surface forcing (Busecke et al., 2022;Gnanadesikan et al., 2007;Luyten et al., 1983a). The efficient supply pathways which maintain high-oxygen water masses weaken as stratification intensifies (Couespel et al., 2019;Oschlies et al., 2018), while OMZ ventilation is sensitive to the strength of turbulent mixing in the ocean (Bahl et al., 2019;Duteil & Oschlies, 2011;Ito et al., 2022;Lévy et al., 2022) and the shape of the OMZ itself (Busecke et al., 2019). Furthermore, the responses of ventilation pathways are not only relevant for studying oxygen dynamics but are also responsible for transporting and storing heat, freshwater and carbon in the ocean (Bronselaer & Zanna, 2020;Mannis et al., 2024; Journal of Geophysical Research: Oceans 10.1029/2025JC022333 âO 2 Öx, y, z, tÜä à τÖx, y, z, tÜ: ÖB1Ü Future oxygen changes in geographic-space can therefore be expressed as follows: d dt O 2 Öx, y, z, tÜ ≈ d dt F pÖO * 2 ,tÜ τO 2 âτÖx, y, z, t à t ref Üä |ÇÇÇÇÇÇÇÇÇÇÇÇÇÇÇÇÇ Ç{zÇÇÇÇÇÇÇÇÇÇÇÇÇÇÇÇÇ Ç} transformation effect á d dt F pÖO * 2 ,tàt ref Ü τO 2 âτÖx, y, z, tÜä |ÇÇÇÇÇÇÇÇÇÇÇÇÇÇÇ{zÇÇÇÇÇÇÇÇÇÇÇÇÇÇÇ} redistribution effect , ÖB2Ü where F pÖO * 2 ,tÜ τO 2 is the inverse of the transfer function in Equation B1, F pÖO * 2 ,tÜ O 2 τ , and t ref is a reference time (the mean over 1950-2014). F pÖO * 2 ,tÜ τO 2 maps from τÖx, y, z, tÜ to O 2 Öx, y, z, tÜ using p O * 2 , t . The transformation effect term in Equation B2 allows the oxygen-percentile relation to evolve in time ÖpÖO 2

* , tÜÜ, but holds τÖx, y, z, t à t ref Ü fixed at its reference state. So, this term captures the evolution of oxygen in a geographic region if the local water masses changed only according to their evolution in ventilation-space, which is aggregated globally (Figure 1a). Conversely, the redistribution effect term holds the oxygen-percentile relation fixed ÖpÖO 2 * , t à t ref ÜÜ, but allows τÖx, y, z, tÜ to evolve in time. Thus, this term captures the evolution of oxygen in a geographic region if local water masses experienced no evolution in ventilation-space (Figure 1b). There is a small residual associated with nonlinearities between the transformation and redistribution terms, but this effect is negligible in our results.

Water masses defined by oxygen-percentiles are spatially discontiguous, so the signal of water mass redistribution includes both a local reshaping of features by changes in advective pathways and remote adjustments to the geographic distribution of oxygen-percentiles. The decomposition outlined above amounts to a similar conception of redistributed properties as the "minimum transformation method" developed by Zika et al. (2021) to describe trends in ocean heat content. Although the methods for computing redistribution differ between that study and this work, redistribution of water masses here may be interpreted similarly as in Zika et al. (2021) with Journal of Geophysical Research: Oceans 10.1029/2025JC022333 comments and suggestions in the review process. The authors acknowledge the World Climate Research Program, which, through its Working Group on Coupled Modeling, coordinated and promoted CMIP6, participating climate modeling groups for producing and publishing their model output, the Earth System Grid Federation (ESGF) for archiving the data and providing access, and the multiple funding agencies who support CMIP6 and ESGF. Journal of Geophysical Research: Oceans 10.1029/2025JC022333