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

# Irreversible loss in marine ecosystem habitability after a temperature overshoot

Yeray Santana-Falcón ID <fnref n="1" />✉, Akitomo Yamamoto ID ², Andrew Lenton ID <fnref n="3" />, Chris D. Jones ID ⁴, Friedrich A. Burger ⁵,⁶, Jasmin G. John ID ⁷,⁸, Jerry Tjiputra ID ⁹, Jörg Schwinger ⁹, Michio Kawamiya ID ², Thomas L. Frölicher ID ⁵,⁶, Tilo Ziehn ID ¹⁰ & Roland Séférian ID <fnref n="1" />

Anthropogenic warming of the oceans and associated deoxygenation are altering marine ecosystems. Current knowledge suggests these changes may be reversible on a centennial timescale at the ocean surface but irreversible at deeper depths even if global warming were to ameliorate. In contrast, the marine ecosystem's response to these persistent changes remains poorly elucidated. Here we explore to what extent global warming may drive alterations in marine habitats by exploring the evolution of a metabolic index that captures marine organisms' ecophysiological response to both temperature and oxygen changes, throughout an idealised ramp-up/ramp-down atmospheric carbon dioxide concentration and an overshoot scenarios. Using a multi-model approach; we find that changes in ocean temperature and oxygen drive a centuries-long irreversible loss in the habitable volume of the upper 1000 m of the world ocean. These results suggest that the combined effect of warming and deoxygenation will have profound and long-lasting impacts on the viability of marine ecosystems, well after global temperatures have peaked.

<footnote n="1">¹CNRM, Université de Toulouse, Météo-France, CNRS, Toulouse, France. ² Japan Agency for Marine-Earth Science and Technology, Yokohama, Japan.</footnote>

<footnote n="3">³Permanent Carbon Locking Future Science Platform, CSIRO, Hobart, TAS, Australia. ⁴Met Office Hadley Centre, Exeter EX1 3PB, UK. ⁵Climate and Environmental Physics, Physics Institute, University of Bern, Bern, Switzerland. ⁶Oeschger Centre for Climate Change Research, University of Bern, Bern, Switzerland. ⁷NOAA/OAR/Atlantic Oceanographic and Meteorological Laboratory, Miami, FL, USA. ⁸NOAA/OAR/Geophysical Fluid Dynamics Laboratory, Princeton, NJ, USA. ⁹NORCE Climate & Environment, Bjerknes Centre for Climate Research, Bergen, Norway. ¹⁰CSIRO Environment, Aspendale, VIC, Australia. ✉email: yeray.santana@meteo.fr</footnote>
Recent trends in global emissions suggest that the cumulative emissions budget compatible with limiting global warming to 1.5 °C, as agreed at Paris COP21 (2015), and reaffirmed at Glasgow COP26 (2021), is likely to be exceeded with corresponding additional warming¹⁻³. Hence, the latest IPCC reports⁴,⁵ include emission pathways for which reductions in greenhouse gas emissions are complemented by technologies that allow direct removal of excess atmospheric carbon dioxide (CO₂)⁶⁻⁹. These carbon dioxide removal (CDR) technologies¹⁰ are anticipated to be available in following decades¹¹,¹². They are needed to counterbalance hard-to-abate emissions in the mid-term and achieve and sustain net negative emissions in the long-term by exceeding annual residual gross emissions⁴,¹³,¹⁴. Therefore, it is essential to understand the response of the Earth system to reduce emissions either artificially or by enhancing the sink strength of the ocean or land biosphere¹⁵.

The Carbon Dioxide Removal Model Intercomparison Project¹⁶ (CDRMIP) provides an idealised 140 year long 1% year⁻¹ CO₂ ramp-down simulation and constant CO₂ thereafter from the end of 140 year long 1% year⁻¹ CO₂ ramp-up simulation. These simulations have allowed the hysteresis of relevant Earth system parameters to be quantified as the difference between their ramp-up and ramp-down paths (for example ref. ¹⁷⁻²²). Hysteresis refers here to the dependence of the Earth system on its past state, not only on its current state. Ramp-up/ramp-down simulations are useful as they provide an idealised framework consistent with the definition of hysteresis, i.e., CO₂ concentrations begin and finish at the same steady-state. However, an experiment with more realistic levels of overshoot is needed to assess the implications of CDR in a plausible socio-economic scenario, for example, the shared socio-economic scenario SSP5-3.4-OS²³⁻²⁵, which has also been used to detect hysteresis in the Earth systems. SSP5-3.4-OS follows a high-emission pathway consistent with the no-mitigation scenario SSP5-8.5 until year 2040, at which point aggressive mitigation is assumed to quickly decrease CO₂ concentrations, resulting in net zero emission levels by year 2080. Then, negative emissions are phased out by 2170. In the ocean, reversibility studies have mainly focused on temperature, dissolved oxygen (DO), and acidity since ocean warming, deoxygenation, and acidification constitute notable risks to life in the world ocean²⁶⁻²⁹. However, the impacts of changes in these potential ocean ecosystem stressors after a temporary overshoot, referred here as a consecutive increase and decrease of CO₂ levels, on the resilience of marine ecosystems have been analysed in isolation so far (for example refs. ¹⁹,³⁰), while their combined effects remain a gap in our knowledge.

Habitability of marine ecosystems depends on the rate of oxygen (O₂) supply that must exceed the resting metabolic demand. As the ocean warms in response to global warming³¹, organisms' metabolic demands are anticipated to increase³²,³³. At the same time, ocean warming reduces the solubility of oxygen in the seawater and increases vertical stratification³⁴, reducing DO concentrations altogether²⁶,³⁵,³⁶. To quantitatively account for the combined effects of warming and oxygen loss, the metabolic index³⁷⁻⁴⁰ (Φ) is defined as the oxygen supply (S) to demand (D) ratio for a marine organism. Oxygen supply increases with ambient O₂ pressure (pO₂), which scales with temperature (T), and with respiratory efficacy, such that S = αₛ(T)B⁶pO₂. Respiratory efficacy results from a per-mass rate of gas transfer between the water and the animal (αₛ), and it scales with body size (B⁶). Resting metabolic demand also scales with body mass (B⁶) and with absolute temperature, such that D = αₛB⁶exp{−E₉/kB(1/T−1/Tᵣₑ)}, where αₛ is a taxon-specific baseline metabolic rate, and kB is Boltzman's constant. The parameters Es and Ed (eV) represent the temperature dependence of the O₂ supply and

the metabolic rate, respectively. Therefore, Φ can be expressed as

$$ \Phi = \frac{S}{D} = \frac{\alpha_S}{\alpha_D} B^e pO_2 \exp \left( \frac{E_d - E_s}{k_B} \left( \frac{1}{T} - \frac{1}{T_{ref}} \right) \right) $$

Using estimates of this key metabolic index computed for 72 different species across the water column (see Methods), we determine the level of hysteresis of marine ecosystems' habitability after an overshoot. For this, we use a suite of seven Earth system model (ESM) simulations of the fifth⁴¹ (CMIP5) and the sixth⁴² (CMIP6) Intercomparison Projects (see Supplementary Note 1 and Supplementary Table 1) (GFDL-ESM4;⁴³ GFDL-ESM2M;⁴⁴ UKESM;⁴⁵ MIROC-ES2L;⁴⁶ CNRM-ESM2-1;⁴⁷ ACCESS-ESM1;⁴⁸ NorESM2-LM⁴⁹,⁵⁰) following the 1pctCO₂ and (four of the ESM simulations following) the SSP5-3.4-OS simulations to explore the response of Φ to a consecutive increase and decrease in atmospheric CO₂ concentrations. The metabolic index has been successfully deployed in biogeography studies³⁸, climate change impacts³⁷,⁵¹, assessments of past mass extinction using fossil records³⁹,⁴⁰, and projecting risks of future mass extinctions of marine organisms⁵². Here, we aim to use it to estimate the loss of habitable ocean volume after levels of CO₂ recover from overshoot.

## Results

**Global ocean response.** Multi-model global-averaged sea surface temperature (SST) responds almost immediately to a 1% increase per year of atmospheric CO₂ (blue line in Fig. 1 and Supplementary Fig. 1). In the 1pctCO₂ experiment, peaking at 4x CO₂ (~1140 ppm) occurs at year 140 of the experiment. However, peak SST occurs after a seven-year lag, reaching a maximum value of ~3.1 °C above the pre-industrial level. When the pre-industrial CO₂ levels are restored at model year 280, SST is still warmer (+0.7 °C) than the pre-industrial value. From model year 280 until the end of the simulation (model year 340), CO₂ concentrations are kept constant at ~285 ppm. SST continues to decrease during this phase, but it is still 0.54 °C warmer than the pre-industrial initial state at the end of the experiment. In the SSP5-3.4-OS, peak CO₂ (~571 ppm) occurs at year 2062, reaching net zero emissions after the overshoot around year 320 from the beginning of the scenario. SST follows the same pattern, though with a slightly longer delay of 10 years from peak CO₂. As in the 1pctCO₂ simulation, SST continues to decrease after reaching net zero, being 0.09 °C warmer at the end of the scenario with respect to year 2170.

