# Emergent constraint on oxygenation of the upper South Eastern Pacific oxygen minimum zone in the twenty-first century

Ivan Almendra, Boris Dewitte, Véronique Garçon, Praxedes Muñoz, Carolina Parada, Ivonne Montes, Olaf Duteil, Aurélien Paulmier, Oscar Pizarro, Marcel Ramos, et al.

► To cite this version:

Ivan Almendra, Boris Dewitte, Véronique Garçon, Praxedes Muñoz, Carolina Parada, et al. Emergent constraint on oxygenation of the upper South Eastern Pacific oxygen minimum zone in the twenty-first century. Communications Earth & Environment, 2024, 5 (1), pp.284. 10.1038/s43247-024-01427-2. insu-04600293

HAL Id: insu-04600293

https://insu.hal.science/insu-04600293

Submitted on 4 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Emergent constraint on oxygenation of the upper South Eastern Pacific oxygen minimum zone in the twenty-first century

Ivan Almendra, Boris Dewitte, Véronique Garçon, Praxedes Muñoz, Carolina Parada, Ivonne Montes, Olaf Duteil, Aurélien Paulmier, Oscar Pizarro, Marcel Ramos, et al.

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

► To cite this version:

Ivan Almendra, Boris Dewitte, Véronique Garçon, Praxedes Muñoz, Carolina Parada, et al. Emergent constraint on oxygenation of the upper South Eastern Pacific oxygen minimum zone in the twenty-first century. Communications Earth & Environment, 2024, 5 (1), pp.284. 10.1038/s43247-024-01427-2. insu-04600293

HAL Id: insu-04600293

https://insu.hal.science/insu-04600293v1

Submitted on 4 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Emergent constraint on oxygenation of the upper South Eastern Pacific oxygen minimum zone in the twenty-first century

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

Check for updates

Ivan Almendra Ⓔ<fnref n="1" />, Boris Dewitte²,³,⁴,⁵✉, Véronique Garçon Ⓔ⁶, Praxedes Muñoz²,³, Carolina Parada<fnref n="1" />, Ivonne Montes Ⓔ⁷, Olaf Duteil⁸, Aurélien Paulmier⁹, Oscar Pizarro<fnref n="1" />, Marcel Ramos²,³, Wolfgang Koeve Ⓔ⁸ & Andreas Oschlies Ⓔ⁸

As a consequence of on-going global warming, the ocean is losing oxygen, which has implications not only in terms of marine resources management and food supply but also in terms of the potentially important feedback on the global carbon cycle and climate. Of particular scrutiny are the extended zones of already low levels of oxygen called the oxygen minimum zones (OMZs) embedded in the subsurface waters of the productive Eastern Boundary Upwelling Systems (EBUS). These OMZs are currently diversely simulated by state-of-the-art Earth System Models (ESM) hampering a reliable projection of ocean deoxygenation on marine ecosystem services in these regions. Here we focus on the most emblematic EBUS OMZs of the planet, that of the South Eastern Pacific (SEP), which is under the direct influence of the El Niño Southern Oscillation (ENSO), the main climate mode on interannual timescales at global scale. We show that, despite the low consensus among ESM long-term projections of oxygen levels, the sensitivity of the depth of the upper margin (oxycline) of the SEP OMZ to El Niño events in an ensemble of ESMs can be used as a predictor of its long-term trend, which establishes an emergent constraint for the SEP OMZ. Because the oxycline along the coast of Peru and Chile deepens during El Niño events, the upper bound of the SEP OMZ is thus likely to deepen in the future climate, therefore oxygenating the SEP OMZ. This has implications not only for understanding the nitrogen and carbon cycles at global scale but also for designing adaptation strategies for regional upper-ocean ecosystem services.

Dissolved oxygen (DO) sustains marine life<fnref n="1" /> and constrains the habitat of almost all marine species². Yet, extended regions of the ocean where oxygen concentrations can fall beyond the detection limit of most sensors, the oxygen-deficient regions or Oxygen Minimum Zones (OMZ), are embedded in very productive regions with high biodiversity³. The largest OMZs of the planet are found in the Eastern Boundary Upwelling Systems of the Pacific Ocean. They are formed from sluggish circulation⁴,⁵ combined with an intense upwelling-induced organic material export⁶,⁷. OMZs regulate global biogeochemical cycles through producing potent greenhouse gases

like nitrous oxide and methane⁸,⁹. They also constrain the distributions and abundances of many zooplankton species and some mesopelagic species through forming a respiratory barrier¹⁰⁻¹². Current projections suggest that OMZs might expand under global warming (Stramma et al. ¹³,¹⁴), which has raised concern on the potential alteration of the structure and functioning of the food-web in these regions, with consequences on food security³. It has urged the scientific community to enhance the observing system of the yet poorly sampled essential variable oxygen¹⁵. In fact the scarcity of oxygen data, combined to the persistent difficulty of global Earth System models

<footnote n="1">¹Departamento de Geofísica, Universidad de Concepción, Concepción, Chile. ²Núcleo Milenio de Ecología y Manejo Sustentable (ESMOI), Facultad de Ciencias del Mar, Departamento de Biología Marina, Universidad Católica del Norte, Coquimbo, Chile. ³Departamento de Biología, Facultad de Ciencias del Mar, Universidad Católica del Norte, Coquimbo, Chile. ⁴Centro de Estudios Avanzado en Zonas Áridas (CEAZA), La Serena, Chile. ⁵CECI, Université de Toulouse, CERFACS/CNRS, Toulouse, France. ⁶Institut de Physique du Globe de Paris, IPGP/CNRS, Paris, France. ⁷Instituto Geofísico del Peru (IGP), Lima, Peru. ⁸GEOMAR – Helmholtz Centre for Ocean Research Kiel, Düstembrooker Weg. 20, 24103 Kiel, Germany. ⁹Laboratoire d'Etudes en Géophysique et Océanographie Spatiales, Université de Toulouse, LEGOS(CNES/CNRS/IRD/UPS), Toulouse, France. ✉ e-mail: boris.dewitte@ceaza.cl</footnote>
<figure><img src="image_3.png" /></figure>

# Emergent constraint on oxygenation of the upper South Eastern Pacific oxygen minimum zone in the twenty-first century

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

Check for updates

Ivan Almendra Ⓔ<fnref n="1" />, Boris Dewitte²,³,⁴,⁵✉, Véronique Garçon Ⓔ⁶, Praxedes Muñoz²,³, Carolina Parada<fnref n="1" />, Ivonne Montes Ⓔ⁷, Olaf Duteil⁸, Aurélien Paulmier⁹, Oscar Pizarro<fnref n="1" />, Marcel Ramos²,³, Wolfgang Koeve Ⓔ⁸ & Andreas Oschlies Ⓔ⁸

As a consequence of on-going global warming, the ocean is losing oxygen, which has implications not only in terms of marine resources management and food supply but also in terms of the potentially important feedback on the global carbon cycle and climate. Of particular scrutiny are the extended zones of already low levels of oxygen called the oxygen minimum zones (OMZs) embedded in the subsurface waters of the productive Eastern Boundary Upwelling Systems (EBUS). These OMZs are currently diversely simulated by state-of-the-art Earth System Models (ESM) hampering a reliable projection of ocean deoxygenation on marine ecosystem services in these regions. Here we focus on the most emblematic EBUS OMZs of the planet, that of the South Eastern Pacific (SEP), which is under the direct influence of the El Niño Southern Oscillation (ENSO), the main climate mode on interannual timescales at global scale. We show that, despite the low consensus among ESM long-term projections of oxygen levels, the sensitivity of the depth of the upper margin (oxycline) of the SEP OMZ to El Niño events in an ensemble of ESMs can be used as a predictor of its long-term trend, which establishes an emergent constraint for the SEP OMZ. Because the oxycline along the coast of Peru and Chile deepens during El Niño events, the upper bound of the SEP OMZ is thus likely to deepen in the future climate, therefore oxygenating the SEP OMZ. This has implications not only for understanding the nitrogen and carbon cycles at global scale but also for designing adaptation strategies for regional upper-ocean ecosystem services.

Dissolved oxygen (DO) sustains marine life<fnref n="1" /> and constrains the habitat of almost all marine species². Yet, extended regions of the ocean where oxygen concentrations can fall beyond the detection limit of most sensors, the oxygen-deficient regions or Oxygen Minimum Zones (OMZ), are embedded in very productive regions with high biodiversity³. The largest OMZs of the planet are found in the Eastern Boundary Upwelling Systems of the Pacific Ocean. They are formed from sluggish circulation⁴,⁵ combined with an intense upwelling-induced organic material export⁶,⁷. OMZs regulate global biogeochemical cycles through producing potent greenhouse gases

like nitrous oxide and methane⁸,⁹. They also constrain the distributions and abundances of many zooplankton species and some mesopelagic species through forming a respiratory barrier¹⁰⁻¹². Current projections suggest that OMZs might expand under global warming (Stramma et al. ¹³,¹⁴), which has raised concern on the potential alteration of the structure and functioning of the food-web in these regions, with consequences on food security³. It has urged the scientific community to enhance the observing system of the yet poorly sampled essential variable oxygen¹⁵. In fact the scarcity of oxygen data, combined to the persistent difficulty of global Earth System models

<footnote n="1">¹Departamento de Geofísica, Universidad de Concepción, Concepción, Chile. ²Núcleo Milenio de Ecología y Manejo Sustentable (ESMOI), Facultad de Ciencias del Mar, Departamento de Biología Marina, Universidad Católica del Norte, Coquimbo, Chile. ³Departamento de Biología, Facultad de Ciencias del Mar, Universidad Católica del Norte, Coquimbo, Chile. ⁴Centro de Estudios Avanzado en Zonas Áridas (CEAZA), La Serena, Chile. ⁵CECI, Université de Toulouse, CERFACS/CNRS, Toulouse, France. ⁶Institut de Physique du Globe de Paris, IPGP/CNRS, Paris, France. ⁷Instituto Geofísico del Peru (IGP), Lima, Peru. ⁸GEOMAR – Helmholtz Centre for Ocean Research Kiel, Düstembrooker Weg. 20, 24103 Kiel, Germany. ⁹Laboratoire d'Etudes en Géophysique et Océanographie Spatiales, Université de Toulouse, LEGOS(CNES/CNRS/IRD/UPS), Toulouse, France. ✉ e-mail: boris.dewitte@ceaza.cl</footnote>
(ESMs) to realistically simulate OMZs<sup>16,17</sup> (Supplementary Fig. S1), has hampered a proper understanding of the processes driving oxygen variability at a variety of timescales. While the solubility effect dominantly controls global warming driven oxygen variability near the surface, a quantitative understanding of contributions from other mechanisms is still lacking<sup>18,19</sup>.

In the Pacific ocean, climate models still poorly simulate certain aspects of the subthermocline circulation, including the Intermediate Deep Water (IDW), the Equatorial Intermediate Current System<sup>20</sup>, and the Antarctic Intermediate Water<sup>16,21,22</sup> that modulate the OMZ waters at interannual to decadal frequencies<sup>23</sup>. Within the thermocline and above, details in the equatorial eastern Pacific circulation not well accounted for by climate models influence the transport of deoxygenated waters to the OMZs. These include the narrow eastward Tsuchiya jets located a few degrees on either side of the equator<sup>24,25</sup> that serve as a conduit by which tropical variability can modulate the OMZs at seasonal to decadal timescales<sup>20,26</sup>. Models have thus exhibited a large spread in OMZ properties<sup>16</sup> (See Supplementary Fig. S1) and their projected future climate response (see Fig. 1a for the SEP OMZ volume). Such a spread on climatic timescales questions the extent to which it correlates to variations on shorter timescales, which could provide an emergent constraint<sup>27</sup> to pinpoint observational predictors for oxygen changes and in the meantime guide the design of observational network for this essential variable<sup>15</sup>.

In the SEP, potentially observable OMZ metrics include the oxycline depth along the coast at a few sites where long-term monitoring has been performed<sup>28-30</sup>, although these coastal sites mostly sample the water column above the OMZ and may monitor anthropogenic signals from different sources difficult to interpret. Regarding relevant timescales for oxygen variability in SEP, it appears obvious to consider the interannual timescale as the El Niño Southern Oscillation (ENSO) is the main mode of interannual variability at global scale. On the one hand, the strong ENSO oceanic

teleconnection along the west coast of South America through the propa-
gation of planetary waves<sup>31-36</sup> (among others) implies a large signal-to-noise
ratio in interannual fluctuations of oxygen due to wave-induced advection
processes. On the other hand, the OMZ sensitivity to ENSO can operate
through many other mechanisms<sup>19</sup>. For instance, while El Niño-induced
deepening of the thermocline along the coast of Peru and Chile<sup>30,33,35</sup> may
foster oxygenation through advection and diffusive processes, DO transport
by the equatorial subsurface currents (i.e., the Equatorial UnderCurrent
(EUC) and the two Tsuchiya Jets (hereafter pSSCC and sSSCC for primary
and secondary subsurface undercurrents respectively)) into the Peru-Chile
Undercurrent (PCUC) may also be modulated by ENSO through altering
their strength<sup>26</sup> and/or the DO concentration in the equatorial region<sup>37,38</sup>.
Montes et al.<sup>26</sup> indicate in particular that the PCUC is mainly fed by the
pSSCC and the EUC during La Niña, while the sSSCC is a main contributor
of the PCUC waters during El Niño. Since these currents are important
contributors to the oxygen balance off Peru<sup>26</sup>, their modulation by ENSO
could impact the ventilation of the OMZ. A slowdown of the EUC is in
general conducive of a more intense OMZ<sup>39,40</sup>. However, the SEP OMZ is
very sensitive to the details of these currents (vertical structure, intensity,
position)<sup>26,39</sup>. Thermocline deepening during El Niño event, by reducing the
nutrient supply to the euphotic layer and associated primary production,
also reduces the downward flux of organic matter, which decreases oxygen
consumption by respiration<sup>41,42</sup>. Enhanced eddy activity during El Niño off
Peru<sup>43</sup> also tends to ventilate the nearshore OMZ<sup>44</sup>. Altogether, these
mechanisms yield a positive relationship between the DO levels in the upper
OMZ (oxycline depth) and ENSO<sup>29,37,38,44-46</sup>. At longer timescales (multi-
decadal), the relationship between the SEP OMZ and tropical variability is
more ubiquitous since it is mediated by the meridional DO transport by the
subtropical cells connecting the subtropics to the tropics in the upper
thermocline. Duteil et al.<sup>47</sup> found in particular that the positive phase of the
PDO is associated with an expansion of the suboxic regions of the eastern

<figure><img src="image_4.png" /><figcaption>**Fig. 1 | Evolution of OMZ volume anomalies and ENSO variance for the period 1850–2100.** a Evolution of the OMZ volume anomalies (in m³) in the models. The OMZ volume is defined as the volume of waters with oxygen concentration below 45 µmol L⁻¹ over the region (5°N–35°S; 160°E–70°W; 0–1000 m). Anomalies are relative to the mean over the period 1950–2005 shown for each model in Supplementary Fig. S1. Colors are for each individual model (see color code on the right-hand side panel, see also Table S1 for model details); the black curve is for the</figcaption></figure>

ensemble mean, while the shading in purple stands for ± the standard deviation
amongst the ensemble. 16 models are considered out of the 19 models analyzed in
this study because three of them have a OMZ volume near zero with the 45 µmol L⁻¹
threshold (see model names in red in Table S1). **b** 21-year running variance of the E
index (see Methods) (December–January–February mean) for the ensemble mean.
The purple stands for ± the standard deviation amongst models. The horizontal
dashed line in red corresponds to the mean value.
(ESMs) to realistically simulate OMZs<sup>16,17</sup> (Supplementary Fig. S1), has hampered a proper understanding of the processes driving oxygen variability at a variety of timescales. While the solubility effect dominantly controls global warming driven oxygen variability near the surface, a quantitative understanding of contributions from other mechanisms is still lacking<sup>18,19</sup>.

In the Pacific ocean, climate models still poorly simulate certain aspects of the subthermocline circulation, including the Intermediate Deep Water (IDW), the Equatorial Intermediate Current System<sup>20</sup>, and the Antarctic Intermediate Water<sup>16,21,22</sup> that modulate the OMZ waters at interannual to decadal frequencies<sup>23</sup>. Within the thermocline and above, details in the equatorial eastern Pacific circulation not well accounted for by climate models influence the transport of deoxygenated waters to the OMZs. These include the narrow eastward Tsuchiya jets located a few degrees on either side of the equator<sup>24,25</sup> that serve as a conduit by which tropical variability can modulate the OMZs at seasonal to decadal timescales<sup>20,26</sup>. Models have thus exhibited a large spread in OMZ properties<sup>16</sup> (See Supplementary Fig. S1) and their projected future climate response (see Fig. 1a for the SEP OMZ volume). Such a spread on climatic timescales questions the extent to which it correlates to variations on shorter timescales, which could provide an emergent constraint<sup>27</sup> to pinpoint observational predictors for oxygen changes and in the meantime guide the design of observational network for this essential variable<sup>15</sup>.

In the SEP, potentially observable OMZ metrics include the oxycline depth along the coast at a few sites where long-term monitoring has been performed<sup>28-30</sup>, although these coastal sites mostly sample the water column above the OMZ and may monitor anthropogenic signals from different sources difficult to interpret. Regarding relevant timescales for oxygen variability in SEP, it appears obvious to consider the interannual timescale as the El Niño Southern Oscillation (ENSO) is the main mode of interannual variability at global scale. On the one hand, the strong ENSO oceanic

teleconnection along the west coast of South America through the propa-
gation of planetary waves<sup>31-36</sup> (among others) implies a large signal-to-noise
ratio in interannual fluctuations of oxygen due to wave-induced advection
processes. On the other hand, the OMZ sensitivity to ENSO can operate
through many other mechanisms<sup>19</sup>. For instance, while El Niño-induced
deepening of the thermocline along the coast of Peru and Chile<sup>30,33,35</sup> may
foster oxygenation through advection and diffusive processes, DO transport
by the equatorial subsurface currents (i.e., the Equatorial UnderCurrent
(EUC) and the two Tsuchiya Jets (hereafter pSSCC and sSSCC for primary
and secondary subsurface undercurrents respectively)) into the Peru-Chile
Undercurrent (PCUC) may also be modulated by ENSO through altering
their strength<sup>26</sup> and/or the DO concentration in the equatorial region<sup>37,38</sup>.
Montes et al.<sup>26</sup> indicate in particular that the PCUC is mainly fed by the
pSSCC and the EUC during La Niña, while the sSSCC is a main contributor
of the PCUC waters during El Niño. Since these currents are important
contributors to the oxygen balance off Peru<sup>26</sup>, their modulation by ENSO
could impact the ventilation of the OMZ. A slowdown of the EUC is in
general conducive of a more intense OMZ<sup>39,40</sup>. However, the SEP OMZ is
very sensitive to the details of these currents (vertical structure, intensity,
position)<sup>26,39</sup>. Thermocline deepening during El Niño event, by reducing the
nutrient supply to the euphotic layer and associated primary production,
also reduces the downward flux of organic matter, which decreases oxygen
consumption by respiration<sup>41,42</sup>. Enhanced eddy activity during El Niño off
Peru<sup>43</sup> also tends to ventilate the nearshore OMZ<sup>44</sup>. Altogether, these
mechanisms yield a positive relationship between the DO levels in the upper
OMZ (oxycline depth) and ENSO<sup>29,37,38,44-46</sup>. At longer timescales (multi-
decadal), the relationship between the SEP OMZ and tropical variability is
more ubiquitous since it is mediated by the meridional DO transport by the
subtropical cells connecting the subtropics to the tropics in the upper
thermocline. Duteil et al.<sup>47</sup> found in particular that the positive phase of the
PDO is associated with an expansion of the suboxic regions of the eastern

<figure><img src="image_5.png" /><figcaption>**Fig. 1 | Evolution of OMZ volume anomalies and ENSO variance for the period 1850–2100.** a Evolution of the OMZ volume anomalies (in m³) in the models. The OMZ volume is defined as the volume of waters with oxygen concentration below 45 µmol L⁻¹ over the region (5°N–35°S; 160°E–70°W; 0–1000 m). Anomalies are relative to the mean over the period 1950–2005 shown for each model in Supplementary Fig. S1. Colors are for each individual model (see color code on the right-hand side panel, see also Table S1 for model details); the black curve is for the</figcaption></figure>

ensemble mean, while the shading in purple stands for ± the standard deviation
amongst the ensemble. 16 models are considered out of the 19 models analyzed in
this study because three of them have a OMZ volume near zero with the 45 µmol L⁻¹
threshold (see model names in red in Table S1). **b** 21-year running variance of the E
index (see Methods) (December–January–February mean) for the ensemble mean.
The purple stands for ± the standard deviation amongst models. The horizontal
dashed line in red corresponds to the mean value.
tropical Pacific OMZ despite the thermocline deepening along the coast of Peru and Chile<sup>48</sup>. In fact, current global ocean models still have limitations in reproducing measured changes in concentrations and their distribution of low oxygen waters over the past 50 years<sup>14,16,49</sup> and current understanding of decadal oxygen variability remains limited<sup>50</sup>.

Recent studies have also pointed out the existence of distinct ENSO regimes that manifest as different peak Sea Surface Temperature (SST) anomaly patterns along the equatorial Pacific as well as different magnitudes<sup>51</sup>. In particular, strong El Niño events that are projected to increase in amplitude and frequency in the future climate<sup>52,53</sup> have their center of action in the Eastern Pacific yielding a marked and rapid deepening of the thermocline along the coast of Peru and Chile through geostrophic adjustments<sup>33-35</sup>, as well as significant fluctuations in the turbulent flow along the coast<sup>43</sup>, which has the potential to modulate the SEP OMZ<sup>29</sup>.

# Results

## Emergent constraint on oxygenation

Here, considering the current uncertainty in the future SEP upper OMZ trend, we use projections from the Coupled Climate Carbon Cycle Model Intercomparison Project (CMIP, Phases 5 and 6) to look for an emergent constraint linking the sensitivity of the upper OMZ in the SEP along the coast of Peru and Chile to El Niño events. The focus is on strong El Niño events that peak in the Eastern Pacific<sup>54</sup>, hereafter referred to as Eastern Pacific El Niño (EP EN), considering their large imprint on the OMZ<sup>29,45</sup>. Figure 2 represents the values of the regression coefficient between an EP EN index (E index, see method section for definition) and simulated oxygen concentration for an ensemble of ESMs from CMIP5 and CMIP6 (see Method section) in a map at 100 m. Vertical zonal sections at 12°S and 30°S, two key upwelling centers along the coast of Peru and Chile, are provided in the Supplementary Material (Fig. S2b, c). Positive values of the regression coefficients indicate increasing oxygen concentrations during EN, which is almost everywhere along the upper margin of the OMZ, that is in the vicinity of the oxycline. At 100 m depth, the values of regression coefficients decline southward and westward. The maximum values are found near the coast above the mean 45 µmol L<sup>-1</sup> isopleths (rectangles in Fig. S2b, c). They also correspond to the regions with positive long-term trends in oxygen (Fig. 2b, S3b, c) with, however, a low consensus amongst the models (cf. Fig. S4) except at 12°S.

Given the emergent constraint framework<sup>55-59</sup>, we speculate that ENSO-induced upper-ocean oxygen fluctuations in these regions can be used as predictors of long-term changes in oxygen. This is grounded in the expectation that, collectively, processes driving fluctuations in oxygen in the upper OMZ operate similarly on interannual timescales (i.e., ENSO

timescales) and on longer timescales of climate change, which is arguably debatable<sup>60</sup>. The by now widely used statistical approach of emergent constraints has, however, value for narrowing uncertainties in the projections of climate change. Overall Fig. 2 indicates that the oxygen climate change pattern within the OMZ at the oxycline level is comparable to that of its sensitivity to EN, that is a tendency for oxygenation off the coasts of Peru and Chile that extends off-shore further off Peru than off Chile. There is a low consensus amongst the models in terms of the trend over most of the OMZ extent at 100 m (white stippling in Fig. 2b) except off central Peru in a localized area. The regions of maximum positive trend correspond to those where the dispersion amongst models is largest (Supplementary Fig. S4). This large uncertainty associated with future oxygen changes across ESM ensembles<sup>61,62</sup> motivates the determination of an emergent relationship.

The Figure 3a, b displays the scatter plot of the climate sensitivity against the ENSO sensitivity (E index) in terms of oxygen averaged along the zonal sections at 12°S and 30°S in a coastal box (see rectangles in Fig. S2b, c) for the model ensemble. The linear fit shows a significant positive value of the correlation (c = 0.59 at 12°S and c = 0.62 at 30°S), indicating that there is a linear relationship between the sensitivity of oxygen to ENSO-driven SST variability and the long-term sensitivity of oxygen to climate change in the twenty-first century in the two main upwelling centers along the coasts of South America. Such a relationship holds over most of the coast of Peru and northern-central Chile (Fig. 3). However, the latitudinal range of significant correlation depends on where the peak SST warming during EN takes place. For the E index, the correlation is significant between 5°S and 17°S and between 26°S and 30°S (blue curve in Fig. 3c), while for the C index (see Methods), it is significant in the latitudinal range between 17°S and 27°S (green curve in Fig. 3c). This is interpreted as resulting from the seasonal evolution of strong EN events. In particular, during the development of strong EN, SST warming takes place in the central equatorial Pacific during Austral summer (C region) and, the SST anomalies of the concurrent EN subsequently peak in the eastern equatorial Pacific (E region) the next year in Austral summer. The yearly mean E and C indices thus can both grasp variability associated with the same strong EN. However, by construction, the C index also grasps variability associated with moderate CP EN and La Niña events<sup>54</sup>, which has distinct teleconnections along the west coast of South America<sup>36,63</sup>. This explains why the relationship between ENSO and climate sensitivities can vary as a function of latitude according to the E (blue curve in Fig. 3c) and C (green curve in Fig. 3c) indices.

Overall, however, the results indicate that EN events, whether it is through the influence of the C or E index, are associated with oxygenation of the upper SEP OMZ and that the spread amongst models is correlated to that of the oxygen climate sensitivity from the equatorial region up to ~31°S

**Fig. 2 |** Comparison between the ENSO-oxygen relationship and the long-term trend at 100 m. Ensemble mean amongst the models of the a regression coefficient between dissolved oxygen (DO) and the E index (in µmol L⁻¹ °C⁻¹) over the period 1920–2014 and of the b long-term linear trend of DO over the period 2015–2100 (in µmol L⁻¹ per decade). The contours in blue and red indicate the limit of the OMZ (i.e., 45 µmol L⁻¹ isopleths) at 300 m for the present (1920–2014) and future (2015–2100) climates, respectively. The zonal lines at 12°S and 30°S indicate the latitudes of the vertical sections for the same quantities shown in Supplementary Fig. S2b, c. Black stippling in a and white stippling in b indicates grid points where <80% of the models agrees on the sign.

<figure><img src="image_5.png" /></figure>
tropical Pacific OMZ despite the thermocline deepening along the coast of Peru and Chile<sup>48</sup>. In fact, current global ocean models still have limitations in reproducing measured changes in concentrations and their distribution of low oxygen waters over the past 50 years<sup>14,16,49</sup> and current understanding of decadal oxygen variability remains limited<sup>50</sup>.

Recent studies have also pointed out the existence of distinct ENSO regimes that manifest as different peak Sea Surface Temperature (SST) anomaly patterns along the equatorial Pacific as well as different magnitudes<sup>51</sup>. In particular, strong El Niño events that are projected to increase in amplitude and frequency in the future climate<sup>52,53</sup> have their center of action in the Eastern Pacific yielding a marked and rapid deepening of the thermocline along the coast of Peru and Chile through geostrophic adjustments<sup>33-35</sup>, as well as significant fluctuations in the turbulent flow along the coast<sup>43</sup>, which has the potential to modulate the SEP OMZ<sup>29</sup>.

# Results

## Emergent constraint on oxygenation

Here, considering the current uncertainty in the future SEP upper OMZ trend, we use projections from the Coupled Climate Carbon Cycle Model Intercomparison Project (CMIP, Phases 5 and 6) to look for an emergent constraint linking the sensitivity of the upper OMZ in the SEP along the coast of Peru and Chile to El Niño events. The focus is on strong El Niño events that peak in the Eastern Pacific<sup>54</sup>, hereafter referred to as Eastern Pacific El Niño (EP EN), considering their large imprint on the OMZ<sup>29,45</sup>. Figure 2 represents the values of the regression coefficient between an EP EN index (E index, see method section for definition) and simulated oxygen concentration for an ensemble of ESMs from CMIP5 and CMIP6 (see Method section) in a map at 100 m. Vertical zonal sections at 12°S and 30°S, two key upwelling centers along the coast of Peru and Chile, are provided in the Supplementary Material (Fig. S2b, c). Positive values of the regression coefficients indicate increasing oxygen concentrations during EN, which is almost everywhere along the upper margin of the OMZ, that is in the vicinity of the oxycline. At 100 m depth, the values of regression coefficients decline southward and westward. The maximum values are found near the coast above the mean 45 µmol L<sup>-1</sup> isopleths (rectangles in Fig. S2b, c). They also correspond to the regions with positive long-term trends in oxygen (Fig. 2b, S3b, c) with, however, a low consensus amongst the models (cf. Fig. S4) except at 12°S.