As SST increases, surface waters' capability to hold oxygen diminishes²⁶,⁵³,⁵⁴. The sea surface pO₂ time series (orange line in Fig. 1 and Supplementary Fig. 1) mimics that of SST, reaching a minimum after a similar lag from the CO₂ peak. A similar deviation from initial conditions is observed after CO₂ has returned to pre-industrial levels in the case of the 1pctCO₂ simulation, and reaches net zero emissions in the case of SSP5-3.4-OS. In the idealised scenario, this negative hysteresis is represented by a decrease of surface pO₂ at the end of the experiment of about -1.1% with respect to initial values. In the SSP5-3.4-OS, the hysteresis is represented by a decrease of surface pO₂ of about -0.24% at the end of the experiment with respect to year 2170.

For temperature, the delay with respect to the CO₂ peak year increases in the epipelagic (0–200 m depth; Fig. 1 middle panels) and mesopelagic (200–1000 m depth; Fig. 1 bottom panels) waters for both experiments. In the 1pctCO₂ simulation, the temperature peak in epipelagic waters roughly coincides with the peak year for SST but occurs much later (57 years) in mesopelagic waters, indicating a strong vertical gradient in hysteresis³⁰. With increasing depth, some decoupling of pO₂ from temperature
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Time series of atmospheric CO₂, ocean temperature, oxygen partial pressure, and metabolic index for the global ocean over the ramp-up and ramp-down and the SSP5-3.4-OS experiments. Thick lines illustrate multi-model mean. Shading indicates multi-model standard deviation. Vertical dashed lines indicate the time at which atmospheric carbon dioxide peaks (model year 140 and year 2062, for each experiment respectively), and the year it returns to initial values for the 1pctCO₂ experiment (model year 280) and net zero is reached after the overshoot (year 2170) for the SSP5-3.4-OS experiment. Time series are shown for the sea surface, and averaged over epipelagic (0–200 m) and mesopelagic (200–1000 m) waters (top, middle and bottom panels respectively). Metabolic Index corresponds to the mean over 72 species (see Methods).</figcaption></figure>

occurs, as it becomes evident with a slight delay in the minimum of pO₂ with respect to peak temperature in epipelagic waters (1 year after temperature peaks), and an earlier minimum of pO₂ in mesopelagic waters (~45 years before temperature peaks). Furthermore, pO₂ presents slightly positive hysteresis (~1.0% higher than at the beginning of the experiment) at this depth layer consistent with the increase of oxygen found in deep waters of the North Atlantic by Bertini and Tjiputra²¹. In the overshoot scenario, the temperature in mesopelagic waters continues to increase until the end of the experiment, at which point the global ocean is 0.1 °C warmer than the year 2170 when negative emissions are phased out. The minimum in pO₂ also occurs at the end of the scenario, being slightly lower than at the moment net zero emissions are reached.

The Metabolic Index Φ (green lines in Fig. 1 and Supplementary Fig. 1) is positively correlated with pO₂ at the sea surface (hence, negatively correlated with SST) throughout the full experiments. However, Φ hysteresis at the end of the 1pctCO₂ simulation (−5.7% with respect to the initial value) almost doubles that for SST (3.0% with respect to the initial value) and triple that for pO₂ (−1.1% with respect to the initial value) in absolute values. In epipelagic waters (see Supplementary Fig. 2), the level of hysteresis for Φ (−5.8% with respect to initial value) remains similar to the surface, with the trough occurring later (13 years) after the peak of temperature and the pO₂ minimum. In

mesopelagic waters, Φ's level of hysteresis slightly increases (−7.3% with respect to the initial value). The evolution of Φ at this level differs from the changes in the other two variables; Φ minimum occurs 14 years after the pO₂ minimum but 34 years before the peak of temperature, thus reflecting a legacy of concomitant changes in temperature and pO₂ on ocean habitability. In the overshoot scenario, Φ also mimics the pattern displayed by pO₂. In agreement with the higher hysteresis observed in the 1pctCO₂ simulation, Φ is 2.05 and 2.21% lower at the end of the experiment with respect to year 2170 in the sea surface and the epipelagic waters respectively. In the mesopelagic, Φ is 1.53% lower at the end of the experiment than at year 2170, as the response to decrease in emissions is lower than in upper waters.

**Regional patterns after the overshoot.** Differences between the behaviour of the global ocean during the last 30 years of the 1pctCO₂ simulation (After peak) and that of 100 years of the *piControl* simulation before the 1% ramp-up started (Before peak), allow us to explore regional differences in the reversibility of a given parameter. Global maps of the multi-model mean of this difference (After peak—Before peak) for temperature, pO₂ and Φ illustrate their deviation from pre-industrial levels (Fig. 2). In several regions, these parameters deviate from initial pre-
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Changes in the last 30 years of the 1pctCO₂ simulation relative to 100 years-average piControl simulation. Multi-model mean differences of temperature (a, °C), oxygen partial pressure (b, atm), and metabolic index (c) at the ocean surface (top panels), epipelagic (0–200 m; middle panels) and mesopelagic (200–1000 m; bottom panels) depth levels. Stippling in the maps marks disagreement across models on signs of change. Metabolic Index corresponds to the mean over 72 species (see Methods). The metabolic index scale is non-linear.</figcaption></figure>

industrial conditions at the ocean surface, in the epipelagic and, especially, mesopelagic depth layers, indicating the After peak global ocean differs from the global ocean at the beginning of the experiment.

Consistently across the seven ESMs (stippling indicates disagreement in the sign of change), the world ocean will generally be warmer after the overshoot at the sea surface except for small regions close to the Pacific sector of the Southern Ocean and the eastern coast of South Africa (Fig. 2a). Warming is also lower in the centre of the subpolar North Atlantic in the epipelagic and mesopelagic realms, where the so-called “warming hole” has been well documented in both observations and models; typically related to a persistent slowdown of the Atlantic meridional overturning circulation⁵⁵⁻⁵⁸. The general surface pattern is maintained through the epipelagic layer, while the warming with respect to the piControl situation is higher in the mesopelagic. The hysteresis in temperature on centennial time-scales is a consequence of an imbalance between the radiative forcing resulting from changes in atmospheric CO₂ levels and changes in ocean heat fluxes (Supplementary Fig. 3). Around model year 230, the global ocean heat fluxes become negative meaning that the ocean only releases part of the previously absorbed heat by the end of the simulation period⁵⁹,⁶⁰. Similarly, a lag in ocean heat uptake has been observed after emissions stop in recent studies (for example ref. ⁶¹).

Global pO₂ maps display more complex variability (Fig. 2b). At the sea surface, temperate and tropical regions will generally have lower pO₂ at the end of the experiment. In epipelagic waters, the pattern is similar to that displayed at the sea surface. In the mesopelagic, vast areas of the Atlantic and, especially, the Southern Oceans will have lower pO₂. The North Pacific shows the highest increases, in agreement with Li et al. ⁶², which found

increases in DO in this region as a result of changes in ocean circulation and of ventilation, which dominated over biological changes in DO concentrations.

Although Φ depends on temperature and pO₂, it displays some differences in regional variability due to the non-linear relationship that ties these two parameters together (Fig. 2c). At the sea surface, bands of lower Φ are seen in the polar regions of the global ocean resulting from the combined effects of positive (warmer) and negative (less oxygen) hysteresis of SST and pO₂, respectively. Some increases are seen in the Pacific sector of the Southern Ocean and part of the Arctic. This general pattern of habitability loss denoted by the decreasing values of Φ is largely reproduced in the epipelagic and mesopelagic. In the mesopelagic layer, a strong increase of Φ occurs in the subpolar North Pacific associated with the strong increase in pO₂, although model confidence in the mesopelagic Pacific is low, as indicated by the stippling. Globally, however, there is a net loss of habitability that increases with depth, as already seen in Fig. 1.

Supplementary Fig. 4d–f indicate changes between the global ocean after and before CO₂ peaks during the SSP5-3.4-OS. Though both experiments differ in the length and magnitude of the overshoot, regional changes in temperature, pO₂ and Φ bring some similar patterns at the sea surface and the epipelagic layer. The main differences occur in the North Atlantic where upper waters are cooler in the “warming-hole” region after the overshoot in the SSP5-3.4-OS. The pattern of changes in temperature in the mesopelagic broadly agree with those found for the 1pctCO₂, with larger differences in pO₂ and Φ. In these mesopelagic waters, most regions, including the North Pacific, present lower values of pO₂. Similarly, ocean habitability, as informed by Φ, reduces after the overshoot except in some regions of the tropical Pacific, Indian, and North Atlantic.
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Changes in the O₂,sat and AOU components of the metabolic index. Multi-model mean differences of the metabolic index at the ocean surface (top panels), epipelagic (0–200 m; middle panels) and mesopelagic (200–1000 m; bottom panels) depth levels. Metabolic index is decomposed into O₂,sat (a) and AOU components (b). Magenta/green indicates lower/higher metabolic indexes at the end of the experiment. A map of the absolute dominant effect of both components is included (c). Maps indicate differences in the last 30 years of the 1pctCO₂ simulation relative to 100 years-average of piControl simulation. Blue colours indicate the dominance of O₂,sat components. Orange colours indicate the dominance of AOU components. Lighter colours indicate the absolute contribution of both components is similar ([-0.25, 0.25]). Metabolic Index corresponds to the mean over 72 species (see Methods).</figcaption></figure>

**Drivers of changes in Φ.** To better understand the processes driving changes of Φ, we have decomposed it into changes driven by the saturation concentration of oxygen (O₂,sat) and those driven by changes in biological and ocean dynamical processes (see Methods). Here we assume that the oxygen concentration in the ocean is the result of (1) its solubility, which is dependent on water mass temperature and salinity⁶³, and (2) on the biological activity and ocean dynamics that can be approximated by the apparent oxygen utilisation (AOU), which accounts for oxygen consumption in a water parcel since its last contact with the atmosphere. Therefore, Φ can be computed using the O₂,sat partial pressure to extract the component that depends on it (Φsat). Then, by subtracting Φsat from Φ, we obtain the component that depends on AOU (ΦAOU).

Multi-model mean maps of After peak—Before peak differ-
ences during the idealised experiment for these two components
of Φ (Fig. 3a–b) highlight that negative hysteresis of Φ is mainly
driven by the O₂,sat component in most parts of the upper
(surface and epipelagic layers) world ocean, in close agreement
with positive hysteresis of temperature (residual warming). In
contrast, the AOU component is smaller but brings a spatial
imprint of the change in AOU. AOU components of Φ display
negative hysteresis in the subpolar regions of both hemispheres
across depth, and positive hysteresis in temperate and tropical
regions. Positive hysteresis of the AOU components of Φ is also
shown in the coastal region of the Antarctic. The magnitude of
the hysteresis of both components of Φ increases with depth.

Differences of these two components of Φ after and before the
overshoot simulated in the SSP5-3.4-OS simulation roughly

coincide with those simulated in the 1pctCO₂ experiment for the
two components in the three vertical layers considered as
informed in Supplementary Fig. 5d–e. Changes in O₂,sat dominate
positive changes of Φ on the North Atlantic upper waters in
agreement with the cooling of the region after the overshoot (see
Supplementary Fig. 4). In contrast, positive changes in AOU, as
observed for the 1pctCO₂ simulation in most temperate waters in
the mesopelagic, reduce as a function of latitude in the northern
and southern hemispheres.

We have computed the dominant component of Φ throughout
the water column (Fig. 3c and Supplementary Fig. 5c, f) by
computing the differences between their absolute values. If the
differences are smaller than ±0.25, we consider neither compo-
nent dominant (see Methods). In general, O₂,sat components drive
most of the changes in Φ during the 1pctCO₂ simulation, except
for regions where old waters accumulate and, consequently, the
AOU component dominates, such as in the mesopelagic layer of
the North Pacific. Overall, there is no clear dominance in tropical
regions at all depth levels, except for the Atlantic where O₂,sat
components dominate in mesopelagic waters. The patterns for
the dominance of O₂,sat vs AOU components for the SSP5-3.4-OS
(Supplementary Fig. 5f) are almost identical to those observed for
the 1pctCO₂ simulation for the upper layers. In contrast, the
regions with no clear dominance are higher with respect to
1pctCO₂ in the mesopelagic waters. The effects of AOU
components are not dominant in this scenario. Using lower
thresholds to define regions with no clear dominance also makes
the O₂,sat component dominant in tropical regions (see
Supplementary Figs. 6–8).
**Loss of marine ecosystem habitability.** For organisms to carry out aerobic metabolic activities, oxygen supply must meet at minimum the resting metabolic demand. While in theory, balance between demand and supply in oxygen occurs when unity is reached, that is, the resting rate ($\Phi_{\text{rest}}$; $\Phi = 1$), minimum values of $\Phi$ between 2 and 7 are observed in viable marine ecosystems³⁷. These values are thus considered as a critical metabolic index ($\Phi_{\text{crit}}$) that reflects the $\Phi_{\text{rest}}$ plus additional energy to undertake key ecological activities such as feeding, growth and reproduction³⁸,⁵¹. We track the fractional change of the global ocean volume of the epipelagic and mesopelagic waters at which $\Phi$ values are (1) lower than $\Phi_{\text{rest}}$, (2) between $\Phi_{\text{rest}}$ and $\Phi_{\text{crit}}$ and (3) higher than $\Phi_{\text{crit}}$ during the experiments (Fig. 4), and with respect to global warming (Supplementary Fig. 9). A compilation of 72 species (see Methods and ref. ³⁸) is used to compute quantiles of fractional change in ocean habitable volume. This compilation cannot be considered as an exhaustive representation of all marine species, but includes species of fishes, crustaceans, and molluscs, most of which are of economic relevance, living across a wide range of depths. Habitability is thus discussed in reference to this broadly subset of marine species.

We follow the evolution of the global ocean volume at which $\Phi$ is below the resting rate ($\Phi < \Phi_{\text{rest}}$), that is, the volume of the world ocean that is unsuitable for supporting resting aerobic metabolism. We call this zone the ‘non-viable zone’. Further, we follow the evolution of the ‘critical zone’ characterised by $\Phi$ values within the resting and the critical rates ($\Phi_{\text{rest}} \le \Phi < \Phi_{\text{crit}}$). These waters can be considered as transitional since some marine organisms may develop acclimatisation or genetic adaptation strategies to cope with low $O_2$ levels (see for example ref. ⁶⁴). Finally, we follow the evolution of the ‘viable zone’, which corresponds to the volume of the global ocean above the critical rate ($\Phi \ge \Phi_{\text{crit}}$), where most organisms can thrive.

In the upper 200 m of the water column (Fig. 4, top row), the global ocean's volume of non-viable zone waters represented ~8% of epipelagic volume at the beginning of the 1pctCO₂ simulation, and it increases to be 28.8 ± 1.9% larger after the overshoot. As already observed for marine biodiversity heat stress⁶⁵, the increase in the volume of non-viable waters scales rapidly with global warming (Supplementary Fig. 9), while the global cooling flattens during the ramp-down slowing down the recovery. This is possibly caused by the long responding timescale of the ocean to reductions in atmospheric $CO_2$, especially when the overshoot is large and persists for a long period⁶². The critical zone also increases from the initial volume (~36%), showing a hysteresis of 9.1 ± 0.2% at the end of the experiment. The rate of loss of these waters is also slightly larger during the ramp-up phase of the experiment than during the ramp-down phase. The volume of habitable waters, that is, those at which $\Phi$ is higher than $\Phi_{\text{crit}}$, decreases steadily during the ramp-up, reaching a loss of >16% when it reaches a minimum around model year 150. These waters represented ~56% of the total volume of epipelagic waters and declined by about ~5.2% at the end of the experiment with respect to the volume at the start of the experiment (Fig. 4c).

In the mesopelagic waters (Fig. 4, second rows), the mean loss in habitability is smaller than in epipelagic waters. At the beginning of the experiment, the non-viable and the critical zones represented ~13% and ~29% of the mesopelagic volume respectively. At this depth level, there is almost no overall expansion in volume of the non-viable zone (1.2 ± 0.3%), while the initial critical zone increases by 6.0 ± 0.2%. However, the distribution of both ocean volumes increase to up to 150%, almost doubling the initial volume. The viable zone represented ~58% of the volume of mesopelagic waters, with a 4.6 ± 0.1% loss at the end of the experiment. As with non-viable and critical zones, the viable zone distribution indicates much higher

losses of ocean volume -of up to 70%- with respect to the initial volume.

We carry out the same exercise for the SSP5-3.4-OS. However, as $CO_2$ concentrations do not return to initial conditions, we modify the approach slightly, by selecting the period just before the overshoot (years 2035–2055) as reference. In the upper 200 m of the water column (Fig. 4, third rows), the global ocean's volume of non-viable zone waters represented ~9% of epipelagic volume at the baseline period, increasing to be 12.7 ± 2.1% greater at the end of the scenario. This is also the case for the critical zone that increases 1.1 ± 0.2% from the baseline volume (~38%). Habitable waters reduce their baseline volume (~53%) by about 1.2 ± 0.1%. Interestingly, this loss occurs after the volume of the viable zone returns to about pre-overshoot conditions as decreasing $CO_2$ levels.

In the mesopelagic (200–1000 m, Fig. 4 fourth panels), habitability volume loss is slightly higher than in upper waters as contrasted with observed for the 1pctCO₂ simulation since the timescales of both overshoots are different. The volume of the non-viable and critical zones represented ~14% and ~31% of the mesopelagic volume, and increase by 12.2 ± 0.3% and 5.9 ± 0.1% respectively at the end of the simulation. This loss in habitability is also observed by a decrease in the viable zone that represented ~55% and declined by -4.9 ± 0.1% at the end of the experiment. The mesopelagic time series until year 2170 demonstrates similar behaviour as that observed during 1pctCO₂, suggesting the overshoot signal takes long to propagate into deep waters. In contrast, the volume of viable waters in the epipelagic realm follows more closely the $CO_2$ concentration.