Given the emergent constraint framework<sup>55-59</sup>, we speculate that ENSO-induced upper-ocean oxygen fluctuations in these regions can be used as predictors of long-term changes in oxygen. This is grounded in the expectation that, collectively, processes driving fluctuations in oxygen in the upper OMZ operate similarly on interannual timescales (i.e., ENSO

timescales) and on longer timescales of climate change, which is arguably debatable<sup>60</sup>. The by now widely used statistical approach of emergent constraints has, however, value for narrowing uncertainties in the projections of climate change. Overall Fig. 2 indicates that the oxygen climate change pattern within the OMZ at the oxycline level is comparable to that of its sensitivity to EN, that is a tendency for oxygenation off the coasts of Peru and Chile that extends off-shore further off Peru than off Chile. There is a low consensus amongst the models in terms of the trend over most of the OMZ extent at 100 m (white stippling in Fig. 2b) except off central Peru in a localized area. The regions of maximum positive trend correspond to those where the dispersion amongst models is largest (Supplementary Fig. S4). This large uncertainty associated with future oxygen changes across ESM ensembles<sup>61,62</sup> motivates the determination of an emergent relationship.

The Figure 3a, b displays the scatter plot of the climate sensitivity against the ENSO sensitivity (E index) in terms of oxygen averaged along the zonal sections at 12°S and 30°S in a coastal box (see rectangles in Fig. S2b, c) for the model ensemble. The linear fit shows a significant positive value of the correlation (c = 0.59 at 12°S and c = 0.62 at 30°S), indicating that there is a linear relationship between the sensitivity of oxygen to ENSO-driven SST variability and the long-term sensitivity of oxygen to climate change in the twenty-first century in the two main upwelling centers along the coasts of South America. Such a relationship holds over most of the coast of Peru and northern-central Chile (Fig. 3). However, the latitudinal range of significant correlation depends on where the peak SST warming during EN takes place. For the E index, the correlation is significant between 5°S and 17°S and between 26°S and 30°S (blue curve in Fig. 3c), while for the C index (see Methods), it is significant in the latitudinal range between 17°S and 27°S (green curve in Fig. 3c). This is interpreted as resulting from the seasonal evolution of strong EN events. In particular, during the development of strong EN, SST warming takes place in the central equatorial Pacific during Austral summer (C region) and, the SST anomalies of the concurrent EN subsequently peak in the eastern equatorial Pacific (E region) the next year in Austral summer. The yearly mean E and C indices thus can both grasp variability associated with the same strong EN. However, by construction, the C index also grasps variability associated with moderate CP EN and La Niña events<sup>54</sup>, which has distinct teleconnections along the west coast of South America<sup>36,63</sup>. This explains why the relationship between ENSO and climate sensitivities can vary as a function of latitude according to the E (blue curve in Fig. 3c) and C (green curve in Fig. 3c) indices.

Overall, however, the results indicate that EN events, whether it is through the influence of the C or E index, are associated with oxygenation of the upper SEP OMZ and that the spread amongst models is correlated to that of the oxygen climate sensitivity from the equatorial region up to ~31°S

**Fig. 2 |** Comparison between the ENSO-oxygen relationship and the long-term trend at 100 m. Ensemble mean amongst the models of the a regression coefficient between dissolved oxygen (DO) and the E index (in µmol L⁻¹ °C⁻¹) over the period 1920–2014 and of the b long-term linear trend of DO over the period 2015–2100 (in µmol L⁻¹ per decade). The contours in blue and red indicate the limit of the OMZ (i.e., 45 µmol L⁻¹ isopleths) at 300 m for the present (1920–2014) and future (2015–2100) climates, respectively. The zonal lines at 12°S and 30°S indicate the latitudes of the vertical sections for the same quantities shown in Supplementary Fig. S2b, c. Black stippling in a and white stippling in b indicates grid points where <80% of the models agrees on the sign.

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

**Fig. 3** | Relationship between the sensitivity of dissolved oxygen (DO) to ENSO and to global warming (emergent constraint). The long-term (2015–2100) sensitivity of DO to SST warming in the Eastern equatorial Pacific (average over the E mode region) against the interannual sensitivity (1920–2014) of DO to the E index at a 12°S and b 30°S for the 18 ESMs (Mi dots, where *i* = 1,...,19). The straight line (black) in a and b corresponds to the linear fit to the Mi points. Pearson correlation value is indicated in the top left corner of each panel. The horizontal and vertical black segments on the point labels provide an estimate of the errors. The size of the black segment represents ± the standard deviation amongst 1000 estimates of the ENSO and climate sensitivities using a bootstrap method for which 80 years are chosen randomly. The colored segments on the point corresponding to CESM-LENS (M10) correspond to the errors associated with natural variability estimated as ± the standard deviation of the ENSO and climate sensitivities amongst the 34 members. The vertical line in blue with purple shading indicates the value of ENSO sensitivity of the regional biogeochemical model simulation (see Methods). The

shading stands for the 25% and 75% percentiles of the distribution of DO ENSO sensitivity obtained from a bootstrap method applied to the regional model data (1958–2008). 1000 values of the DO ENSO sensitivity were estimated based on 40 years of the simulation selected randomly amongst the 51 years. The red vertical line indicates an estimate of the observed value at 12°S. The value is estimated from Figure 6a of Graco et al.²⁹ taking the value at 100 m as an estimate of the mean between 25 m and 200 m like for the models. This value is further divided by 0.29 °C that corresponds to the norm of the observed E mode pattern (see Methods).

c Correlation amongst the models between the DO ENSO sensitivity and the DO climate sensitivity as a function of latitude between 5°S–39°S for the E index (blue) and the C index (green). DO is averaged in a region corresponding to the upper oxycline: a i.e., (25 m–200 m; 82°W–77°W), b (25 m–200 m; 75°W–72°W) (see rectangles in Supplementary Fig. S2b, c), and c 25 m–200 m in a coastal fringe of 3° in longitudinal extent. The thicker portions of the curves indicate where the correlations are significant at the 95% level based on a student's *t* test.

(the maximum correlation value between the two curves reaches 0.78 on average between 5°S and 31°S, Fig. 3c). Such an emergent relationship is also evidenced when considering a proxy of the oxycline depth anomalies (see Supplementary Fig. S5). In summary, the results indicate that models with higher interannual oxygen sensitivity to SST anomalies of the equatorial Pacific in historical simulations consistently project greater oxygenation within the oxycline (upper OMZ) per degree warming in the twenty-first century along the coast of Peru and Northern-central Chile.

## The forcing of the SEP OMZ by a changing ENSO variability

An implication of our results regards the mechanistic understanding of the sensitivity of the SEP upper OMZ to tropical forcing at low-frequency timescales. First, the relationship between the SEP OMZ and the tropical variability evidenced here is consistent with modeling studies focused on different timescales⁶⁴ and mean climate conditions⁶⁵⁻⁶⁸. From a statistical perspective, our results suggest that the oxygenation of the SEP upper OMZ during strong EN events may have a residual on the long-term mean because cool events (La Niña episode) and CP EN events that are weaker in magnitude (as measured as the amplitude of the C index, see Methods) than EP EN events⁵⁴ have a weaker imprint on the OMZ. In other words, the ENSO asymmetry is transmitted to the SEP upper OMZ through a number

of processes (e.g., radiation of downwelling extra-tropical Rossby waves, downwelling Kelvin wave propagation, enhanced mesoscale activity, etc. see Pitcher et al.¹⁹). This is evidenced by the map of skewness of dissolved oxygen anomalies at 100 m for the model simulating a realistic ENSO asymmetry (see Methods, Fig. 4a). It shows an elongated pattern of positive skewness extending from the equatorial region (~10°S) to along the coast of Peru and Chile. A region of peak positive skewness is also found along the sections at 12°S and 30°S that matches with regions of influence of ENSO (rectangles in Fig. S2b, c) but with a more off-shore extension along the oxycline. Retaining all the models yields comparable results albeit with a positive skewness peaking more off-shore below 100 m (Fig. S6).

As El Niño events are expected to increase in frequency and amplitude⁵²,⁶⁹ as well as to last longer⁷⁰,⁷¹ in the warmer climate, we may thus expect a long-term increase in oxygen concentration in the SEP upper OMZ as greenhouse gases in the atmosphere increase. This does not preclude decadal variations as the increase in ENSO variance is not steady⁷⁰,⁷² and other ventilation processes can take place at long timescales¹⁸ (Cabré et al.¹⁶,²⁰,⁷³,⁷⁴). Consistently, we find that our ensemble model evidences a significant long-term mean increase in the variance of the E index (Fig. 1b; the slope of the linear fit of the 21-year running variance of the E index reaches 0.04 °C century⁻¹ which is significant at the 95% confidence level based on a
<figure><img src="image_7.png" /></figure>

**Fig. 3 | Relationship between the sensitivity of dissolved oxygen (DO) to ENSO and to global warming (emergent constraint).** The long-term (2015–2100) sensitivity of DO to SST warming in the Eastern equatorial Pacific (average over the E mode region) against the interannual sensitivity (1920–2014) of DO to the E index at **a** 12°S and **b** 30°S for the 18 ESMs (Mi dots, where *i* = 1,...,19). The straight line (black) in **a** and **b** corresponds to the linear fit to the Mi points. Pearson correlation value is indicated in the top left corner of each panel. The horizontal and vertical black segments on the point labels provide an estimate of the errors. The size of the black segment represents ± the standard deviation amongst 1000 estimates of the ENSO and climate sensitivities using a bootstrap method for which 80 years are chosen randomly. The colored segments on the point corresponding to CESM-LENS (M10) correspond to the errors associated with natural variability estimated as ± the standard deviation of the ENSO and climate sensitivities amongst the 34 members. The vertical line in blue with purple shading indicates the value of ENSO sensitivity of the regional biogeochemical model simulation (see Methods). The

shading stands for the 25% and 75% percentiles of the distribution of DO ENSO sensitivity obtained from a bootstrap method applied to the regional model data (1958–2008). 1000 values of the DO ENSO sensitivity were estimated based on 40 years of the simulation selected randomly amongst the 51 years. The red vertical line indicates an estimate of the observed value at 12°S. The value is estimated from Figure 6a of Graco et al.²⁹ taking the value at 100 m as an estimate of the mean between 25 m and 200 m like for the models. This value is further divided by 0.29 °C that corresponds to the norm of the observed E mode pattern (see Methods).

**c** Correlation amongst the models between the DO ENSO sensitivity and the DO climate sensitivity as a function of latitude between 5°S–39°S for the E index (blue) and the C index (green). DO is averaged in a region corresponding to the upper oxycline: **a** i.e., (25 m–200 m; 82°W–77°W), **b** (25 m–200 m; 75°W–72°W) (see rectangles in Supplementary Fig. S2b, c), and **c** 25 m–200 m in a coastal fringe of 3° in longitudinal extent. The thicker portions of the curves indicate where the correlations are significant at the 95% level based on a student's *t* test.

(the maximum correlation value between the two curves reaches 0.78 on average between 5°S and 31°S, Fig. 3c). Such an emergent relationship is also evidenced when considering a proxy of the oxycline depth anomalies (see Supplementary Fig. S5). In summary, the results indicate that models with higher interannual oxygen sensitivity to SST anomalies of the equatorial Pacific in historical simulations consistently project greater oxygenation within the oxycline (upper OMZ) per degree warming in the twenty-first century along the coast of Peru and Northern-central Chile.

## The forcing of the SEP OMZ by a changing ENSO variability

An implication of our results regards the mechanistic understanding of the sensitivity of the SEP upper OMZ to tropical forcing at low-frequency timescales. First, the relationship between the SEP OMZ and the tropical variability evidenced here is consistent with modeling studies focused on different timescales⁶⁴ and mean climate conditions⁶⁵⁻⁶⁸. From a statistical perspective, our results suggest that the oxygenation of the SEP upper OMZ during strong EN events may have a residual on the long-term mean because cool events (La Niña episode) and CP EN events that are weaker in magnitude (as measured as the amplitude of the C index, see Methods) than EP EN events⁵⁴ have a weaker imprint on the OMZ. In other words, the ENSO asymmetry is transmitted to the SEP upper OMZ through a number

of processes (e.g., radiation of downwelling extra-tropical Rossby waves, downwelling Kelvin wave propagation, enhanced mesoscale activity, etc. see Pitcher et al.¹⁹). This is evidenced by the map of skewness of dissolved oxygen anomalies at 100 m for the model simulating a realistic ENSO asymmetry (see Methods, Fig. 4a). It shows an elongated pattern of positive skewness extending from the equatorial region (~10°S) to along the coast of Peru and Chile. A region of peak positive skewness is also found along the sections at 12°S and 30°S that matches with regions of influence of ENSO (rectangles in Fig. S2b, c) but with a more off-shore extension along the oxycline. Retaining all the models yields comparable results albeit with a positive skewness peaking more off-shore below 100 m (Fig. S6).

As El Niño events are expected to increase in frequency and amplitude⁵²,⁶⁹ as well as to last longer⁷⁰,⁷¹ in the warmer climate, we may thus expect a long-term increase in oxygen concentration in the SEP upper OMZ as greenhouse gases in the atmosphere increase. This does not preclude decadal variations as the increase in ENSO variance is not steady⁷⁰,⁷² and other ventilation processes can take place at long timescales¹⁸ (Cabré et al.¹⁶,²⁰,⁷³,⁷⁴). Consistently, we find that our ensemble model evidences a significant long-term mean increase in the variance of the E index (Fig. 1b; the slope of the linear fit of the 21-year running variance of the E index reaches 0.04 °C century⁻¹ which is significant at the 95% confidence level based on a
<figure><img src="image_7.png" /><figcaption>**Fig. 4** | Asymmetry of the SEP OMZ. Ensemble mean weighted skewness of dissolved oxygen anomalies (a) at 100 m, (b, c) along sections at 12°S and 30°S in the historical simulations (1920–2014). The contours in blue indicate the limit of the OMZ (i.e., 45 µmol/L isopleths) (at 300 m in a)). The linear trend has been removed</figcaption></figure>

prior to the calculation of weighted skewness (see Method). In this figure, only the models that simulate a realistic ENSO asymmetry, i.e., with a value of the α parameter lower than −0.15 were considered (see Methods for the definition of α and Table 1 of the supplementary material for the value of α).

student's t-test; see also Fig. S7a) associated with a mean increase in the frequency of occurrence of strong El Niño events from the present climate to the future climate (Fig. S7c). Considering the complex of processes involved in the SEP OMZ low-frequency variability, our findings may thus help interpreting proxy record of past SEP OMZ states⁷⁵⁻⁷⁹ noting that paleo OMZ studies mostly rely on mechanisms that involve low-frequency changes in the mean ocean circulation (i.e., changes in upwelling-driven biological export production and changes in ventilation associated with water formation in the ocean interior) but not in the high-frequency variability (i.e., ENSO). For instance, during the Miocene period, a period warmer than today's average temperatures, Hess et al.⁷⁸ showed that the eastern tropical Pacific was certainly well-oxygenated, suggesting a contracted SEP OMZ, which is attributed to a weaker Walker circulation. This period was also associated with enhanced ENSO variability⁸⁰, which could have also participated to oxygenate the SEP OMZ during that period. Paleoclimate reconstructions of the ENSO phenomenon covering the past ~10,000 years evidence an intensification of the ENSO variance in the late 20th century relative to other pre-industrial periods⁸¹,⁸². This could have also compensated the global deoxygenation trend due to the rapid rate of global warming since the mid 20th century in the SEP upper OMZ.

## Concluding remarks

The large dispersion in the simulation of the sensitivity of the OMZ to ENSO and tropical Pacific warming evidenced here calls for further investigating if its source is embedded in the physical or/and biological components of the ESMs⁸³, or if it results from errors propagation⁸⁴. Global ESMs have inherent limitations associated with their too low horizontal resolution to realistically account for coastal upwelling non-linear dynamics (Gruber et al.⁸⁵). Regional high-resolution coupled modeling and observational studies evidence in particular that mesoscale to submesoscale dynamics, not accounted for in the current generation of ESMs, are certainly important for explaining oxygen variability⁸⁶⁻⁹⁰. Global ESMs also currently simulate ENSO patterns and amplitude diversely⁹¹, which contributes to the spread of the oxygen sensitivity to ENSO-driven SST variability. The representation of the microbial loop driving the organic matter fate and the associated O₂ consumption, and thus part of the OMZs' dynamics⁹², is still also highly simplified in the current-generation biogeochemical models. More research is

thus needed to reduce the structural errors in the models, in order to enhance our confidence in the probabilistic interpretation of multi-model ensembles.

Given the low confidence in model projections, the existence of the emergent relationship evidenced here has thus practical value since it may help constrain the SEP upper OMZ sensitivity to climate change from observations. Currently, this is difficult to achieve observationally owing to the paucity of long-term oxygen data sets in this region of the world (Stramma et al.¹³,¹⁵). Off the coast of Peru (12°S), a monthly timeseries from 1996 is available²⁹ providing to some extent a benchmark for the sensitivity of DO to ENSO (see vertical red line in Fig. 3a; see Method in the caption text of Fig. 3a). However, limitations are associated with the record sampling (i.e., only one strong EP EN event (1997/1998)) and its location near the coast exposing it to natural variability associated with mesoscale processes and eutrophication effects. High-resolution long-term regional model simulations that may represent more realistically the OMZ than global ESMs can provide a clue on where the oxygen sensitivity to ENSO stands (see vertical line in blue color in Fig. 3 estimated from a regional biogeochemical coupled model simulation, see Methods), although confidence in their realism remains also dependent on the availability of observations. Previous studies have emphasized the sensitivity of the estimate of the observed DO trends along the coast of Peru and Chile to the sampling and period⁴¹,⁴⁴. Here we find that this is also the case for the DO ENSO sensitivity and that the estimate from Graco et al.'s²⁹ data may be too conservative (Fig. S8). Assuming that the regional model simulation provides a more realistic value of where the DO ENSO sensitivity stands, the DO climate sensitivity at 12°S and 30°S could reach 9.58 µmol L⁻¹ °C⁻¹ and 4.87 µmol L⁻¹ °C⁻¹, respectively (See Table S2 for error estimates), meaning that the upper margin of the OMZ would oxygenate twice faster off central Peru than off central Chile. More observations are needed to confirm this prediction. The current status of the observational record in the SEP in terms of the oxygen variable⁹³,⁹⁴ still makes it difficult to constrain the upper OMZ sensitivity to climate change using the emergent relationship established here. International programs and networks (TPOS2020, GO₂NE,
<figure><img src="image_8.png" /><figcaption>**Fig. 4** | Asymmetry of the SEP OMZ. Ensemble mean weighted skewness of dissolved oxygen anomalies (a) at 100 m, (b, c) along sections at 12°S and 30°S in the historical simulations (1920–2014). The contours in blue indicate the limit of the OMZ (i.e., 45 µmol/L isopleths) (at 300 m in a)). The linear trend has been removed</figcaption></figure>