Changes in the volume of waters that can or cannot hold marine habitats to the 72 species listed here differ regionally. For the 1pctCO₂ experiment, we have computed the global ocean multi-model After peak—Before peak (as defined for Figs. 2, 3) differences in volume for the epipelagic and mesopelagic depth layers for each of the three types of waters (Fig. 5 upper two rows). In epipelagic waters, the volume of non-viable waters increase both at the tropical and temperate regions. The gain in critical waters occurs mainly in temperate regions, compensated by losses in the tropics. The map showing the differences in the viable zone displays an overall loss in habitability across the global ocean except for regions such as the so-called warming hole (see above) and the North Pacific. In the mesopelagic, the pattern of changes in viable volumes are amplified. For example, the tropical Atlantic will lose habitable volume by the increase of non-viable waters and the decrease of the viable zone, but the North Pacific will gain habitable volume as $pO_2$ is higher (see Fig. 2) by the end of the experiment. The volume of non-viable waters strongly increase in the subtropical Pacific. However, a concomitant increase in the critical zone volume may indicate that those organisms capable of coping with low levels of oxygen may still thrive in this region after the simulation period.

For the SSP5-3.4-OS, we have computed regional differences in the habitable volume between after and before overshoot, i.e., years 2080–2100 minus years 2035–2055 (Fig. 5 lower two rows). Regional patterns coincide with those observed for the 1pctCO₂ simulation in epipelagic waters, though relative changes are weaker. In mesopelagic waters, main differences occur in the North Pacific as the increase in oxygen supply after recovering from the overshoot is not simulated under the SSP5-3.4-OS scenario (see Supplementary Fig. 4), leading to net loss in viable waters in this region.

**Reversibility and long-term commitment.** We can now discuss the level of (ir)-reversibility of temperature, $pO_2$, and $\Phi$ by tracking the deviation from initial values and the trends over last
<figure><img src="image_5.png" /><figcaption>**Fig. 4 Change in ocean habitability volume.** Model-mean time series of the fractional change of global ocean volume where $\Phi < \Phi_{\text{rest}}$ (non-viable zone; **a**, **d**), $\Phi_{\text{rest}} \ge \Phi < \Phi_{\text{crit}}$ (critical zone; **b**, **e**), and $\Phi \ge \Phi_{\text{crit}}$ (viable zone; **c**, **f**). Time series for 1pctCO₂/SSP5-3.4-OS experiments are shown in the upper/lower panels. For 1pctCO₂, the fractional change is computed relative to a 100 year-average of piControl simulation. For SSP5-3.4-OS, the fractional change is computed relative to the years 2035 to 2055, just before the peak in CO₂ concentrations. Time series are shown for epipelagic (0–200 m; top panels) and mesopelagic (200–1000 m; bottom panels) depth layers. Red solid (dashed) lines illustrate the mean (median) over 72 species and seven (1pctCO₂)/four (SSP5-3.4-OS) ESMs. Red shading indicates the multimodel’s 72 species’ full distribution: 1, 2.5, 5, 10, 17, 25, 75, 83, 90, 95, 97.5 and 99th percentiles. ‘Vol’ indicates the volume that each type of water (non-viable, critical and viable zones) represents to the total volume of the epipelagic and mesopelagic layers at the reference periods. ‘$\Delta vol$’ indicates the level of hysteresis depicted as the difference between ocean habitability volume at the last 10 years of the simulations relative to the reference period for each simulation.</figcaption></figure>

decades of the 1pctCO₂ experiment for each parameter (see Methods). We have combined in a diagram deviation-slope pairs for each ESM and depth level (Fig. 6). According to recent studies (for example ref. ⁶²,⁶⁶), several processes linked with ocean warming and deoxygenation are largely reversible (within

centuries) at the sea surface. At the same time, they will take centuries to millennia to recover at depth²¹,⁶⁷⁻⁷¹. To determine whether this is the case in our study, we also estimate the same deviation-slope pairs arising from internal climate variability. A change is considered as irreversible at centennial timescales for a
<figure><img src="image_6.png" /></figure>

**Fig. 5** Regional changes on ocean habitability volume. Multi-model mean differences of the global ocean volume where $\Phi < \Phi_{\text{rest}}$ (non-viable zone; **a**, **d**), $\Phi_{\text{rest}} \ge \Phi < \Phi_{\text{crit}}$ (critical zone; **b**, **e**), and $\Phi \ge \Phi_{\text{crit}}$ (viable zone; **c**, **f**) at epipelagic (0 - 200 m; top panels) and mesopelagic (200 - 1000 m; bottom panels). For 1pctCO₂ (upper two rows), maps indicate differences in the last 30 years of the experiment relative to 100 years-average of piControl simulation. For SP5-3.4-OS (lower two rows), maps indicate differences after CO₂ peak (years 2080 to 2100) relative to before CO₂ peak (years 2035 to 2055). Metabolic Index corresponds to the mean over 72 species (see Methods).

given ESM, depth level (and species in the case of $\Phi$) if the slope of the last 30 years of the model simulations falls within the range of the internal variability whereas its deviation doesn't (light red shading boxes). A change is considered as reversible if the deviation falls within internal variability (light green shading boxes).

In epipelagic waters, ocean temperature does not return to initial values after the temperature overshoot within the experiment timescale. Ocean temperatures are still adjusting in mesopelagic waters, indicated by both trend and deviation parameters being far from the reference values defined by piControl. This indicates that the time to recover in deeper water is proportional to the strength of the overshoot, consistent with literature³⁰,⁶². Similarly, pO₂ trend-deviation pairs indicate irreversible behaviour for pO₂ at both surface and epipelagic waters. Multi-model uncertainty on pO₂⁷² results in positive and negative deviations at the mesopelagic realm. Nonetheless, ESMs in which global mesopelagic pO₂ increases or decreases after the experiment show irreversibility at centennial timescales except for MIROC-ES2L, for which the trend and deviation fall within natural variability.

Figure 6c shows each ESM's slope/deviation pairs for the 72 species median $\Phi$. In general, all ESMs show irreversible behaviour at the sea surface and the upper 200 m depth. In addition, in some models, particularly CNRM-ESM2-1 and

UKESM-1-0-LL, $\Phi$ is still adjusting. In mesopelagic waters, the number of ESMs at which $\Phi$ is still adjusting increases, as seen in ACCESS-ESM1-5 and NorESM2-LM.

Unlike the 1pctCO₂ simulation in which CO₂ levels return to initial pre-industrial conditions, the SSP5-3.4-OS scenario simulation proposes a different state after recovering from the overshoot. We make use of the extension of the SSP5-3.4-OS after the overshoot to give insights into the regional reversibility behaviour of the ocean's habitability. Figure 7 follows the regional changes in ocean habitable volume in the last 20 years of the experiment with respect to after the CO₂ overshoot (years 2080-2100). In epipelagic waters, the volume of habitable waters generally increases after the overshoot, as opposed to the decrease in habitability observed when comparing after and before the overshoot (Fig. 5). However, even after ~200 years post overshoot, the general pattern towards recovering ocean habitability is relatively low. Changes are stronger at depth. In the North Pacific mesopelagic waters, the pattern shifts from losing to gaining habitability by the increase of viable waters and the decrease on non-viable waters. A similar finding was observed for the After peak—Before peak differences of the idealised experiment (Fig. 5), in which an increase in pO₂ occurred in the region driving the increase in habitability. In the tropical regions, the volume of the critical zone decreases. In contrast, the volume of non-viable waters increases in the Southern ocean,
<figure><img src="image_7.png" /><figcaption>**Fig. 6 Trend-deviation diagram.** Last thirty year trend (years 311-340, x axis) and deviation from the initial value (y axis) of the 1pctCO₂ experiment for T (a, d, g), pO₂ (b, e, h) and metabolic index (c, f, i) at the sea surface (top panels), epipelagic (0-200 m; middle panels) and mesopelagic (200-1000 m; bottom panels) depth levels for each ESM (bold dots). Trend-deviation pairs are also computed for 25 randomly sampled 30 year periods from the *piControl* simulations of each ESM. Light green rectangles indicate reversibility on centennial timescales as trends fall within *piControl* internal variability. Light red rectangles indicate irreversibility at longer timescales as the slope is within *piControl* internal variability, while deviation falls outside *piControl* internal variability. The median over 72 species is displayed (bold dots) for each ESM for metabolic index. Slope and deviation for pO₂ are in 10⁻³ atm.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7 Regional changes in ocean habitable volume after the overshoot simulated in SSP5-3.4-OS.** Multimodel mean differences of the global ocean volume where $\Phi < \Phi_{\text{rest}}$ (non-viable zone; a), $\Phi_{\text{rest}} \ge \Phi < \Phi_{\text{crit}}$ (critical zone; b), and $\Phi \ge \Phi_{\text{crit}}$ (viable zone; c) at epipelagic (0-200 m; upper panels) and mesopelagic (200-1000 m; lower panels) waters. Maps indicate differences in the last 20 years of the experiment relative to years 2080 to 2100, immediately post the peak of CO₂ concentrations. Metabolic Index corresponds to the mean over 72 species (see Methods).</figcaption></figure>
limiting habitability in that important region. Altogether, these maps suggest that the volume of habitable waters increases when CO₂ levels recover from the overshoot. However, though the response to the atmospheric CO₂ forcing is more rapid in upper waters, this recovery may take centuries to accomplish. Furthermore, many deeper ocean regions will continue to lose habitability, reflecting the legacy of the global warming overshoot propagating into deep waters.

**Summary and implications.** The experiments presented here allow us to study the ocean system's reversibility after hypothetical global warming overshoot levels in the Earth system. Through the evolution of a metabolic index ($\Phi$), this study provides insights into the marine organisms' ecophysiological response to changes in ocean temperature and pO₂. Using a suite of ESMs, we have shown that the level of hysteresis in the upper layers of the ocean may be higher for $\Phi$ than for its drivers (Figs. 1, 6), potentially affecting the resilience of marine habitats if CO₂ levels were to decline. As pO₂ is directly related to temperature for a given O₂ concentration, the limited solubility of O₂ as the ocean remains warmer than during the pre-industrial period after the overshoot (Fig. 2) mainly drives this negative hysteresis (Fig. 3).

Moving into the ocean interior, the lagged response of temperature increase results in a warmer global ocean after the overshoot (Fig. 2). Though the global hysteresis of pO₂ reverses in the ocean interior during the 1pctCO₂ simulation (Supplementary Fig. 2), there are regional differences: in the Indian and Pacific Oceans, increases in pO₂ are simulated in the mesopelagic realm, while in the Atlantic and Southern Oceans, pO₂ decreases after the overshoot. In these deep waters (200–1000 m), values of $\Phi$ decrease in the global ocean, except in the North Pacific. Changes in the sea ice cover and vertical stratification due to warming (Fig. 2) may amplify the vertical mixing, thereby ventilating in oxygen the deeper layer of the ocean. The hysteresis of pO₂ is less obvious for the SSP5-3.4-OS scenario (Supplementary Fig. 4) as the overshoot duration and magnitude is shorter; hence, the warming intensity is not enough to produce these such large perturbations.

The decrease of $\Phi$ after the temperature overshoot suggests the capacity of the ocean to hold viable marine ecosystems also declines. Our analysis indicates that the global ocean will lose on average 4–5% of habitable waters after the overshoot above 1000 m depth during the 1pctCO₂ simulation, and 1 up to 5% on average during the SSP5-3.4-OS. Furthermore, the volume of waters whose level of oxygen is insufficient to maintain life increases significantly (~13–29% for the SSP5-3.4-OS and the 1pctCO₂ respectively) in the upper 200 m where most of the biota live<sup>73</sup>, limiting altogether the resilience of shallow marine ecosystems. In these upper waters, habitable waters will be lost across the global ocean, and marine organisms' resilience will be impacted. In mesopelagic waters, even though the magnitude of habitable volume loss is higher in the Atlantic, it is offset by strong gains in the North Pacific (Fig. 5). The magnitude of these gains is lower in the SSP5-3.4-OS, as the ocean state does not come back to initial conditions but moves into a new quasi-stable state for at least two hundreds of years. (Fig. 7).

As dependent on temperature and pO₂ extracted from ESM simulations, it is worth to specify that values of $\Phi$ are sensitive to model biases. Nonetheless, and irrespective of the ESM considered, the multi-model response across the metabolic index computed for 72 different species provides a robust and consistent feature that brings evidence on the loss in ocean habitability after an overshoot.

These global warming induced alterations of the habitability of marine ecosystems are likely to require centuries to recover after

the peak in CO₂ concentrations even if the anthropogenic climate forcing is removed, that is, by returning CO₂ concentrations back to pre-industrial levels (Fig. 6). These timescales may also be longer than the time needed for temperature and pO₂ to recover back to pre-industrial levels in the upper layers of the ocean. As our analysis considers an idealised overshoot scenario consisting of a rapid increase and decrease in atmospheric CO₂, the overshoot and implementation (if any) of negative emissions following CDR deployment would be expected to take much longer. Under a more comprehensive overshoot scenario, the loss of habitability starts to halt after the overshoot, yet the recovery takes a long time even though the magnitude and extension of the atmospheric CO₂ overshoot are lower. Together, this means that the hysteresis and recovery times might be considerable longer than that estimated in this study and are likely underestimated.

Using an ecophysiological index that combines aerobic metabolism and cold temperature tolerances of global marine biodiversity, Penn and Deutsch<sup>52</sup> quantified the risk of extinction for marine organisms as global warming continues during the present century. As equivalent to the habitability volume evaluated here, they assessed local extirpations as wherever $\Phi$ falls below $\Phi_{crit}$ and found larger impacts in the tropics. Our results agree suggesting that the volume of non-viable waters after passing an overshoot will increase in the tropical oceans where species are living closer to their limits (see also ref. <sup>74</sup>). Furthermore, a mechanistic model developed by Deutsch et al.<sup>51</sup> suggests that reduction in $\Phi$ may be compensated by a reduction in body size for those organisms for which their metabolic demand increases with temperature. This implies these species may still be capable of occupying those regions at which $\Phi$ falls below critical values in detriment of their body size. In addition, application of novel indices that are not based on experimentally derived temperature-dependent O₂ thresholds such as the Aerobic Growth Index<sup>75,76</sup>, or the use of ecophysiological traits defined to represent marine biodiversity (see ref. <sup>39,52</sup>) has the potential to expand our framework to a wider range of marine species. Nonetheless, the set of available ecophysiological parameters used here provide a novel mapping for economically-relevant species that can be found in the global ocean and across a wide range of depths.

Notwithstanding with the above-mentioned limitations, we show that the combined effect of ocean warming and deoxygenation after a temperature overshoot will impact ocean habitability by increasing the volume of waters that are unsuitable for supporting marine ecosystems at the expense of those waters that can hold them. We found this emerging picture is robust across a range of idealised and comprehensive overshoot simulations. Even if temperature returned to pre-industrial levels, as simulated in the idealised experiment, these habitability loss effects will likely be felt for centuries. The lagged effects of global warming on marine ecosystems after CO₂ levels recover from overshoot found here underlines that mitigation actions should be promptly implemented as any delay in stabilising climate are highly likely to have an irreversible impact on marine ecosystems and the vital ecosystem services they provide.

## Methods

**Simulations.** The CDRMIP cdr-reversibility experiment (1pctCO₂-cdr) is fully described in Keller et al.<sup>16</sup>, and only the main characteristics are briefly described here. It corresponds to a highly idealised experiment that investigates the reversibility of the climate system by leveraging the prescribed 1% yearly increase in atmospheric CO₂ concentrations (ramp-up) from pre-industrial levels (~285 ppm), until quadrupling of CO₂ concentrations as designed for the 1pctCO₂ standard DECK run for CMIP6<sup>42,77</sup>. At this concentration, it prescribes a symmetric 1% yearly decline in
atmospheric CO₂ concentrations (ramp-down) until reaching initial conditions (pre-industrial level), at which point concentrations are held constant for an additional 60 years. Peak atmospheric CO₂ concentration of ~1140 ppm is reached after 140 years of experiment. Therefore, 280 years of the combined ramp-up—ramp-down experiments are considered in addition to 60 years of constant concentrations.

The socio-economic shared scenario SSP5-3.4-OS explores a peak and decline in CO₂ concentrations during the 21st century. Carbon emissions during this scenario initially starts by following a high emission pathway consistent with SSP5-8.5 until 2040, at which point strong mitigation policy including CDRs is undertaken, reducing emissions to zero by around 2080. Then, negative emissions continue until 2140. Net zero is again reached by 2170²³,²⁴. SSP5-3.4-OS is extended to year 2300. As CO₂ concentrations do not return to initial values, reversibility is assessed by comparing a period after the peak in CO₂ (2080–2100, when CO₂ levels are ~500 ppm) with respect to 20 years before the peak on CO₂ that occurs around year 2062 considering the mean across the four ESMs considered in which CO₂ concentrations are similar to those after the overshoot. These years correspond to 2035–2055 when CO₂ levels vary between 475 and 564 ppm.

Outputs from a suite of six CMIP6 generation and one CMIP5 generation ESM models participating in the CDRMIP project are used in this study. They are listed in Supplementary Table 1, and described in Supplementary Note 1. All simulations provide monthly 3D outputs. Here, only those vertical layers describing the first 1000 m depth of the water column are used. Multi-model mean and uncertainty are displayed throughout the study. Following previous assessments²⁸,⁷⁸, inter-model uncertainty is estimated as the inter-model standard deviation (SD).