prior to the calculation of weighted skewness (see Method). In this figure, only the models that simulate a realistic ENSO asymmetry, i.e., with a value of the α parameter lower than −0.15 were considered (see Methods for the definition of α and Table 1 of the supplementary material for the value of α).

student's t-test; see also Fig. S7a) associated with a mean increase in the frequency of occurrence of strong El Niño events from the present climate to the future climate (Fig. S7c). Considering the complex of processes involved in the SEP OMZ low-frequency variability, our findings may thus help interpreting proxy record of past SEP OMZ states⁷⁵⁻⁷⁹ noting that paleo OMZ studies mostly rely on mechanisms that involve low-frequency changes in the mean ocean circulation (i.e., changes in upwelling-driven biological export production and changes in ventilation associated with water formation in the ocean interior) but not in the high-frequency variability (i.e., ENSO). For instance, during the Miocene period, a period warmer than today's average temperatures, Hess et al.⁷⁸ showed that the eastern tropical Pacific was certainly well-oxygenated, suggesting a contracted SEP OMZ, which is attributed to a weaker Walker circulation. This period was also associated with enhanced ENSO variability⁸⁰, which could have also participated to oxygenate the SEP OMZ during that period. Paleoclimate reconstructions of the ENSO phenomenon covering the past ~10,000 years evidence an intensification of the ENSO variance in the late 20th century relative to other pre-industrial periods⁸¹,⁸². This could have also compensated the global deoxygenation trend due to the rapid rate of global warming since the mid 20th century in the SEP upper OMZ.

## Concluding remarks

The large dispersion in the simulation of the sensitivity of the OMZ to ENSO and tropical Pacific warming evidenced here calls for further investigating if its source is embedded in the physical or/and biological components of the ESMs⁸³, or if it results from errors propagation⁸⁴. Global ESMs have inherent limitations associated with their too low horizontal resolution to realistically account for coastal upwelling non-linear dynamics (Gruber et al.⁸⁵). Regional high-resolution coupled modeling and observational studies evidence in particular that mesoscale to submesoscale dynamics, not accounted for in the current generation of ESMs, are certainly important for explaining oxygen variability⁸⁶⁻⁹⁰. Global ESMs also currently simulate ENSO patterns and amplitude diversely⁹¹, which contributes to the spread of the oxygen sensitivity to ENSO-driven SST variability. The representation of the microbial loop driving the organic matter fate and the associated O₂ consumption, and thus part of the OMZs' dynamics⁹², is still also highly simplified in the current-generation biogeochemical models. More research is

thus needed to reduce the structural errors in the models, in order to enhance our confidence in the probabilistic interpretation of multi-model ensembles.

Given the low confidence in model projections, the existence of the emergent relationship evidenced here has thus practical value since it may help constrain the SEP upper OMZ sensitivity to climate change from observations. Currently, this is difficult to achieve observationally owing to the paucity of long-term oxygen data sets in this region of the world (Stramma et al.¹³,¹⁵). Off the coast of Peru (12°S), a monthly timeseries from 1996 is available²⁹ providing to some extent a benchmark for the sensitivity of DO to ENSO (see vertical red line in Fig. 3a; see Method in the caption text of Fig. 3a). However, limitations are associated with the record sampling (i.e., only one strong EP EN event (1997/1998)) and its location near the coast exposing it to natural variability associated with mesoscale processes and eutrophication effects. High-resolution long-term regional model simulations that may represent more realistically the OMZ than global ESMs can provide a clue on where the oxygen sensitivity to ENSO stands (see vertical line in blue color in Fig. 3 estimated from a regional biogeochemical coupled model simulation, see Methods), although confidence in their realism remains also dependent on the availability of observations. Previous studies have emphasized the sensitivity of the estimate of the observed DO trends along the coast of Peru and Chile to the sampling and period⁴¹,⁴⁴. Here we find that this is also the case for the DO ENSO sensitivity and that the estimate from Graco et al.'s²⁹ data may be too conservative (Fig. S8). Assuming that the regional model simulation provides a more realistic value of where the DO ENSO sensitivity stands, the DO climate sensitivity at 12°S and 30°S could reach 9.58 µmol L⁻¹ °C⁻¹ and 4.87 µmol L⁻¹ °C⁻¹, respectively (See Table S2 for error estimates), meaning that the upper margin of the OMZ would oxygenate twice faster off central Peru than off central Chile. More observations are needed to confirm this prediction. The current status of the observational record in the SEP in terms of the oxygen variable⁹³,⁹⁴ still makes it difficult to constrain the upper OMZ sensitivity to climate change using the emergent relationship established here. International programs and networks (TPOS2020, GO₂NE,
ARGO) have recently made recommendations and provided pathways for filling in that gap¹⁵,⁹⁵⁻⁹⁷. Our results thus offer an incentive for pursuing observational efforts deployed in the framework of regional programs (cf. the SEPICAF (South Eastern Pacific Circulation from Argo floats) project, a contribution to the International BGC-Argo program for the SEP that emerged from the recommendations of the TPOS2020 report). More importantly, it also provides the foundation of an alternative storyline for the evolution of marine ecosystems in this region (i.e., ocean oxygenation), which may guide the development of adaptation strategies. In particular, while much attention has been drawn to ocean deoxygenation at a global scale and its negative effects on marine biota, the fate of the SEP upper OMZ is more likely to be dominated by oxygenation.

# Methods

## Observations

Climatological oxygen concentrations are from the CARS (CSIRO Atlas of Regional Seas) 2009 data set (Ridgway et al.⁹⁸). SST data (1870–2019) are from HadISST v1.1 (Hadley Centre Sea Ice and Sea Surface Temperature version 1.1) (Rayner et al.⁹⁹).

## Earth system models

Oxygen model outputs of 19 models belonging to the Coupled Model Intercomparison Project phase5 (CMIP5¹⁰⁰) and 6 (CMIP6¹⁰¹) were used (See lists in Table S1). Additionally, we used 34 members of the CESM Large Ensemble¹⁰² to evaluate the sensitivity of the result to natural variability and to temporal resolution (monthly mean versus yearly mean) of the data since most models provide only yearly mean data for oxygen (see Supplementary Fig. S9). The historical (1950–2014) scenario is used as a surrogate of the present climate, whereas the Representative Concentration Pathway (RCP) 8.5 (RCP8.5) and the Socioeconomic Pathway (SSP) 8.5 (SSP85), which corresponds to a high greenhouse gas emissions pathway for which a radiative forcing of 8.5 W m⁻² is reached at the end of the century¹⁰³, is used as future climate scenario (2015–2100). Model data were bilinearly interpolated horizontally on a similar grid with resolution of 1.5° × 1.5° and linearly interpolated on the vertical using 60 vertical levels (native vertical grid of the CESM1-BCG model).

## Regional model simulation

Oxygen data of a long-term regional biogeochemical coupled model simulation are also used as a benchmark for evaluating the global ESMs in terms of their sensitivity to ENSO (Fig. 3a, b). The simulation covers the period 1958–2008. The model is based on the biogeochemical model of Gutknecht et al.¹⁰⁴,¹⁰⁵ tuned for realistically simulating the Peruvian OMZ²⁶,⁸⁷ coupled to the hydro-dynamical solution of Dewitte et al.³³. The simulation was extensively validated from both satellite and in situ observations, and was shown to account realistically for the changes in oceanic circulation off Peru and Chile during the 3 intense El Niño that it samples (1972/73, 1982/ 83, 1997/1998)³⁵,⁴³. It has been used to study the OMZ dynamics off Peru⁸⁶ and off Chile⁸⁸,¹⁰⁶.

## ENSO indices

We use objectively defined indices for Eastern Pacific (EP) El Niño and Central Pacific (CP) El Niño occurrence and amplitude in order to reduce uncertainties associated with the models simulating a large diversity in ENSO patterns⁵². These indices, named the E and C indices, are defined as in Takahashi et al.⁵⁴, i.e., based on the first two Empirical Orthogonal Functions (EOFs) of SST anomalies over the region (10°N–10°S; 120°E–90°W). $E = \frac{PC1-PC2}{\sqrt{2}}$ and $C = \frac{PC1+PC2}{\sqrt{2}}$ where PC1 and PC2 are the first two principal components of the EOF modes that have been normalized by their corresponding variance. The E and C indices have no unit unless they have been dimensionalized (in °C) by the norm of their corresponding patterns ($E_{xy}$ and $C_{xy}$). $E_{xy}$ and $C_{xy}$ are obtained by bilinearly regressing the SST anomalies onto the E and C indices. We define their norm as their root mean

square, i.e.: $||E_{xy}|| = \sqrt{\frac{1}{L_x L_y} \iint E_{xy}^2(x, y) \cdot dxdy}$ (resp. $C_{xy}$)), where $L_x$ and $L_y$ are the zonal and meridional extents of the domain used to calculate the mode patterns (i.e., (10°N–10°S; 120°E–90°W)). The E index accounts for the amplitude of SST interannual anomalies during EP El Niño events which center of action is strongest in the far eastern equatorial Pacific and along the coast of Peru, while the C index, accounts for Central Pacific El Niño (CP EN) and La Niña events (Supplementary Fig. S10). These two indices are independent (correlation between them is zero) so that they can be conveniently used for inferring the contribution of ENSO variability to a particular variable through bilinear regression. Yearly mean data are considered for consistency with the oxygen data analyses based on yearly mean data. The data are linearly detrended over the future period (2015–2100). For the calculation of the regression coefficient between DO and the ENSO indices (so-called oxygen ENSO sensitivity), the E and C indices are dimensionized in order to have unit in µmol L⁻¹ °C⁻¹. Since the E and C indices are derived separately for the two periods (present: 1920–2014 and future: 2015–2100), the indices of the future need also to be scaled so as to take into account the changes in amplitude and “shape” of their corresponding patterns from the historical period. This consists in projecting the mode patterns of the future period onto the mode patterns of the present period. This scaling factor is estimated as follows: $\frac{1}{L_x L_y} \iint E_{xy}^{\text{future}}(x, y) \cdot E_{xy}^{\text{present}}(x, y) \cdot dxdy$ where $L_x$ and $L_y$ are the zonal and meridional extent of the domain used to calculate the mode patterns (i.e., (10°N–10°S; 120°E–90°W)) and the overline indicates that the mode pattern was normalized by its corresponding root mean square, i.e., $\bar{E}_{xy} = E_{xy} / \sqrt{\frac{1}{L_x L_y} \iint E_{xy}^2(x, y) \cdot dxdy}$. This procedure follows Carréric et al.⁷⁰. Note that this method yields comparable results if the E and C indices (resp. $E_{xy}$, $C_{xy}$) are estimated based on the EOF decomposition of monthly mean SST anomalies performed over the whole period (1920–2100) as in Cai et al.⁵². The Figs. 1b and S7 that are based on monthly mean SST data are derived using the Cai et al.’s⁵² approach. A difference with Cai et al.⁵² is the removal of a linear trend for the period 2015–2100 instead of a quadratic trend over the entire record, which was checked to have little impact on the results. Monthly mean anomalies are relative to the climatology calculated over the period 1920–2020.