A pre-industrial control (without anthropogenic forcing) one-member 500 yearlong simulation (*piControl*) is used for each ESM to extrapolate the level of global warming as the difference between atmospheric temperature during the experiment and the long-term *piControl* mean. Likewise, 100 years of these *piControl* experiments are used as the Before Peak period (see Figs. 2, 3, and 5).

To facilitate intercomparison, model outputs were re-gridded from their model grids to a regular 1° × 1° horizontal grid using distance weighted average remapping (climate data operators; *remapdis*). Likewise, vertical coordinates were all interpolated to the 5, 10, 20, 30, 50, 75, 100, 125, 150, 200, 250, 300, 400, 500, 600, 700, 800, 900 and 1000 m depth levels (climate data operators; *intlevel*).

**Estimation of metabolic index.** The supply of oxygen (S) to a marine species should be higher than the resting metabolic oxygen demand (D) in order to be metabolically viable. As a result of global warming, the solubility of oxygen in upper waters is reduced as warmer waters saturate quicker, while the vertical stratification increases reducing mixing of surface oxygen-rich waters with subsurface oxygen-poor waters⁷⁹. These processes reduce the concentration of DO in the upper ocean²⁶. In addition, warming increases the metabolism of organisms⁸⁰,⁸¹ thus making them more oxygen demanding. Therefore, there is a need for a metric that makes a link between global warming perturbations and their effects on the ocean distribution of life. In this regard, the aerobic energy balance can be represented by a Metabolic Index³⁷,³⁸ (*Φ*), which represents the ratio of oxygen supply to resting demand.

$$ \Phi = \frac{S}{D} $$

The rate of oxygen supply depends on ambient partial pressure of oxygen (pO₂) and on respiratory efficacy⁸², which is equal to a

per-mass rate of gas transfer between water and the organism (α₈), and its scaling with body mass (B⁶). The denominator of metabolic index depends as well on body mass (B⁶) and on the absolute temperature (T). Also, it depends on species specific terms of metabolic rate (αₔ), its temperature dependence (E₉), and to the Boltzmann's constant⁸³ (kₐ), which captures the temperature dependence often described by a Q10 factor⁸³.

$$ S = \alpha_{s}(T)B^{\sigma}pO_{2} $$

$$ D = \alpha_{D}B^{\delta} \exp \left( \frac{-E_{d}}{k_{B}} \left( \frac{1}{T} - \frac{1}{T_{ref}} \right) \right) $$

Because diffusive gas (as O₂) fluxes are governed by physicochemical kinetics, the temperature dependence of the oxygen supply (α₈(T)) can be approximated by the Arrhenius function that reflects standard gas exchange across a diffusive boundary layer model⁸⁴.

$$ \alpha_{s}(T) = \alpha_{s} \exp \left( \frac{-E_{s}}{k_{B}} \left( \frac{1}{T} - \frac{1}{T_{ref}} \right) \right) $$

Therefore, the metabolic index includes both environmental and physiological variables as

$$ \Phi = \frac{S}{D} = A_{o}B^{\epsilon}pO_{2} \exp \left( \frac{E_{o}}{k_{B}} \left( \frac{1}{T} - \frac{1}{T_{ref}} \right) \right) $$

Where A₀ is the ratio of rate coefficients for oxygen supply and metabolic demand (α₈/αₔ), and ε is the difference between the allometric scales of body size and mass (σ – δ). E₀ (E₉—E₈) describes the effect of temperature on the critical pO₂, which corresponds to the minimum pO₂ required for maintaining a metabolic rate equals to 1, that is, the resting metabolic demand. To take into account the uncertainties in the computation of the evolution of Φ caused by species-specific ecological parameters (E₀, A₀), Φ is computed for a total of 72 marine species compiled by Deutsch et al. ³⁸. For most marine organisms, the allometric scales of body size (ε) tend to 0, thus B⁶ is considered to be 1.

For environmental parameters, we use ocean temperature and DO from the simulations described above to estimate Φ. pO₂ is calculated as the ratio between the oxygen concentration and the Henry's law coefficient for oxygen saturation following Garcia and Gordon⁶³. Maps of Φ computed using observations from the World Ocean Atlas 2018 and the seven ESMs *piControl* simulations are presented in Supplementary Fig. 10. World Ocean Atlas 2018 observations can be accessed via https://www.ncei.noaa.gov/access/world-ocean-atlas-2018/.

**Decomposition of the metabolic index.** In the ocean, O₂ concentrations are the result of the solubility of O₂, which depends on the physical properties of water masses (temperature and salinity), and the effects of ocean dynamics (ventilation and stratification) and biological activity that can be approximated by changes in AOU. The solubility of O₂ can be obtained by changes in O₂ saturation (O₂,sat), which indicates the maximum O₂ concentration a parcel of water can contain. O₂,sat is computed from ESMs temperature and salinity outputs, and represents the effect of changes in oxygen solubility on dissolved O₂ concentration. Therefore, O₂ can be decomposed as O₂ = O₂,sat—AOU.

Following the decomposition of O₂, the metabolic index can thus be decomposed in a term mainly affected by O₂,sat, and a
term driven by AOU.

$$ \Phi = f(pO_2, T) = A_o B^e pO_2 \exp\left(\frac{E_o}{k_B}\left(\frac{1}{T} - \frac{1}{T_{ref}}\right)\right) $$

Given that

$$ pO_2 = \frac{O_2}{sP_{dep}} = \frac{O_{2,sat} - AOU}{sP_{dep}} $$

Being $sP_{dep}$ the pressure correction which is temperature dependent, and that $B^e \sim 1$, then:

$$ \Phi = \frac{A_o O_{2,sat} - A_o AOU}{sP_{dep}} \exp\left(\frac{E_o}{k_B}\left(\frac{1}{T} - \frac{1}{T_{ref}}\right)\right) \\ \Phi = \frac{A_o O_{2,sat}}{sP_{dep}} \exp\left(\frac{E_o}{k_B}\left(\frac{1}{T} - \frac{1}{T_{ref}}\right)\right) - \frac{A_o AOU}{sP_{dep}} \exp\left(\frac{E_o}{k_B}\left(\frac{1}{T} - \frac{1}{T_{ref}}\right)\right) \\ \Phi = \Phi_{sat} - \Phi_{AOU} $$

We thus compute $\Phi_{sat}$ by using $O_{2,sat}$ partial pressure. Finally, we subtract $\Phi_{sat}$ from $\Phi$ to obtain $\Phi_{AOU}$.

We have also computed the relative dominance of each of $\Phi$ components throughout depth levels by computing their difference in absolute values. If the differences are above a given threshold in absolute values we consider a component to be dominant. If these are below the thresholds in absolute values, we consider none of the components dominant. We have used the median of the multi model SD of $\Phi$ at the sea surface, that roughly corresponds to 0.25, as the threshold. In order to assess the effect of different thresholds, we also computed the dominance approach using 0.5, 0.1 and 0.05 as thresholds (Supplementary Figs. 6–8). Increasing the threshold to consider no clear dominance between components of $\Phi$ mainly affects the dominant role of the $O_{2,sat}$ term in tropical regions. In these regions, no clear dominance between $O_{2,sat}$ and AOU terms dominate widely as thresholds increase.

**Determination of missing $\Phi_{crit}$**. In this study, a compilation of ecological parameters ($\Phi_{crit}$, $V_h$ and $E_o$) provided by Deutsch et al. 38, for 72 different species is used. Nonetheless, $\Phi_{crit}$ is only available for 58 species, thus it needs to be determined for 14 species. To approximate the values of $\Phi_{crit}$ for these latter species, we have computed the linear regression between $\Phi_{crit}$ and the species-specific hypoxia vulnerability (1/V$_h$; inverse of $A_o$) of the rest of the species. Then, the slope and intercept of those linear regression were used to derive missing $\Phi_{crit}$.

**Estimation of ocean volume for each type of waters.** For each model (seven for 1pctCO$_2$ and four for SSP5-3.4-OS simulations) and species (72), we compute, at epipelagic and mesopelagic depth levels, the ocean volume at which $\Phi$ is below $\Phi_{rest}$, that is, non-viable zone; the ocean volume at which $\Phi$ is between $\Phi_{rest}$ and $\Phi_{crit}$, that is, critical zone; and the ocean volume at which $\Phi$ is above $\Phi_{crit}$, that is, viable zone. Once the volume is estimated for each ocean zone, we compute -during the 1pctCO$_2$ simulation- its fraction with respect to total volume of the depth level at the beginning of the experiment, and the fractional change with respect to piControl simulation. For the SSP5-3.4-OS, we compute its fraction with respect to total volume of the depth level just before the overshoot (years 2035–2055), and the fractional change with respect to this reference period. Being $V\Phi_x$ the ocean volume at each depth level for the three (x) types of waters at a given time and $V\Phi_{x_0}$ the reference period mean value of $V\Phi_x$, the fractional change can be estimated as: $(V\Phi_x - V\Phi_{x_0})/V\Phi_{x_0}$. We multiply by 100 to express it in percentage.