The ENSO asymmetry is diagnosed as in Karamperidou et al.¹⁰⁷, that is from the $\alpha$ parameter that corresponds to the first coefficient of the quadratic fit of PC2 onto PC1, i.e., PC2\_approx(t) = $\alpha \cdot PC1^2(t) + \beta \cdot PC1(t) + \gamma$ where PC1 and PC2 are the timeseries associated with the first two EOF mode of SST anomalies in the tropical Pacific. The value of $\alpha$ is provided in Table S1 for the models and the observations. The smaller value of $\alpha$, the stronger the ENSO asymmetry⁵².

## Oxygen climate sensitivity

In order to estimate the sensitivity of the upper bound of the SEP OMZ to tropical warming, we calculate the ratio of the long-term trend (2014–2100) in dissolved oxygen concentration in the vicinity of the oxycline (average between 25 m and 200 m, and from the coast to 3° off-shore, see rectangles in Fig. S2b, c) and the long-term trend (2014–2100) of SST projected onto the E mode pattern ($E_{xy}$). The latter is motivated by the fact that the E mode pattern matches the region of peak long-term SST warming, which has a so-called El Niño-like pattern (i.e., stronger SST warming in the tropical eastern Pacific relative to the western Pacific; see Fig. S11). The SST in the E region for each model is then calculated as follows: $\frac{1}{L_x L_y} \iint SST(x, y) \cdot \bar{E}_{xy}(x, y) \cdot dxdy$ where $L_x$ and $L_y$ are the zonal and meridional extents of the domain used to calculate the E mode pattern (10°N–10°S; 120°E–90°W) and $\bar{E}_{xy}$ is the normalized E mode pattern (see above).

## Oxygen skewness

The skewness is a normalized third statistical moment¹⁰⁸. Since a small standard deviation may cause large skewness rather than the normalized
ARGO) have recently made recommendations and provided pathways for filling in that gap¹⁵,⁹⁵⁻⁹⁷. Our results thus offer an incentive for pursuing observational efforts deployed in the framework of regional programs (cf. the SEPICAF (South Eastern Pacific Circulation from Argo floats) project, a contribution to the International BGC-Argo program for the SEP that emerged from the recommendations of the TPOS2020 report). More importantly, it also provides the foundation of an alternative storyline for the evolution of marine ecosystems in this region (i.e., ocean oxygenation), which may guide the development of adaptation strategies. In particular, while much attention has been drawn to ocean deoxygenation at a global scale and its negative effects on marine biota, the fate of the SEP upper OMZ is more likely to be dominated by oxygenation.

# Methods

## Observations

Climatological oxygen concentrations are from the CARS (CSIRO Atlas of Regional Seas) 2009 data set (Ridgway et al.⁹⁸). SST data (1870–2019) are from HadISST v1.1 (Hadley Centre Sea Ice and Sea Surface Temperature version 1.1) (Rayner et al.⁹⁹).

## Earth system models

Oxygen model outputs of 19 models belonging to the Coupled Model Intercomparison Project phase5 (CMIP5¹⁰⁰) and 6 (CMIP6¹⁰¹) were used (See lists in Table S1). Additionally, we used 34 members of the CESM Large Ensemble¹⁰² to evaluate the sensitivity of the result to natural variability and to temporal resolution (monthly mean versus yearly mean) of the data since most models provide only yearly mean data for oxygen (see Supplementary Fig. S9). The historical (1950–2014) scenario is used as a surrogate of the present climate, whereas the Representative Concentration Pathway (RCP) 8.5 (RCP8.5) and the Socioeconomic Pathway (SSP) 8.5 (SSP85), which corresponds to a high greenhouse gas emissions pathway for which a radiative forcing of 8.5 W m⁻² is reached at the end of the century¹⁰³, is used as future climate scenario (2015–2100). Model data were bilinearly interpolated horizontally on a similar grid with resolution of 1.5° × 1.5° and linearly interpolated on the vertical using 60 vertical levels (native vertical grid of the CESM1-BCG model).

## Regional model simulation

Oxygen data of a long-term regional biogeochemical coupled model simulation are also used as a benchmark for evaluating the global ESMs in terms of their sensitivity to ENSO (Fig. 3a, b). The simulation covers the period 1958–2008. The model is based on the biogeochemical model of Gutknecht et al.¹⁰⁴,¹⁰⁵ tuned for realistically simulating the Peruvian OMZ²⁶,⁸⁷ coupled to the hydro-dynamical solution of Dewitte et al.³³. The simulation was extensively validated from both satellite and in situ observations, and was shown to account realistically for the changes in oceanic circulation off Peru and Chile during the 3 intense El Niño that it samples (1972/73, 1982/ 83, 1997/1998)³⁵,⁴³. It has been used to study the OMZ dynamics off Peru⁸⁶ and off Chile⁸⁸,¹⁰⁶.

## ENSO indices

We use objectively defined indices for Eastern Pacific (EP) El Niño and Central Pacific (CP) El Niño occurrence and amplitude in order to reduce uncertainties associated with the models simulating a large diversity in ENSO patterns⁵². These indices, named the E and C indices, are defined as in Takahashi et al.⁵⁴, i.e., based on the first two Empirical Orthogonal Functions (EOFs) of SST anomalies over the region (10°N–10°S; 120°E–90°W). $E = \frac{PC1-PC2}{\sqrt{2}}$ and $C = \frac{PC1+PC2}{\sqrt{2}}$ where PC1 and PC2 are the first two principal components of the EOF modes that have been normalized by their corresponding variance. The E and C indices have no unit unless they have been dimensionalized (in °C) by the norm of their corresponding patterns ($E_{xy}$ and $C_{xy}$). $E_{xy}$ and $C_{xy}$ are obtained by bilinearly regressing the SST anomalies onto the E and C indices. We define their norm as their root mean

square, i.e.: $||E_{xy}|| = \sqrt{\frac{1}{L_x L_y} \iint E_{xy}^2(x, y) .dxdy}$ (resp. $C_{xy}$)), where $L_x$ and $L_y$ are the zonal and meridional extents of the domain used to calculate the mode patterns (i.e., (10°N–10°S; 120°E–90°W)). The E index accounts for the amplitude of SST interannual anomalies during EP El Niño events which center of action is strongest in the far eastern equatorial Pacific and along the coast of Peru, while the C index, accounts for Central Pacific El Niño (CP EN) and La Niña events (Supplementary Fig. S10). These two indices are independent (correlation between them is zero) so that they can be conveniently used for inferring the contribution of ENSO variability to a particular variable through bilinear regression. Yearly mean data are considered for consistency with the oxygen data analyses based on yearly mean data. The data are linearly detrended over the future period (2015–2100). For the calculation of the regression coefficient between DO and the ENSO indices (so-called oxygen ENSO sensitivity), the E and C indices are dimensionized in order to have unit in µmol L⁻¹ °C⁻¹. Since the E and C indices are derived separately for the two periods (present: 1920–2014 and future: 2015–2100), the indices of the future need also to be scaled so as to take into account the changes in amplitude and “shape” of their corresponding patterns from the historical period. This consists in projecting the mode patterns of the future period onto the mode patterns of the present period. This scaling factor is estimated as follows: $\frac{1}{L_x L_y} \iint \frac{E_{xy}^{\text{future}}(x, y) \cdot E_{xy}^{\text{present}}(x, y)}{\sqrt{\frac{1}{L_x L_y} \iint E_{xy}^2(x, y) \cdot dxdy}}$. This procedure follows Carréric et al.⁷⁰. Note that this method yields comparable results if the E and C indices (resp. $E_{xy}$, $C_{xy}$) are estimated based on the EOF decomposition of monthly mean SST anomalies performed over the whole period (1920–2100) as in Cai et al.⁵². The Figs. 1b and S7 that are based on monthly mean SST data are derived using the Cai et al.’s⁵² approach. A difference with Cai et al.⁵² is the removal of a linear trend for the period 2015–2100 instead of a quadratic trend over the entire record, which was checked to have little impact on the results. Monthly mean anomalies are relative to the climatology calculated over the period 1920–2020.

The ENSO asymmetry is diagnosed as in Karamperidou et al.¹⁰⁷, that is from the $\alpha$ parameter that corresponds to the first coefficient of the quadratic fit of PC2 onto PC1, i.e., PC2\_approx(t) = $\alpha$.PC1²(t) + $\beta$PC1(t) + $\gamma$ where PC1 and PC2 are the timeseries associated with the first two EOF mode of SST anomalies in the tropical Pacific. The value of $\alpha$ is provided in Table S1 for the models and the observations. The smaller value of $\alpha$, the stronger the ENSO asymmetry⁵².

## Oxygen climate sensitivity

In order to estimate the sensitivity of the upper bound of the SEP OMZ to tropical warming, we calculate the ratio of the long-term trend (2014–2100) in dissolved oxygen concentration in the vicinity of the oxycline (average between 25 m and 200 m, and from the coast to 3° off-shore, see rectangles in Fig. S2b, c) and the long-term trend (2014–2100) of SST projected onto the E mode pattern ($E_{xy}$). The latter is motivated by the fact that the E mode pattern matches the region of peak long-term SST warming, which has a so-called El Niño-like pattern (i.e., stronger SST warming in the tropical eastern Pacific relative to the western Pacific; see Fig. S11). The SST in the E region for each model is then calculated as follows: $\frac{1}{L_x L_y} \iint SST(x, y) \cdot \bar{E}_{xy}(x, y) \cdot dxdy$ where $L_x$ and $L_y$ are the zonal and meridional extents of the domain used to calculate the E mode pattern (10°N–10°S; 120°E–90°W) and $\bar{E}_{xy}$ is the normalized E mode pattern (see above).

## Oxygen skewness

The skewness is a normalized third statistical moment¹⁰⁸. Since a small standard deviation may cause large skewness rather than the normalized
skewness, we used the weighted skewness, m3/m2 where

$$m_k = \sum_{i=1}^{N} (x_i - \bar{x})^k / N$$

$x_i$ is the ith data point, N is the number of data (years) and $\bar{x}$ is the mean.

### Statistical significance test

Regression coefficient's significances are estimated based on a Student's t test. For RCP8.5 (SSP5-8.5) simulations, we calculated linear trends over the period 2006–2100 (2015–2100). In order to estimate the significance of the ensemble mean long-term linear trend values, we set the requirement to have at least 80% of the models to yield the same sign of the trend values¹⁰⁹,¹¹⁰.

### Data availability

The CMIP5 models is available at: https://aims2.lnl.gov/search/cmip5/ The CMIP6 models is available at: https://aims2.lnl.gov/search/cmip6/ SST data from HadISST is available at: https://www.metoffice.gov.uk/hadobs/ hadisst/data/download.html Dissolved Oxygen data from CARS2009 is available at: https://catalogue-imos.aodn.org.au/geonetwork/srv/eng/ catalog.search#/metadata/d9302a48-57b1-41c2-a0dc-78bd00dd5e4b.

Received: 28 August 2023; Accepted: 29 April 2024;
Published online: 28 May 2024

### References

1. Gallo, N. D. & Levin, L. A. Fish ecology and evolution in the world's oxygen minimum zones and implications of ocean deoxygenation. *Adv. Mar. Biol.* **74**, 117–198 (2016).

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

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

4. Wyrtki, K. Physical oceanography of the Southeast Asian waters (Vol. 2). University of California, Scripps Institution of Oceanography. (1961).

5. Karstensen, J., Stramma, L. & Visbeck, M. Oxygen minimum zones in the eastern tropical Atlantic and Pacific oceans. *Prog. Oceanogr.* **77**, 331–350 (2008).

6. Helly, J. J. & Levin, L. A. Global distribution of naturally occurring marine hypoxia on continental margins. *Deep Sea Res. Part I: Oceanogr. Res. Pap.* **51**, 1159–1168 (2004).

7. Paulmier, A. & Ruiz-Pino, D. Oxygen minimum zones (OMZs) in the modern ocean. *Prog. Oceanogr.* **80**, 113–128 (2009).

8. Kalvelage, T. et al. Aerobic microbial respiration in oceanic oxygen minimum zones. *PloS One* **10**, e0133526 (2015).

9. Lam, P. et al. Revising the nitrogen cycle in the Peruvian oxygen minimum zone. *Proc. Natl Acad. Sci.* **106**, 4752–4757 (2009).

10. Wishner, K. F., Seibel, B. & Outram, D. Ocean deoxygenation and copepods: coping with oxygen minimum zone variability. *Biogeosciences* **17**, 2315–2339 (2020).

11. Bertrand, A. et al. Oxygen: a fundamental property regulating pelagic ecosystem structure in the coastal southeastern tropical Pacific. *PloS One* **6**, e29558 (2011).

12. Paulmier, A. et al. High-sustained concentrations of organisms at very low oxygen concentration indicated by acoustic profiles in the oxygen deficit region off Peru. *Front. Mar. Sci.* **8**, 723056 (2021).

13. Stramma, L., Johnson, G. C., Sprintall, J. & Mohrholz, V. Expanding oxygen-minimum zones in the tropical oceans. *Science* **320**, 655–658 (2008).

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

15. Grégoire, M. et al. A global ocean oxygen database and Atlas for assessing and predicting deoxygenation and ocean health in the open and coastal ocean. *Front. Mar. Sci.* **1638** (2021).

16. Cabré, A., Marinov, I., Bernardello, R. & Bianchi, D. Oxygen minimum zones in the tropical Pacific across CMIP5 models: mean state differences and climate change trends. *Biogeosciences* **12**, 5429–5454 (2015).