**Reversibility timescale analysis.** To estimate the time for recovering initial values for temperature, $pO_2$ and the metabolic index, we track the slope of the linear regression of the last 30 years of the 1pctCO$_2$ simulation, and plot against the difference between the mean value over these last 30 years of the simulation and the initial value of the simulation. To discern whether these parameters display irreversibility, we also estimate the same quantities arising from internal climate variability by using 25 randomly selected 30 years-long subsamples of the piControl simulation for each model. A change is considered as irreversible at these centennial timescales for a given depth level and model (and species in the case of $\Phi$) when the slope is within internal variability whereas its deviation doesn't. In contrast, a change is considered to be reversible at these centennial timescales when both slope and deviation fall within the internal variability as diagnosed from the piControl.

## Data availability

Data presented in the paper can be accessed via Zenodo at https://doi.org/10.5281/zenodo.8279001.

## Code availability

All code used in the current study is available from the corresponding author upon reasonable request.

Received: 20 December 2022; Accepted: 14 September 2023; Published online: 10 October 2023

## References

1. Friedlingstein, P. et al. Persistent growth of CO$_2$ emissions and implications for reaching climate targets. *Nat. Geosci.* **7**, 709–715 (2014).

2. Peters, G. P. et al. Carbon dioxide emissions continue to grow amidst slowly emerging climate policies. *Nat. Clim. Change* **10**, 3–6 (2020).

3. Matthews, H. D. & Wynes, S. Current global efforts are insufficient to limit warming to 1.5 °C. *Science* **376**, 1404–1409 (2022).

4. Rogelj, J. et al. (2018). Chapter 2: Mitigation Pathways Compatible with 1.5 °C in the Context of Sustainable Development. in *Global Warming of 1.5 °C: An IPCC special report on the impacts of global warming of 1.5 °C above pre-industrial levels and related global greenhouse gas emission*.

5. Lee, J.-Y., et al, Future global climate: scenario-based projections and near-term information. In *Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the 6th Assessment Report of the Intergovernmental Panel on Climate Change* [eds Masson-Delmotte, V., P. Zhai, A. Pirani, S.L. Connors, C. Péan, S. Berger, N. Caud, Y. Chen, L. Goldfarb, M.I. Gomis, M. Huang, K. Leitzell, E. Lonnoy, J.B.R. Matthews, T.K. Maycock, T. Waterfield, O. Yelekçi, R. Yu, and B. Zhou]. 553–672 (Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 2021) https://doi.org/10.1017/9781009157896.006.

6. Gasser, T., Guivarch, C., Tachiiri, K., Jones, C. D. & Ciais, P. Negative emissions physically needed to keep global warming below 2. *Nat. Commun.* **6**, 7958 (2015).

7. Rockström, J. et al. The world's biggest gamble. *Earth's Future* **4**, 465–470 (2016).

8. Krieger, E. et al. Pathways limiting warming to 1.5 °C: a tale of turning around in no time? *Philos. Trans. R. Soc. A Math. Phys. Eng. Sci.* **376**, 20160457 (2018).

9. Lawrence, M. G. et al. Evaluating climate geoengineering proposals in the context of the Paris Agreement temperature goals. *Nat. Commun.* **9**, 3734 (2018).

10. Shepherd, J. G. *Geoengineering the climate: science, governance and uncertainty*. https://eprints.soton.ac.uk/156647/ (2009).

11. Scott, V., Haszeldine, R. S., Tett, S. F. B. & Oschlies, A. Fossil fuels in a trillion tonne world. *Nat. Clim. Change* **5**, 419–423 (2015).

12. Smith, P. Soil carbon sequestration and biochar as negative emission technologies. *Glob. Change Biol.* **22**, 1315–1324 (2016).

13. van Vuuren, D. P. et al. Carbon budgets and energy transition pathways. *Environ. Res. Lett.* **11**, 075002 (2016).

14. Minx, J. C. et al. Negative emissions—Part 1: Research landscape and synthesis. *Environ. Res. Lett.* **13**, 063001 (2018).
15. Jones, C. D. Chapter 3. Negative emissions: the role and response of the climate system. 27–56 (2022). https://doi.org/10.1039/9781839165245-00027.

16. Keller, D. P. et al. The Carbon Dioxide Removal Model Intercomparison Project (CDRMIP): rationale and experimental protocol for CMIP6. Geosci. Model Dev. 11, 1133–1160 (2018).

17. Frölicher, T. L. & Joos, F. Reversible and irreversible impacts of greenhouse gas emissions in multi-century projections with the NCAR global coupled carbon cycle-climate model. Clim. Dyn. 35, 1439–1459 (2010).

18. Schwinger, J. & Tjiputra, J. Ocean carbon cycle feedbacks under negative emissions. Geophys. Res. Lett. 45, 5062–5070 (2018).

19. Jeltsch-Thömmes, A., Stocker, T. F. & Joos, F. Hysteresis of the Earth system under positive and negative CO 2 emissions. Environ. Res. Lett. 15, 124026 (2020).

20. Kug, J.-S. et al. Hysteresis of the intertropical convergence zone to CO₂ forcing. Nat. Clim. Change. 12, 47–53 (2022).

21. Bertini, L. & Tjiputra, J. Biogeochemical timescales of climate change onset and recovery in the North Atlantic interior under rapid atmospheric CO₂ forcing. J. Geophys. Res. Ocean. 127, e2021JC017929 (2022).

22. Boucher, O. et al. Reversibility in an Earth System model in response to CO₂ concentration changes. Environ. Res. Lett. 7, 024013 (2012).

23. O'Neill, B. C. et al. The Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6. Geosci. Model Dev. 9, 3461–3482 (2016).

24. Meinshausen, M. et al. The shared socio-economic pathway (SSP) greenhouse gas concentrations and their extensions to 2500. Geosci. Model Dev. 13, 3571–3605 (2020).

25. Tebaldi, C. et al. Climate model projections from the Scenario Model Intercomparison Project (ScenarioMIP) of CMIP6. Earth Syst. Dyn. 12, 253–293 (2021).

26. Keeling, R. F., Körtzinger, A. & Gruber, N. Ocean deoxygenation in a warming world. Ann. Rev. Mar. Sci. 2, 199–229 (2010).

27. Gruber, N. Warming up, turning sour, losing breath: ocean biogeochemistry under global change. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 369, 1980–1996 (2011).

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

29. Bijma, J., Pörtner, H.-O., Yesson, C. & Rogers, A. D. Climate change and the oceans—What does the future hold? Mar. Pollut. Bull. 74, 495–505 (2013).

30. John, J. G., Stock, C. A. & Dunne, J. P. A more productive, but different, ocean after mitigation. Geophys. Res. Lett. 42, 9836–9845 (2015).

31. Levitus, S. et al. Global ocean heat content 1955–2008 in light of recently revealed instrumentation problems. Geophys. Res. Lett. 36, L07608 (2009).

32. Cheung, W. W. L. et al. Shrinking of fishes exacerbates impacts of global ocean changes on marine ecosystems. Nat. Clim. Change 3, 254–258 (2013).

33. Eaton, K. M., Hallaj, A., Stoeckel, J. A. & Bernal, M. A. Ocean warming leads to increases in aerobic demand and changes to gene expression in the pinfish (Lagodon rhomboides). Front. Ecol. Evol. 9, https://doi.org/10.3389/fevo.2021.809375 (2022).

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

35. Breitburg, D. et al. Declining oxygen in the global ocean and coastal waters. Science. 359, eaam7240 (2018).

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

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

38. Deutsch, C., Penn, J. L. & Seibel, B. Metabolic trait diversity shapes marine biogeography. Nature 585, 557–562 (2020).

39. Penn, J. L., Deutsch, C., Payne, J. L. & Sperling, E. A. Temperature-dependent hypoxia explains biogeography and severity of end-Permian marine mass extinction. Science 362, eaat1327 (2018).

40. Stockey, R. G., Pohl, A., Ridgwell, A., Finnegan, S. & Sperling, E. A. Decreasing Phanerozoic extinction intensity as a consequence of Earth surface oxygenation and metazoan ecophysiology. Proc. Natl. Acad. Sci. 118, e2101900118 (2021).

41. Taylor, K. E., Stouffer, R. J. & Meehl, G. A. An overview of CMIP5 and the experiment design. Bull. Am. Meteorol. Soc. 93, 485–498 (2012).

42. Eyring, V. et al. Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization. Geosci. Model Dev. 9, 1937–1958 (2016).

43. Dunne, J. P. et al. The GFDL Earth system model Version 4.1 (GFDL-ESM 4.1): overall coupled model description and simulation characteristics. J. Adv. Model. Earth Syst. 12, e2019MS002015 (2020).

44. Dunne, J. P. et al. GFDL's ESM2 global coupled climate-carbon ocean system models. Part I: physical formulation and baseline simulation characteristics. J. Clim. 25, 6646–6665 (2012).

45. Sellar, A. A. et al. UKESM1: description and evaluation of the U.K. Earth system model. J. Adv. Model. Earth Syst. 11, 4513–4558 (2019).

46. Hajima, T. et al. Development of the MIROC-ES2L Earth system model and the evaluation of biogeochemical processes and feedbacks. Geosci. Model Dev. 13, 2197–2244 (2020).

47. Séférían, R. et al. Evaluation of CNRM Earth system model, CNRM-ESM2-1: role of Earth system processes in present-day and future climate. J. Adv. Model. Earth Syst. 11, 4182–4227 (2019).

48. Ziehn, T. et al. The Australian Earth system model: ACCESS-ESM1.5. J. South. Hemisph. Earth Syst. Sci. 70, 193 (2020).

49. Seland, Ø. et al. Overview of the Norwegian Earth System Model (NorESM2) and key climate response of CMIP6 DECK, historical, and scenario simulations. Geosci. Model Dev. 13, 6165–6200 (2020).

50. Tjiputra, J. F. et al. Ocean biogeochemistry in the Norwegian Earth System Model version 2 (NorESM2). Geosci. Model Dev. 13, 2393–2431 (2020).

51. Deutsch, C. et al. Impact of warming on aquatic body sizes explained by metabolic scaling from microbes to macrofauna. Proc. Natl. Acad. Sci. 119, e2201345119 (2022).

52. Penn, J. L. & Deutsch, C. Avoiding ocean mass extinction from climate warming. Science 376, 524–526 (2022).

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

54. Frölicher, T. L., Joos, F., Plattner, G.-K., Steinacher, M. & Doney, S. C. Natural variability and anthropogenic trends in oceanic oxygen in a coupled carbon cycle-climate model ensemble. Glob. Biogeochem. Cycles 23, 1–15 (2009).

55. Drijfhout, S., Van Oldenborgh, G. J. & Cimatoribus, A. Is a Decline of AMOC Causing the Warming Hole above the North Atlantic in Observed and Modeled Warming Patterns? J. Clim. 25, 8373–8379 (2012).

56. Menary, M. B. & Wood, R. A. An anatomy of the projected North Atlantic warming hole in CMIP5 models. Clim. Dyn. 50, 3063–3080 (2018).

57. Swingedouw, D. et al. On the risk of abrupt changes in the North Atlantic subpolar gyre in CMIP6 models. Ann. N. Y. Acad. Sci. 1504, 187–201 (2021).

58. Schwinger, J., Asaadi, A., Goris, N. & Lee, H. Possibility for strong northern hemisphere high-latitude cooling under negative emissions. Nat. Commun. 13, 1095 (2022).

59. Ehlert, D. & Zickfeld, K. What determines the warming commitment after cessation of CO₂ emissions? Environ. Res. Lett. 12, 015002 (2017).

60. Williams, R. G., Roussenov, V., Frölicher, T. L. & Goodwin, P. Drivers of Continued Surface Warming After Cessation of Carbon Emissions. Geophys. Res. Lett. 44, 10,633–10,642 (2017).

61. MacDougall, A. H. et al. Is there warming in the pipeline? A multi-model analysis of the Zero Emissions Commitment from CO₂. Biogeosciences 17, 2987–3016 (2020).

62. Li, X., Zickfeld, K., Mathesius, S., Kohfeld, K. & Matthews, J. B. R. Irreversibility of marine climate change impacts under carbon dioxide removal. Geophys. Res. Lett. 47, e2020GL088507 (2020).

63. Garcia, H. E. & Gordon, L. I. Oxygen solubility in seawater: better fitting equations. Limnol. Oceanogr. 37, 1307–1312 (1992).

64. Salvanes, A. G. V. & Gibbons, M. J. Adaptation to hypoxic environments; bearded gobies Sufflogobius bibarbatus in the Benguela upwelling ecosystem. J. Fish Biol. 92, 752–772 (2018).

65. Meyer, A. L. S., Bentley, J., Odoulami, R. C., Pigot, A. L. & Trisos, C. H. Risks to biodiversity from temperature overshoot pathways. Philos. Trans. R. Soc. B Biol. Sci. 377, 20210394 (2022).

66. Mathesius, S., Hofmann, M., Caldeira, K. & Schellnhuber, H. J. Long-term response of oceans to CO₂ removal from the atmosphere. Nat. Clim. Change 5, 1107–1113 (2015).

67. Garcia, H. E. & Gordon, L. I. Oxygen solubility in seawater: better fitting equations. Limnol. Oceanogr. 37, 1307–1312 (1992).

68. Gillett, N. P., Arora, V. K., Zickfeld, K., Marshall, S. J. & Merryfield, W. J. Ongoing climate change following a complete cessation of carbon dioxide emissions. Nat. Geosci. 4, 83–87 (2011).

69. Oschlies, A. A committed fourfold increase in ocean oxygen loss. Nat. Commun. 12, 2307 (2021).

70. Frölicher, T. L. et al. Contrasting Upper and Deep Ocean Oxygen Response to Protracted Global Warming. Global Biogeochem. Cycles 34, e2020GB006601 (2020).

71. Yamamoto, A. et al. Global deep ocean oxygenation by enhanced ventilation in the Southern Ocean under long-term global warming. Glob. Biogeochem. Cycles 29, 1801–1815 (2015).

72. Frölicher, T. L., Rodgers, K. B., Stock, C. A. & Cheung, W. W. L. Sources of uncertainties in 21st century projections of potential ocean ecosystem stressors. Glob. Biogeochem. Cycles 30, 1224–1243 (2016).

73. Sutton, T. T. Vertical ecology of the pelagic ocean: classical patterns and new perspectives. J. Fish Biol. 83, 1508–1527 (2013).

74. Bauer, N. et al. Exploring risks and benefits of overshooting a 1.5 °C carbon budget over space and time. Environ. Res. Lett. 18, 054015 (2023).
75. Clarke, T. M. et al. Aerobic growth index (AGI): an index to understand the impacts of ocean warming and deoxygenation on global marine fisheries resources. *Prog. Oceanogr.* **195**, 102588 (2021).

76. Morée, A., Clarke, T. M., Cheung, W. W. L. & Frölicher, T. L. Impact of deoxygenation and warming on global marine species in the 21st century. *Biogeosciences* **20**, 2425–2454 (2023).

77. Meehl, G. A. et al. Climate Model Intercomparisons: preparing for the next phase. *Eos, Trans. Am. Geophys. Union* **95**, 77–78 (2014).

78. Kwiatkowski, L. et al. Twenty-first century ocean warming, acidification, deoxygenation, and upper-ocean nutrient and primary production decline from CMIP6 model projections. *Biogeosciences* **17**, 3439–3470 (2020).

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

80. Clarke, A. & Fraser, K. P. Why does metabolism scale with temperature? *Funct. Ecol.* **18**, 243–251 (2004).

81. Seebacher, F., White, C. R. & Franklin, C. E. Physiological plasticity increases resilience of ectothermic animals to climate change. *Nat. Clim. Change* **5**, 61–66 (2015).

82. Piiper, J., Dejours, P., Haab, P. & Rahn, H. Concepts and basic quantities in gas exchange physiology. *Respir. Physiol.* **13**, 292–304 (1971).

83. Gillooly, J. F., Brown, J. H., West, G. B., Savage, V. M. & Charnov, E. L. Effects of size and temperature on metabolic rate. *Science* **293**, 2248–2251 (2001).

84. Malte, H. & Weber, R. E. A mathematical model for gas exchange in the fish gill based on non-linear blood gas equilibrium curves. *Respir. Physiol.* **62**, 359–374 (1985).

## Acknowledgements

This work was supported by the European Union's Horizon 2020 research and innovation programme with the COMFORT project under the grant agreement No 820989 (Y.S.F., R.S., J.T., J.S. and T.L.F.), the ESM2025 project under the grant agreement No 101003536 (R.S.), the OceanNETs project under the grant agreement No 869357 (J.S.), and the PROVIDE project under the grant agreement No 101003687 (T.L.F.). C.D.J. was supported by the Joint UK BEIS/Defra Met Office Hadley Centre Climate Programme (GA01101). T.Z. received funding from the Australian Government under the National Environmental Science Programme. J.S. was supported by the Research Council of Norway (project IMPOSE with grant no 294930). John Dunne is acknowledged for his comments and suggestions on an internal review of the manuscript. We acknowledge the World Climate Research Programme, which, through its Working Group on Coupled Modelling, coordinated and promoted CMIP5 and CMIP6. The authors thank the climate modelling groups for producing and making available their model output, the Earth System Grid Federation (ESGF) for archiving the data and providing access, and the multiple funding agencies who support CMIP5, CMIP6 and ESGF.

## Author contributions

Y.S.F. and R.S. conceived the study, performed the analysis and wrote the manuscript with contributions from A.Y., A.L., C.D.J., F.A.B., J.G.J., J.T., J.S., M.K., T.L.F. and T.Z.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-023-01002-1.

**Correspondence** and requests for materials should be addressed to Yeray Falcón.

**Peer review information** *Communications Earth & Environment* thanks Andreas Oschlies and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Christopher Cornwall, Alienor Lavergne, Clare Davis. A peer review file is available.

**Reprints and permission information** is available at http://www.nature.com/reprints

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

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

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

© The Author(s) 2023