17. Busecke, J. J. M., Resplandy, L., Ditkovsky, S. J. & John, J. G. Diverging fates of the Pacific Ocean oxygen minimum zone and its core in a warming world. *AGU Adv.* **3**, e2021AV000470 (2022).

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

19. Pitcher G.C. et al. System controls of coastal and open ocean oxygen depletion. *Prog. Oceanogr.* https://doi.org/10.1016/j. pocean.2021.102613. (2021)

20. Duteil, O., Frenger, I. & Getzlaff, J. The riddle of eastern tropical Pacific ocean oxygen levels: the role of the supply by intermediate-depth waters. *Ocean Sci.* **17**, 1489–1507 (2021).

21. Zhu, C., Liu, Z. & Gu, S. Model bias for South Atlantic Antarctic intermediate water in CMIP5. *Clim. Dyn.* **50**, 3613–3624 (2018).

22. Sloyan, B. M. & Kamenkovich, I. V. Simulation of subantarctic mode and Antarctic intermediate waters in climate models. *J. Clim.* **20**, 5061–5080 (2007).

23. Fuenzalida, R., Schneider, W., Garcés-Vargas, J., Bravo, L. & Lange, C. Vertical and horizontal extension of the oxygen minimum zone in the eastern South Pacific Ocean. *Deep Sea Res. Part II: Top. Stud. Oceanogr.* **56**, 992–1003 (2009).

24. Johnson, G. C. & Moore, D. W. The Pacific subsurface countercurrents and an inertial model. *J. Phys. Oceanogr.* **27**, 2448–2459 (1997).

25. Tsuchiya, M. A subsurface north equatorial countercurrent in the eastern Pacific Ocean. *J. Geophys. Res.* **77**, 5981–5986 (1972).

26. Montes, I. et al. High-resolution modeling of the Eastern Tropical Pacific oxygen minimum zone: sensitivity to the tropical oceanic circulation. *J. Geophys. Res.: Oceans* **119**, 5515–5532 (2014).

27. Hall, A. & Qu, X. Using the current seasonal cycle to constrain snow albedo feedback in future climate change. *Geophys. Res. Lett.* **33**, L03502 (2006).

28. De La Maza, L. & Farías, L. The intensification of coastal hypoxia off central Chile: long term and high frequency variability. *Front. Earth Sci.* **10**, 929271 (2023).

29. Graco, M. I. et al. The OMZ and nutrient features as a signature of interannual and low-frequency variability in the Peruvian upwelling system. *Biogeosciences* **14**, 4601–4617 (2017).

30. Morales, C., Hormazabal, S. & Blanco, J. L. Interannual variability in the mesoscale distribution of the depth of the upper boundary of the oxygen minimum layer off northern Chile (18°–24°S): Implications for the pelagic system and biogeochemical cycling. *J. Mar. Res.* **57**, 909–932 (1999).

31. Vergara, O., B. Dewitte, M. Ramos and O. Pizarro. Vertical energy flux at ENSO time scales in the subthermocline of the Southeastern Pacific. *J. Geophys. Res. Oceans* **122**, https://doi.org/10.1002/ 2016JC012614. (2017).

32. Clarke, A. J. & Van Gorder, S. On ENSO coastal currents and sea levels. *J. Phys. Oceanogr.* **24**, 661–680 (1994).

33. Dewitte, B. et al. Change in El Niño flavours over 1958–2008: implications for the long-term trend of the upwelling off Peru. *Deep Sea Res. Part II: Top. Stud. Oceanogr.* **77**, 143–156 (2012).

34. Pizarro, O., Clarke, A. J. & Van Gorder, S. El Niño sea level and currents along the South American coast: comparison of observations with theory. *J. Phys. Oceanogr.* **31**, 1891–1903 (2001).

35. Pizarro, O., Shaffer, G., Dewitte, B. & Ramos, M. Dynamics of seasonal and interannual variability of the Peru-Chile undercurrent. *Geophys. Res. Lett.* **29**, 22–1 (2002).
skewness, we used the weighted skewness, m3/m2 where

$$m_k = \sum_{i=1}^{N} (x_i - \bar{x})^k / N$$

$x_i$ is the ith data point, N is the number of data (years) and $\bar{x}$ is the mean.

## Statistical significance test

Regression coefficient's significances are estimated based on a Student's t test. For RCP8.5 (SSP5-8.5) simulations, we calculated linear trends over the period 2006–2100 (2015–2100). In order to estimate the significance of the ensemble mean long-term linear trend values, we set the requirement to have at least 80% of the models to yield the same sign of the trend values¹⁰⁹,¹¹⁰.

## Data availability

The CMIP5 models is available at: https://aims2.llnl.gov/search/cmip5/ The CMIP6 models is available at: https://aims2.llnl.gov/search/cmip6/ SST data from HadISST is available at: https://www.metoffice.gov.uk/hadobs/ hadisst/data/download.html Dissolved Oxygen data from CARS2009 is available at: https://catalogue-imos.aodn.org.au/geonetwork/srv/eng/ catalog.search#/metadata/d9302a48-57b1-41c2-a0dc-78bd00dd5e4b.

Received: 28 August 2023; Accepted: 29 April 2024;
Published online: 28 May 2024

## References

1. Gallo, N. D. & Levin, L. A. Fish ecology and evolution in the world's oxygen minimum zones and implications of ocean deoxygenation. *Adv. Mar. Biol.* **74**, 117–198 (2016).

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

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

4. Wyrtki, K. Physical oceanography of the Southeast Asian waters (Vol. 2). University of California, Scripps Institution of Oceanography. (1961).

5. Karstensen, J., Stramma, L. & Visbeck, M. Oxygen minimum zones in the eastern tropical Atlantic and Pacific oceans. *Prog. Oceanogr.* **77**, 331–350 (2008).

6. Helly, J. J. & Levin, L. A. Global distribution of naturally occurring marine hypoxia on continental margins. *Deep Sea Res. Part I: Oceanogr. Res. Pap.* **51**, 1159–1168 (2004).

7. Paulmier, A. & Ruiz-Pino, D. Oxygen minimum zones (OMZs) in the modern ocean. *Prog. Oceanogr.* **80**, 113–128 (2009).

8. Kalvelage, T. et al. Aerobic microbial respiration in oceanic oxygen minimum zones. *PloS One* **10**, e0133526 (2015).

9. Lam, P. et al. Revising the nitrogen cycle in the Peruvian oxygen minimum zone. *Proc. Natl Acad. Sci.* **106**, 4752–4757 (2009).

10. Wishner, K. F., Seibel, B. & Outram, D. Ocean deoxygenation and copepods: coping with oxygen minimum zone variability. *Biogeosciences* **17**, 2315–2339 (2020).

11. Bertrand, A. et al. Oxygen: a fundamental property regulating pelagic ecosystem structure in the coastal southeastern tropical Pacific. *PloS One* **6**, e29558 (2011).

12. Paulmier, A. et al. High-sustained concentrations of organisms at very low oxygen concentration indicated by acoustic profiles in the oxygen deficit region off Peru. *Front. Mar. Sci.* **8**, 723056 (2021).

13. Stramma, L., Johnson, G. C., Sprintall, J. & Mohrholz, V. Expanding oxygen-minimum zones in the tropical oceans. *Science* **320**, 655–658 (2008).

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

15. Grégoire, M. et al. A global ocean oxygen database and Atlas for assessing and predicting deoxygenation and ocean health in the open and coastal ocean. *Front. Mar. Sci.* **1638** (2021).

16. Cabré, A., Marinov, I., Bernardello, R. & Bianchi, D. Oxygen minimum zones in the tropical Pacific across CMIP5 models: mean state differences and climate change trends. *Biogeosciences* **12**, 5429–5454 (2015).

17. Busecke, J. J. M., Resplandy, L., Ditkovsky, S. J. & John, J. G. Diverging fates of the Pacific Ocean oxygen minimum zone and its core in a warming world. *AGU Adv.* **3**, e2021AV000470 (2022).

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

19. Pitcher G.C. et al. System controls of coastal and open ocean oxygen depletion. *Prog. Oceanogr.* https://doi.org/10.1016/j. pocean.2021.102613. (2021)

20. Duteil, O., Frenger, I. & Getzlaff, J. The riddle of eastern tropical Pacific ocean oxygen levels: the role of the supply by intermediate-depth waters. *Ocean Sci.* **17**, 1489–1507 (2021).

21. Zhu, C., Liu, Z. & Gu, S. Model bias for South Atlantic Antarctic intermediate water in CMIP5. *Clim. Dyn.* **50**, 3613–3624 (2018).

22. Sloyan, B. M. & Kamenkovich, I. V. Simulation of subantarctic mode and Antarctic intermediate waters in climate models. *J. Clim.* **20**, 5061–5080 (2007).

23. Fuenzalida, R., Schneider, W., Garcés-Vargas, J., Bravo, L. & Lange, C. Vertical and horizontal extension of the oxygen minimum zone in the eastern South Pacific Ocean. *Deep Sea Res. Part II: Top. Stud. Oceanogr.* **56**, 992–1003 (2009).

24. Johnson, G. C. & Moore, D. W. The Pacific subsurface countercurrents and an inertial model. *J. Phys. Oceanogr.* **27**, 2448–2459 (1997).

25. Tsuchiya, M. A subsurface north equatorial countercurrent in the eastern Pacific Ocean. *J. Geophys. Res.* **77**, 5981–5986 (1972).

26. Montes, I. et al. High-resolution modeling of the Eastern Tropical Pacific oxygen minimum zone: sensitivity to the tropical oceanic circulation. *J. Geophys. Res.: Oceans* **119**, 5515–5532 (2014).

27. Hall, A. & Qu, X. Using the current seasonal cycle to constrain snow albedo feedback in future climate change. *Geophys. Res. Lett.* **33**, L03502 (2006).

28. De La Maza, L. & Farías, L. The intensification of coastal hypoxia off central Chile: long term and high frequency variability. *Front. Earth Sci.* **10**, 929271 (2023).

29. Graco, M. I. et al. The OMZ and nutrient features as a signature of interannual and low-frequency variability in the Peruvian upwelling system. *Biogeosciences* **14**, 4601–4617 (2017).

30. Morales, C., Hormazabal, S. & Blanco, J. L. Interannual variability in the mesoscale distribution of the depth of the upper boundary of the oxygen minimum layer off northern Chile (18°–24°S): Implications for the pelagic system and biogeochemical cycling. *J. Mar. Res.* **57**, 909–932 (1999).

31. Vergara, O., B. Dewitte, M. Ramos and O. Pizarro. Vertical energy flux at ENSO time scales in the subthermocline of the Southeastern Pacific. *J. Geophys. Res. Oceans* **122**, https://doi.org/10.1002/ 2016JC012614. (2017).

32. Clarke, A. J. & Van Gorder, S. On ENSO coastal currents and sea levels. *J. Phys. Oceanogr.* **24**, 661–680 (1994).

33. Dewitte, B. et al. Change in El Niño flavours over 1958–2008: implications for the long-term trend of the upwelling off Peru. *Deep Sea Res. Part II: Top. Stud. Oceanogr.* **77**, 143–156 (2012).

34. Pizarro, O., Clarke, A. J. & Van Gorder, S. El Niño sea level and currents along the South American coast: comparison of observations with theory. *J. Phys. Oceanogr.* **31**, 1891–1903 (2001).

35. Pizarro, O., Shaffer, G., Dewitte, B. & Ramos, M. Dynamics of seasonal and interannual variability of the Peru-Chile undercurrent. *Geophys. Res. Lett.* **29**, 22–1 (2002).
36. Sprintall, J., S. Cravatte, B. Dewitte, Y. Du and A. S. Gupta. Oceanic teleconnections, chapter 15 in "El Niño in a changing climate" AGU Book, ISBN: 978-1-119-54816-4, 528 (2020).

37. Köhn, E. E., Münnich, M., Vogt, M., Desmet, F. & Gruber, N. Strong habitat compression by extreme shoaling events of hypoxic waters in the Eastern Pacific. *J. Geophys. Res.: Oceans* **127**, e2022JC018429 (2022).

38. Leung, S., Thompson, L., McPhaden, M. J. & Mislan, K. A. S. ENSO drives near-surface oxygen and vertical habitat variability in the tropical Pacific. *Environ. Res. Lett.* **14**, 060402 (2019).

39. Busecke, J. J. M., Resplandy, L. & Dunne, J. P. The equatorial undercurrent and the oxygen minimum zone in the Pacific. *Geophys. Res. Lett.* **46**, 6716–6725 (2019).

40. Shigemitsu, M., Yamamoto, A., Oka, A. & Yamanaka, Y. One possible uncertainty in CMIP5 projections of low-oxygen water volume in the Eastern Tropical Pacific. *Glob. Biogeochem. Cycles* **31**, 804–820 (2017).

41. Ito, T. & Deutsch, C. Variability of the Oxygen Minimum Zone in the Tropical North Pacific during the late twentieth century. *Glob. Biogeochem. Cycles* **27**, 1119–1128 (2013).

42. José, Y. S., Stramma, L., Schmidtko, S., & Oschlies, A. (2019). ENSO-driven fluctuations in oxygen supply and vertical extent of oxygen-poor waters in the oxygen minimum zone of the Eastern Tropical South Pacific. *Biogeosci. Discuss.* 1–20 (2019).

43. Conejero, C., Dewitte, B., Garçon, V., Sudre, J. & Montes, I. ENSO diversity driving low-frequency change in mesoscale activity off Peru and Chile. *Sci. Rep.* **10**, 17902 (2020).

44. Espinoza-Morriberón, D. et al. Oxygen variability during ENSO in the tropical South Eastern Pacific. *Front. Mar. Sci.* **5**, 526 (2019).

45. Gutierrez, D. et al. Oxygenation episodes on the continental shelf of central Peru: remote forcing and benthic ecosystem response. *Prog. Oceanogr.* **79**, 177–189 (2008).

46. Mogollón, R. & Calil, P. H. R. On the effects of ENSO on ocean biogeochemistry in the Northern Humboldt Current System (NHCS): a modeling study. *J. Mar. Syst.* **172**, 137–159 (2017).

47. Duteil, O., Oschlies, A. & Böning, C. W. Pacific decadal oscillation and recent oxygen decline in the eastern tropical Pacific Ocean. *Biogeosciences (BG)* **15**, 7111–7126 (2018).

48. Pizarro, O., & Montecinos, A. Interdecadal variability of the thermocline along the west coast of South America. *Geophys. Res. Lett.* **31** (2004).

49. Oschlies, A. et al. Patterns of deoxygenation: sensitivity to natural and anthropogenic drivers. *Philos. Trans. R. Soc. A* **375**, 20160325 (2017).

50. Ito, T. & Deutsch, C. A conceptual model for the temporal spectrum of oceanic oxygen variability. *Geophys. Res. Lett.* **37**, L03601 (2010).

51. Capotondi, A., Wittenberg, A. T., Kug, J. S., Takahashi, K., & McPhaden, M. J. ENSO diversity. El Niño Southern Oscillation in a changing climate. **4**, 65–86 (2020).

52. Cai, W. et al. Increased variability of eastern Pacific El Niño under greenhouse warming. *Nature* **564**, 201–206 (2018).

53. Cai, W. et al. Changing El Niño–Southern Oscillation in a warming climate. *Nat. Rev. Earth Environ.* **2**, 628–644 (2021).

54. Takahashi, K., Montecinos, A., Goubanova, K., & Dewitte, B. ENSO regimes: reinterpreting the canonical and Modoki El Niño. *Geophys. Res. Lett.* **38** (2011).

55. Wenzel, S., Cox, P. M., Eyring, V. & Friedlingstein, P. Emergent constraints on climate-carbon cycle feedbacks in the CMIP5 Earth system models. *J. Geophys. Res.: Biogeosci.* **119**, 794–807 (2014).

56. Cox, P. M. Emergent constraints on climate-carbon cycle feedbacks. *Curr. Clim. Change Rep.* **5**, 275–281 (2019).

57. Hall, A., Cox, P., Huntingford, C. & Klein, S. Progressing emergent constraints on future climate change. *Nat. Clim. Change* **9**, 269–278 (2019).

58. Kwiatkowski, L. et al. Emergent constraints on projections of declining primary production in the tropical oceans. *Nat. Clim. Change* **7**, 355–358 (2017).

59. Terhaar, J., Kwiatkowski, L. & Bopp, L. Emergent constraint on Arctic Ocean acidification in the twenty-first century. *Nature* **582**, 379–383 (2020).

60. Brient, F. Reducing uncertainties in climate projections with emergent constraints: concepts, examples and prospects. *Adv. Atmos. Sci.* **37**, 1–15 (2020).

61. Cocco, V. et al. Oxygen and indicators of stress for marine life in multi-model global warming projections. *Biogeosciences* **10**, 1849–1868 (2013).

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

63. Cai, W. et al. Climate impacts of the El Niño–Southern Oscillation on South America. *Nat. Rev. Earth Environ.* **1**, 215–231 (2020a).

64. Duteil, O., Böning, C. W. & Oschlies, A. Variability in subtropical–tropical cells drives oxygen levels in the tropical Pacific Ocean. *Geophys. Res. Lett.* **41**, 8926–8934 (2014).

65. Koutavas, A., Demenocal, P. B., Olive, G. C. & Lynch-Stieglitz, J. Mid-Holocene El Niño–Southern Oscillation (ENSO) attenuation revealed by individual foraminifera in eastern tropical Pacific sediments. *Geology* **34**, 993e996 (2006).

66. Mollier-Vogel, E. et al. Mid-Holocene deepening of the Southeast Pacific oxycline. *Glob. Planet. Change* **172**, 365–373 (2019).

67. Muñoz, P. et al. Reconstructing past variations in environmental conditions and paleoproductivity over the last ~8000 years off north- central Chile (30°S). *Biogeosciences* **17**, 5763–5785 (2020).

68. Salvatteci, R. et al. Centennial to millennial-scale changes in oxygenation and productivity in the Eastern Tropical South Pacific during the last 25000 years. *Quat. Sci. Rev.* **131**, 102–117 (2016).

69. Shin, N. Y. et al. More frequent central Pacific El Niño and stronger eastern pacific El Niño in a warmer climate. *npj Clim. Atmos. Sci.* **5**, 101 (2022).

70. Carréric, A. et al. Change in strong Eastern Pacific El Niño events dynamics in the warming climate. *Clim. Dyn.* **54**, 901–918 (2020).

71. Lopez, H., Lee, S. K., Kim, D., Wittenberg, A. T. & Yeh, S. W. Projections of faster onset and slower decay of El Niño in the 21st century. *Nat. Commun.* **13**, 1915 (2022).

72. Cai, W. et al. Butterfly effect and a self-modulating El Niño response to global warming. *Nature* **585**, 68–73 (2020b).

73. 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** (2009).

74. Frölicher, T. L. et al. Contrasting upper and deep ocean oxygen response to protracted global warming. *Glob. Biogeochem. Cycles* **34**, e2020GB006601 (2020).

75. Auderset, A. et al. Enhanced ocean oxygenation during Cenozoic warm periods. *Nature* **609**, 77–82 (2022).

76. Cardich, J. et al. Multidecadal changes in marine subsurface oxygenation off central Peru during the last ca. 170 years. *Front. Mar. Sci.* **6**, 270 (2019).

77. Glock, N., Erdem, Z. & Schönfeld, J. The Peruvian oxygen minimum zone was similar in extent but weaker during the Last Glacial Maximum than Late Holocene. *Commun. Earth Environ.* **3**, 307 (2022).

78. Hess, A. V. et al. A well-oxygenated eastern tropical Pacific during the warm Miocene. *Nature*, **619**, 521–525 (2023).

79. Moffitt, S. E. et al. Paleoceanographic insights on recent oxygen minimum zone expansion: lessons for modern oceanography. *PloS one* **10**, e0115246 (2015).
36. Sprintall, J., S. Cravatte, B. Dewitte, Y. Du and A. S. Gupta. Oceanic teleconnections, chapter 15 in "El Niño in a changing climate" AGU Book, ISBN: 978-1-119-54816-4, 528 (2020).

37. Köhn, E. E., Münnich, M., Vogt, M., Desmet, F. & Gruber, N. Strong habitat compression by extreme shoaling events of hypoxic waters in the Eastern Pacific. *J. Geophys. Res.: Oceans* **127**, e2022JC018429 (2022).

38. Leung, S., Thompson, L., McPhaden, M. J. & Mislan, K. A. S. ENSO drives near-surface oxygen and vertical habitat variability in the tropical Pacific. *Environ. Res. Lett.* **14**, 060402 (2019).

39. Busecke, J. J. M., Resplandy, L. & Dunne, J. P. The equatorial undercurrent and the oxygen minimum zone in the Pacific. *Geophys. Res. Lett.* **46**, 6716–6725 (2019).

40. Shigemitsu, M., Yamamoto, A., Oka, A. & Yamanaka, Y. One possible uncertainty in CMIP5 projections of low-oxygen water volume in the Eastern Tropical Pacific. *Glob. Biogeochem. Cycles* **31**, 804–820 (2017).

41. Ito, T. & Deutsch, C. Variability of the Oxygen Minimum Zone in the Tropical North Pacific during the late twentieth century. *Glob. Biogeochem. Cycles* **27**, 1119–1128 (2013).

42. José, Y. S., Stramma, L., Schmidtko, S., & Oschlies, A. (2019). ENSO-driven fluctuations in oxygen supply and vertical extent of oxygen-poor waters in the oxygen minimum zone of the Eastern Tropical South Pacific. *Biogeosci. Discuss.* 1–20 (2019).

43. Conejero, C., Dewitte, B., Garçon, V., Sudre, J. & Montes, I. ENSO diversity driving low-frequency change in mesoscale activity off Peru and Chile. *Sci. Rep.* **10**, 17902 (2020).

44. Espinoza-Morriberón, D. et al. Oxygen variability during ENSO in the tropical South Eastern Pacific. *Front. Mar. Sci.* **5**, 526 (2019).

45. Gutierrez, D. et al. Oxygenation episodes on the continental shelf of central Peru: remote forcing and benthic ecosystem response. *Prog. Oceanogr.* **79**, 177–189 (2008).

46. Mogollón, R. & Calil, P. H. R. On the effects of ENSO on ocean biogeochemistry in the Northern Humboldt Current System (NHCS): a modeling study. *J. Mar. Syst.* **172**, 137–159 (2017).

47. Duteil, O., Oschlies, A. & Böning, C. W. Pacific decadal oscillation and recent oxygen decline in the eastern tropical Pacific Ocean. *Biogeosciences (BG)* **15**, 7111–7126 (2018).

48. Pizarro, O., & Montecinos, A. Interdecadal variability of the thermocline along the west coast of South America. *Geophys. Res. Lett.* **31** (2004).

49. Oschlies, A. et al. Patterns of deoxygenation: sensitivity to natural and anthropogenic drivers. *Philos. Trans. R. Soc. A* **375**, 20160325 (2017).

50. Ito, T. & Deutsch, C. A conceptual model for the temporal spectrum of oceanic oxygen variability. *Geophys. Res. Lett.* **37**, L03601 (2010).

51. Capotondi, A., Wittenberg, A. T., Kug, J. S., Takahashi, K., & McPhaden, M. J. ENSO diversity. El Niño Southern Oscillation in a changing climate. **4**, 65–86 (2020).

52. Cai, W. et al. Increased variability of eastern Pacific El Niño under greenhouse warming. *Nature* **564**, 201–206 (2018).

53. Cai, W. et al. Changing El Niño–Southern Oscillation in a warming climate. *Nat. Rev. Earth Environ.* **2**, 628–644 (2021).

54. Takahashi, K., Montecinos, A., Goubanova, K., & Dewitte, B. ENSO regimes: reinterpreting the canonical and Modoki El Niño. *Geophys. Res. Lett.* **38** (2011).

55. Wenzel, S., Cox, P. M., Eyring, V. & Friedlingstein, P. Emergent constraints on climate-carbon cycle feedbacks in the CMIP5 Earth system models. *J. Geophys. Res.: Biogeosci.* **119**, 794–807 (2014).

56. Cox, P. M. Emergent constraints on climate-carbon cycle feedbacks. *Curr. Clim. Change Rep.* **5**, 275–281 (2019).

57. Hall, A., Cox, P., Huntingford, C. & Klein, S. Progressing emergent constraints on future climate change. *Nat. Clim. Change* **9**, 269–278 (2019).

58. Kwiatkowski, L. et al. Emergent constraints on projections of declining primary production in the tropical oceans. *Nat. Clim. Change* **7**, 355–358 (2017).

59. Terhaar, J., Kwiatkowski, L. & Bopp, L. Emergent constraint on Arctic Ocean acidification in the twenty-first century. *Nature* **582**, 379–383 (2020).

60. Brient, F. Reducing uncertainties in climate projections with emergent constraints: concepts, examples and prospects. *Adv. Atmos. Sci.* **37**, 1–15 (2020).

61. Cocco, V. et al. Oxygen and indicators of stress for marine life in multi-model global warming projections. *Biogeosciences* **10**, 1849–1868 (2013).

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

63. Cai, W. et al. Climate impacts of the El Niño–Southern Oscillation on South America. *Nat. Rev. Earth Environ.* **1**, 215–231 (2020a).

64. Duteil, O., Böning, C. W. & Oschlies, A. Variability in subtropical–tropical cells drives oxygen levels in the tropical Pacific Ocean. *Geophys. Res. Lett.* **41**, 8926–8934 (2014).

65. Koutavas, A., Demenocal, P. B., Olive, G. C. & Lynch-Stieglitz, J. Mid-Holocene El Niño–Southern Oscillation (ENSO) attenuation revealed by individual foraminifera in eastern tropical Pacific sediments. *Geology* **34**, 993e996 (2006).

66. Mollier-Vogel, E. et al. Mid-Holocene deepening of the Southeast Pacific oxycline. *Glob. Planet. Change* **172**, 365–373 (2019).

67. Muñoz, P. et al. Reconstructing past variations in environmental conditions and paleoproductivity over the last ~8000 years off north- central Chile (30°S). *Biogeosciences* **17**, 5763–5785 (2020).

68. Salvatteci, R. et al. Centennial to millennial-scale changes in oxygenation and productivity in the Eastern Tropical South Pacific during the last 25000 years. *Quat. Sci. Rev.* **131**, 102–117 (2016).

69. Shin, N. Y. et al. More frequent central Pacific El Niño and stronger eastern pacific El Niño in a warmer climate. *npj Clim. Atmos. Sci.* **5**, 101 (2022).

70. Carréric, A. et al. Change in strong Eastern Pacific El Niño events dynamics in the warming climate. *Clim. Dyn.* **54**, 901–918 (2020).

71. Lopez, H., Lee, S. K., Kim, D., Wittenberg, A. T. & Yeh, S. W. Projections of faster onset and slower decay of El Niño in the 21st century. *Nat. Commun.* **13**, 1915 (2022).

72. Cai, W. et al. Butterfly effect and a self-modulating El Niño response to global warming. *Nature* **585**, 68–73 (2020b).

73. 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** (2009).

74. Frölicher, T. L. et al. Contrasting upper and deep ocean oxygen response to protracted global warming. *Glob. Biogeochem. Cycles* **34**, e2020GB006601 (2020).

75. Auderset, A. et al. Enhanced ocean oxygenation during Cenozoic warm periods. *Nature* **609**, 77–82 (2022).

76. Cardich, J. et al. Multidecadal changes in marine subsurface oxygenation off central Peru during the last ca. 170 years. *Front. Mar. Sci.* **6**, 270 (2019).

77. Glock, N., Erdem, Z. & Schönfeld, J. The Peruvian oxygen minimum zone was similar in extent but weaker during the Last Glacial Maximum than Late Holocene. *Commun. Earth Environ.* **3**, 307 (2022).

78. Hess, A. V. et al. A well-oxygenated eastern tropical Pacific during the warm Miocene. *Nature*, **619**, 521–525 (2023).

79. Moffitt, S. E. et al. Paleooceanographic insights on recent oxygen minimum zone expansion: lessons for modern oceanography. *PloS one* **10**, e0115246 (2015).
80. Galeotti, S. et al. Evidence for active El Niño Southern Oscillation variability in the Late Miocene greenhouse climate. *Geology* **38**, 419–422 (2010).

81. Cobb, K. M. et al. Highly variable El Niño-Southern Oscillation throughout the holocene. *Science* **339**, 67–70 (2013).

82. McGregor, S., Timmermann, A., England, M. H., Timm, O. E. & Wittenberg, A. T. Inferred changes in El Niño-Southern oscillation variance over the past six centuries. *Clim* **9**, 2269–2284 (2013).

83. Séférían, R. et al. Tracking improvement in simulated marine biogeochemistry between CMIP5 and CMIP6. *Curr. Clim. Change Rep.* **6**, 95–119 (2020).

84. Payne, M. R. et al. Uncertainties in projecting climate-change impacts in marine ecosystems. *ICES J. Mar. Sci.* **73**, 1272–1282 (2016).

85. Gruber, N. et al. Eddy-induced reduction of biological production in eastern boundary upwelling systems. *Nat. Geosci.* **4**, 787–792 (2011).

86. Vergara, O. et al. Seasonal variability of the Oxygen Minimum Zone off Peru in a high-resolution regional coupled model. *Biogeosciences* **13**, 4389–4410 (2016).

87. Bettencourt J. et al. Boundaries of the Oxygen Minimum Zone shaped by coherent mesoscale dynamics. *Nat. Geosci.* https://doi.org/10.1038/NGEO2570. (2015).

88. Pizarro-Koch, M. et al. On the interpretation of changes in OMZ volume off Central Chile during two La Niña events (2001 and 2007). *Front. Mar. Sci.* https://doi.org/10.3389/fmars.2023.1155932 (2023).

89. Thomsen, S. et al. Do submesoscale frontal processes ventilate the oxygen minimum zone off Peru? *Geophys. Res. Lett.* **43**, 8133–8142 (2016).

90. Xiu, P. & Chai, F. Eddies affect subsurface phytoplankton and oxygen distributions in the North Pacific subtropical gyre. *Geophys. Res. Lett.* **47**, e2020GL087037 (2020).

91. Chen, C., Cane, M. A., Wittenberg, A. T. & Chen, D. ENSO in the CMIP5 simulations: life cycles, diversity, and responses to climate change. *J. Clim.* **30**, 775–801 (2017).

92. Long, A. M., Jurgensen, S. K., Petchel, A. R., Savoie, E. R., & Brum, J. R. Microbial ecology of oxygen minimum zones amidst ocean deoxygenation. *Front. Microbiol.* **12**, 748961 (2021).

93. Ito, T. Optimal interpolation of global dissolved oxygen: 1965–2015. *Geosci. Data J.* https://doi.org/10.1002/gdj3.130 (2021).

94. Kessler, W.S., S. Cravatte, & Lead Authors. Final Report of TPOS 2020. GOOS-268, 83 pp. [Available online at https://tropicalocean.org/tpos2020-project-archive/reports/](https://tropicalocean.org/tpos2020-project-archive/reports/) (2021).

95. Garçon, V. et al. Multidisciplinary observing in the world ocean's oxygen minimum zone regions: from climate to fish—the VOICE Initiative. *Front. Mar. Sci.* **6**, 722 (2019).

96. Smith, N. et al. Tropical pacific observing system. *Front. Mar. Sci.* https://doi.org/10.3389/fmars.2019.00031. (2019).

97. Stammer, S. et al. Ocean climate observing requirements in support of climate research and climate information. *Front. Mar. Sci.* https://doi.org/10.3389/fmars.2019.00444. (2019).

98. Ridgway, K. R., Dunn, J. R. & Wilkin, J. L. Ocean interpolation by four-dimensional least squares -Application to the waters around Australia. *J. Atmos. Ocean. Tech.* **19**, 1357–1375 (2002).

99. Rayner, N. A. et al. Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century. *J. Geophys. Res.* **108**, 4407 (2003).

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

101. 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).

102. Kay, J. E. et al. The community earth system model (CESM) large ensemble project: a community resource for studying climate change in the presence of internal climate variability. *Bull. Am. Meteorol. Soc.* **96**, 1333–1349 (2015).

103. O'Neill, B. C. et al. The roads ahead: Narratives for shared socioeconomic pathways describing world futures in the 21st century. *Glob. Environ. Change* **42**, 169–180 (2017).

104. Gutknecht, E. et al. Coupled physical/biogeochemical modeling including O2-dependent processes in the Eastern Boundary Upwelling Systems: application in the Benguela. *Biogeosciences* **10**, 3559–3591 (2013a).

105. Gutknecht, E. et al. Nitrogen transfers off Walvis Bay: a 3-D coupled physical/biogeochemical modeling approach in the Namibian upwelling system. *Biogeosciences* **10**, 4117–4135 (2013b).

106. Pizarro-Koch M. et al. Seasonal variability of the southern tip of the Oxygen Minimum Zone in the Eastern South Pacific (30°–38°S): a modeling study. *J. Geophys. Res.-Oceans.* **124**, 8574–8604 (2019).

107. Karamperidou, C., Jin, F.-F. & Conroy, J. L. The importance of ENSO nonlinearities in tropical pacific response to external forcing. *Clim. Dyn.* **49**, 2695–2704 (2017).

108. White, H. G. Skewness, kurtosis and extreme values of Northern Hemisphere geopotential heights. *Mon. Weather Rev.* **108**, 1446–1455 (1980).

109. 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).

110. Matveeva, T., Gushchina, D. & Dewitte, B. The seasonal relationship between intraseasonal tropical variability and ENSO in CMIP5. *Geosci. Model Dev.* **11**, 2373–2392 (2018).

## Acknowledgements

B.D. acknowledges support from ANID (Concurso de Fortalecimiento al Desarrollo Científico de Centros Regionales 2020-R20F0008-CEAZA, Anillo Eclipse ACT210071 and Centro de Investigación Oceanográfica en el Pacífico Sur-Oriental COPAS COASTAL FB210021, Fondecyt Regular no. 1190276 and no. 1231174). B.D. and V.G. are also supported by the EU H2020 FutureMares project (Theme LC-CLA-06-2019, Grant agreement no. 869300) and the CE2COAST project funded by ANR (FR), BELSPO (BE), FCT (PT), IZM (LV), MI (IE), MIUR (IT), Rannis (IS) and RCN (NO) through the 2019 "Joint Transnational Call on Next Generation Climate Science in Europe for Oceans" initiated by JPI Climate and JPI Oceans. M.R. acknowledges supports from ANID through the Anillo BiodUCCT (ATE220044). O.D. acknowledges support from the DFG (Deutsch Forschungsgemeinschaft) (project "SyVarBio" - 434479332). This work was carried out as part of I.A.'s Master's thesis in Geophysics (University of Concepcion). We acknowledge the World Climate Research Program's Working Group on Coupled Modeling, which is responsible for CMIP, and we thank the climate modeling groups (listed in Table S1 of this paper) for producing and making available their model output. For CMIP the U.S. Department of Energy's Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. We would like to thank the two anonymous reviewers for their constructive comments, which helped improve the presentation of the results.

## Author contributions

B.D. and V.G. conceived the study. I.A. performed the model output analysis and produced the figures, with help from B.D. B.D. wrote the initial manuscript in discussion with I.A., P.M., C.P., and V.G. G.I.M., O.D., A.P., O.P., and A.O. contributed with expertise to the writing. M.R. provided observations, and W.K. provided the CMIP5 model data. All authors commented on and revised the manuscript.

## Competing interests

The authors declare no competing interests.
80. Galeotti, S. et al. Evidence for active El Niño Southern Oscillation variability in the Late Miocene greenhouse climate. *Geology* **38**, 419–422 (2010).

81. Cobb, K. M. et al. Highly variable El Niño-Southern Oscillation throughout the holocene. *Science* **339**, 67–70 (2013).

82. McGregor, S., Timmermann, A., England, M. H., Timm, O. E. & Wittenberg, A. T. Inferred changes in El Niño-Southern oscillation variance over the past six centuries. *Clim* **9**, 2269–2284 (2013).

83. Séférían, R. et al. Tracking improvement in simulated marine biogeochemistry between CMIP5 and CMIP6. *Curr. Clim. Change Rep.* **6**, 95–119 (2020).

84. Payne, M. R. et al. Uncertainties in projecting climate-change impacts in marine ecosystems. *ICES J. Mar. Sci.* **73**, 1272–1282 (2016).

85. Gruber, N. et al. Eddy-induced reduction of biological production in eastern boundary upwelling systems. *Nat. Geosci.* **4**, 787–792 (2011).

86. Vergara, O. et al. Seasonal variability of the Oxygen Minimum Zone off Peru in a high-resolution regional coupled model. *Biogeosciences* **13**, 4389–4410 (2016).

87. Bettencourt J. et al. Boundaries of the Oxygen Minimum Zone shaped by coherent mesoscale dynamics. *Nat. Geosci.* https://doi.org/10.1038/NGEO2570. (2015).

88. Pizarro-Koch, M. et al. On the interpretation of changes in OMZ volume off Central Chile during two La Niña events (2001 and 2007). *Front. Mar. Sci.* https://doi.org/10.3389/fmars.2023.1155932 (2023).

89. Thomsen, S. et al. Do submesoscale frontal processes ventilate the oxygen minimum zone off Peru? *Geophys. Res. Lett.* **43**, 8133–8142 (2016).

90. Xiu, P. & Chai, F. Eddies affect subsurface phytoplankton and oxygen distributions in the North Pacific subtropical gyre. *Geophys. Res. Lett.* **47**, e2020GL087037 (2020).

91. Chen, C., Cane, M. A., Wittenberg, A. T. & Chen, D. ENSO in the CMIP5 simulations: life cycles, diversity, and responses to climate change. *J. Clim.* **30**, 775–801 (2017).

92. Long, A. M., Jurgensen, S. K., Petchel, A. R., Savoie, E. R., & Brum, J. R. Microbial ecology of oxygen minimum zones amidst ocean deoxygenation. *Front. Microbiol.* **12**, 748961 (2021).

93. Ito, T. Optimal interpolation of global dissolved oxygen: 1965–2015. *Geosci. Data J.* https://doi.org/10.1002/gdj3.130 (2021).

94. Kessler, W.S., S. Cravatte, & Lead Authors. Final Report of TPOS 2020. GOOS-268, 83 pp. [Available online at https://tropicalocean.org/tpos2020-project-archive/reports/](https://tropicalocean.org/tpos2020-project-archive/reports/) (2021).

95. Garçon, V. et al. Multidisciplinary observing in the world ocean's oxygen minimum zone regions: from climate to fish—the VOICE Initiative. *Front. Mar. Sci.* **6**, 722 (2019).

96. Smith, N. et al. Tropical pacific observing system. *Front. Mar. Sci.* https://doi.org/10.3389/fmars.2019.00031. (2019).

97. Stammer, S. et al. Ocean climate observing requirements in support of climate research and climate information. *Front. Mar. Sci.* https://doi.org/10.3389/fmars.2019.00444. (2019).

98. Ridgway, K. R., Dunn, J. R. & Wilkin, J. L. Ocean interpolation by four-dimensional least squares -Application to the waters around Australia. *J. Atmos. Ocean. Tech.* **19**, 1357–1375 (2002).

99. Rayner, N. A. et al. Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century. *J. Geophys. Res.* **108**, 4407 (2003).

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

101. 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).

102. Kay, J. E. et al. The community earth system model (CESM) large ensemble project: a community resource for studying climate change in the presence of internal climate variability. *Bull. Am. Meteorol. Soc.* **96**, 1333–1349 (2015).

103. O'Neill, B. C. et al. The roads ahead: Narratives for shared socioeconomic pathways describing world futures in the 21st century. *Glob. Environ. Change* **42**, 169–180 (2017).

104. Gutknecht, E. et al. Coupled physical/biogeochemical modeling including O2-dependent processes in the Eastern Boundary Upwelling Systems: application in the Benguela. *Biogeosciences* **10**, 3559–3591 (2013a).

105. Gutknecht, E. et al. Nitrogen transfers off Walvis Bay: a 3-D coupled physical/biogeochemical modeling approach in the Namibian upwelling system. *Biogeosciences* **10**, 4117–4135 (2013b).

106. Pizarro-Koch M. et al. Seasonal variability of the southern tip of the Oxygen Minimum Zone in the Eastern South Pacific (30°–38°S): a modeling study. *J. Geophys. Res.-Oceans.* **124**, 8574–8604 (2019).

107. Karamperidou, C., Jin, F.-F. & Conroy, J. L. The importance of ENSO nonlinearities in tropical pacific response to external forcing. *Clim. Dyn.* **49**, 2695–2704 (2017).

108. White, H. G. Skewness, kurtosis and extreme values of Northern Hemisphere geopotential heights. *Mon. Weather Rev.* **108**, 1446–1455 (1980).

109. 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).

110. Matveeva, T., Gushchina, D. & Dewitte, B. The seasonal relationship between intraseasonal tropical variability and ENSO in CMIP5. *Geosci. Model Dev.* **11**, 2373–2392 (2018).

## Acknowledgements

B.D. acknowledges support from ANID (Concurso de Fortalecimiento al Desarrollo Científico de Centros Regionales 2020-R20F0008-CEAZA, Anillo Eclipse ACT210071 and Centro de Investigación Oceanográfica en el Pacífico Sur-Oriental COPAS COASTAL FB210021, Fondecyt Regular no. 1190276 and no. 1231174). B.D. and V.G. are also supported by the EU H2020 FutureMares project (Theme LC-CLA-06-2019, Grant agreement no. 869300) and the CE2COAST project funded by ANR (FR), BELSPO (BE), FCT (PT), IZM (LV), MI (IE), MIUR (IT), Rannis (IS) and RCN (NO) through the 2019 "Joint Transnational Call on Next Generation Climate Science in Europe for Oceans" initiated by JPI Climate and JPI Oceans. M.R. acknowledges supports from ANID through the Anillo BiodUCCT (ATE220044). O.D. acknowledges support from the DFG (Deutsch Forschungsgemeinschaft) (project "SyVarBio" - 434479332). This work was carried out as part of I.A.'s Master's thesis in Geophysics (University of Concepcion). We acknowledge the World Climate Research Program's Working Group on Coupled Modeling, which is responsible for CMIP, and we thank the climate modeling groups (listed in Table S1 of this paper) for producing and making available their model output. For CMIP the U.S. Department of Energy's Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. We would like to thank the two anonymous reviewers for their constructive comments, which helped improve the presentation of the results.

## Author contributions

B.D. and V.G. conceived the study. I.A. performed the model output analysis and produced the figures, with help from B.D. B.D. wrote the initial manuscript in discussion with I.A., P.M., C.P., and V.G. G.I.M., O.D., A.P., O.P., and A.O. contributed with expertise to the writing. M.R. provided observations, and W.K. provided the CMIP5 model data. All authors commented on and revised the manuscript.

## 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-024-01427-2.

**Correspondence** and requests for materials should be addressed to Boris Dewitte.

**Peer review information** *Communications Earth & Environment* thanks the anonymous reviewers for their contribution to the peer review of this work. Primary Handling Editors: Michael Stukel, Clare Davis and Aliénor Lavergne. A peer review file is available

**Reprints and permissions 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.

**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) 2024
**Additional information**

**Supplementary information** The online version contains supplementary material available at
https://doi.org/10.1038/s43247-024-01427-2.

**Correspondence** and requests for materials should be addressed to Boris Dewitte.

**Peer review information** *Communications Earth & Environment* thanks the anonymous reviewers for their contribution to the peer review of this work. Primary Handling Editors: Michael Stukel, Clare Davis and Aliénor Lavergne. A peer review file is available

**Reprints and permissions 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.

**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) 2024