{'fraction_non_alphanumeric': 0.03282752080010346, 'fraction_numerical': 0.02819330085786955, 'mean_word_length': 4.7712402040054736, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Contrasting Upper and Deep Ocean Oxygen Response to Protracted Global Warming

It is well established that the ocean is currently losing dissolved oxygen (O2) in response to ocean warming, but the long‐term, equilibrium response of O2 to a warmer climate is neither well quantified nor understood. Here we use idealized multimillennial global warming simulations with a comprehensive Earth system model to show that the equilibrium response in ocean O2 differs fundamentally from the ongoing transient response. After physical equilibration of the model (>4,000 years) under a two times preindustrial CO2 scenario, the deep ocean is better ventilated and oxygenated compared to preindustrial conditions, even though the deep ocean is substantially warmer. The recovery and overshoot of deep convection in the Weddell Sea and especially the Ross Sea after ~720 years causes a strong increase in deep ocean O2 that overcompensates the solubility‐driven decrease in O2. In contrast, O2 in most of the upper tropical ocean is substantially depleted owing to the warming‐induced O2 decrease dominating over changes in ventilation and biology. Our results emphasize the millennial‐scale impact of global warming on marine life, with some impacts emerging many centuries or even millennia after atmospheric CO2 has stabilized.

Global warming is the main driver of open ocean deoxygenation, but it acts in two different ways. First, the solubility of O 2 decreases in warmer waters (Garcia & Gordon, 1992). Second, the warming of the upper ocean combined with the climate change-induced freshening of the subpolar regions leads to a general strengthening of upper ocean stratification (Collins et al., 2013;Manabe et al., 1991;Sarmiento et al., 1998). This reduces the downward supply of O 2 -rich surface waters (Breitburg et al., 2018;Keeling et al., 2010;Levin, 2018;Oschlies et al., 2018), depleting O 2 in the ocean's interior. The contribution of changes in the biological consumption of O 2 in the ocean's interior tends to differ between the low latitudes, where the enhanced upper ocean stratification reduces the resupply of nutrients to the productive surface ocean curtailing primary and export production (Laufkötter et al., 2015(Laufkötter et al., , 2016Steinacher et al., 2010), and the high latitudes where most models predict an increase in export and hence O 2 consumption at depth. These three aspects reinforce each other in the temperate to high latitudes, leading to very consistently projected decreases in O 2 under transient global warming. In the upper ocean of the low latitudes, the changes in ventilation/biology and O 2 solubility tend to compensate each other, resulting in smaller projected O 2 decreases (Bopp et al., 2017;Gnanadesikan et al., 2012) or even increases.
Only a few studies have documented ocean O 2 changes beyond 2100. Recently, Fu et al. (2018) suggested that the current trend of gradually declining O 2 levels in the tropical thermocline may reverse after 2100, resulting in shrinking volumes of the OMZs through the end of their simulation in the Year 2300. They attributed this potential reoxygenation of the tropical ocean thermocline to intensified ventilation that kicked in their simulation after 2150 and to a decline in tropical export production leading to reduced biological O 2 consumption in the underlying thermocline. In the deep ocean, several studies suggested a recovery of deep ocean O 2 levels under multimillennial global warming scenarios (Battaglia & Joos, 2018;Ridgwell & Schmidt, 2010;Schmittner et al., 2008;Yamamoto et al., 2015). In all of these models, the recovery was caused by a millennial-scale emergence of a more vigorous deep ocean ventilation under protracted global warming. However, all but one of these studies were conducted with reduced complexity models that have no dynamic atmosphere or employ an ocean model with very coarse ocean resolution. This is primarily a consequence of comprehensive ESMs generally being too expensive to run for millennia. And the one exception, that is, the study by Yamamoto et al. (2015), employed an offline approach with circulation fields from a CMIP3-class climate model. It is thus imperative to revisit the deep ocean reoxygenation with a state-of-theart ESM of the CMIP5 class and to investigate also the development of O 2 in the upper ocean, especially in the tropical thermocline. There, the evolution of the O 2 -depleted water masses depends on the subtle balance between changes in O 2 -solubility and ocean ventilation (including O 2 consumption) Ito & Deutsch, 2013). It remains, therefore, unclear whether these typically opposing processes in the low latitudes will evolve simultaneously on multimillennial timescales, when the climate system approaches a new equilibrium under elevated atmospheric CO 2 concentrations.
Here we use millennial-scale simulations of a comprehensive ESM under an idealized global warming scenario to show that the deep ocean becomes better ventilated and oxygenated and that, conversely, the (sub) tropical shallow subsurface ocean is projected to lose O 2 in an equilibrated warm climate.

GFDL-ESM2M
The Geophysical Fluid Dynamics Laboratory (GFDL)-ESM2M is a comprehensive global ESM developed at the GFDL . The GFDL-ESM2M consists of atmosphere, ocean, sea ice, and land modules including biogeochemistry. The ocean component of the GFDL-ESM2M is the Modular Ocean Model version 4p1 (MOM4p1; Griffies, 2012), with a nominal 1°horizontal resolution with higher resolution near the equator and 50 unevenly spaced vertical depth levels. Ocean mesoscale eddies are not explicitly resolved but are parameterized using the GM-Redi schemes (Griffies, 1998). The atmospheric module (AM2; Anderson et al., 2004) has a horizontal resolution of 2°latitude × 2.5°longitude and 24 vertical levels. In contrast to earlier multimillennial simulations focusing on oceanic O 2 and using models of intermediate complexity (Battaglia & Joos, 2018;Schmittner et al., 2008), the GFDL-ESM2M includes changes in the hydrological cycle, cloud cover, wind fields, and sea ice (Winton, 2000). Similar to other CMIP5-type ESMs, the GFDL-ESM2M lacks a land ice sheet model and, therefore, freshwater fluxes from the Greenland and Antarctic ice sheets as well as mountain glaciers.
The ocean biogeochemical and ecological component is version two of the Tracers of Ocean Productivity with Allometric Zooplankton (TOPAZv2) module that parameterizes the cycling of carbon, nitrogen, phosphorus, silicon, iron, oxygen, alkalinity, and lithogenic material as well as surface sediment calcite . TOPAZv2 includes three explicit phytoplankton groups, namely, small and large phytoplankton as well as diazotrophs, and one implicit zooplankton group. A fraction of the sinking organic material is associated with denser ballast material (biogenic opal, calcium carbonate, and lithogenic material) and remineralizes only when these materials dissolve. At well oxygenated conditions, the unprotected fraction is remineralized with a constant depth scale of 188 m. Remineralization is slower at low O 2 and suppressed in the absence of O 2 (<1 μmol kg −1 ) (Gnanadesikan et al., 2012). The remineralization and sinking speed of organic material is insensitive to temperature.

Global Biogeochemical Cycles
FRÖLICHER ET AL.
The GFDL-ESM2M simulates the large-scale physical  and biogeochemical  patterns with competitive fidelity among the CMIP5 suite of models (Bopp et al., 2013). Analyzing a 1,500 year preindustrial control simulation (see section 2.2), the model reproduces observation-based estimates of the formation and circulation of North Atlantic Deep Water . But the model forms Antarctic Bottom Water by open ocean convection (Heuzé et al., 2013), while in reality bottom water is generally formed as dense water spills off the Antarctic continental shelf. This feature is a common bias in current ESMs, whose implications will be discussed in our caveat section below. The simulated spatial distribution of O 2 matches the observations well, both in the upper and the deeper ocean (supporting information Figure S1). The mean O 2 concentration simulated by the GFDL-ESM2M in the upper (sub)tropical ocean (30°S-30°N) between 200 and 1,000 m of 95 mmol m −3 is somewhat lower than the observed concentration of 113 mmol m −3 . This negative bias in O 2 is also reflected in the volume of O 2 depleted waters. For a threshold of 20 mmol m −3 , the volume of O 2 -depleted waters in GFDL-ESM2M is by a factor of 3.8 higher than the observational-based estimate (Bianchi et al., 2012;Fu et al., 2018). In the global deep ocean below 2,000 m, the model simulates also a slightly lower O 2 concentrations of 174 mmol m −3 compared to observations (189 mmol m −3 ), with most of the negative bias located in the eastern Pacific and Atlantic (supporting information Figure S1b). Also taking into account previous assessments of simulated O 2 by the GFDL-ESM2M (Bopp et al., 2013;Cabré et al., 2015;Dunne et al., 2013;Fu et al., 2018;Gnanadesikan et al., 2012), we conclude that the model performs relatively well against some key performance metrics and, therefore, appears to be adequate to assess changes in open ocean oxygen on multicentennial to millennial timescales.

Multimillennial Global Warming Simulation
We analyze an idealized 4,500 year-long CO 2 doubling simulation (Paynter et al., 2018) initialized from a multimillennial 1860 preindustrial simulation, wherein atmospheric CO 2 was held fixed at 287.1 ppm. From this value, atmospheric CO 2 was then prescribed to increase by 1% per year during a period of 70 years until it reached a CO 2 doubling (574.2 ppm). Thereafter, atmospheric CO 2 is kept fixed (Figure 1a). All radiative forcing agents other than CO 2 were kept constant at their preindustrial levels. We also performed a second ensemble member for the period ranging from Years 673 to 822 branching off the first ensemble simulation to assess whether the onset of the deepwater formation in the Ross Sea was stochastically driven. In addition, an accompanying control simulation with constant preindustrial radiative forcing was run for the first 1,500 years in order to determine the influence of model drift.
The general climate response of the GFDL-ESM2M to a doubling of atmospheric CO 2 is described in detail elsewhere (Paynter et al., 2018;Rugenstein et al., 2020). Here we nonetheless provide a brief summary of the most important physical characteristics. During the first 70 years (average period 61 to 80 years), the global mean surface air temperature increases by 1.3°C (blue line in Figure 1b). This transient climate response is at the lower end of the CMIP5 model range of 1.1°C to 2.5°C (Forster et al., 2013). The initial rate of surface air  (Savitzky & Golay, 1964). Note that the x axis is discontinuous, with finer temporal resolution of the first 70 and 720 years.
temperature increase poleward of 55°S is slow, as it is dampened by the large heat capacity of the ocean and the continuing Southern Ocean heat uptake, but maintains a higher warming rate than the global mean after CO 2 is stabilized in Year 70 and is abruptly increased by an additional 0.7°C around Year 720 (orange line in Figure 1b). The equilibrium temperature response of 3.4°C to a doubling of atmospheric CO 2 is similar to the CMIP5 multimodel mean of 3.2°C (Forster et al., 2013) and to the last deglacial global average temperature change (Shakun et al., 2012). The physical model can be considered to be in equilibrium with the elevated radiative forcing at the end of the 4,500 year simulation, that is, the global ocean heat uptake is virtually zero (−0.07 Wm −2 averaged over last 200 years of the simulation; Paynter et al., 2018). The surface ocean temperature change equilibrates at 2.6°C (green line in Figure 1b) and the global ocean temperature change at 1.7°C (red line in Figure 1b).

Analysis
To assess changes in ocean O 2 , we decompose the total changes in oxygen ( (1) ΔO 2 sol is determined by the O 2 solubility in seawater, which is a function of ocean temperature and to a lesser extent salinity (Garcia & Gordon, 1992). Warming and salinification of the ocean decreases the solubility and therefore ΔO 2 sol . The term ΔO 2 circ/rem reflects changes in organic matter remineralization and physical O 2 supply by the large-scale ocean circulation and mixing. A decrease in biological export production would lead to an increase in ΔO 2 circ/rem owing to reduced O 2 demand from remineralization. In contrast, a decrease in ventilation (increase in stratification) would lead to a decrease in ΔO 2 circ/rem , as the longer residence time increases cumulative O 2 demand in subsurface waters. ΔO 2 diseq is the disequilibrium of the surface ocean O 2 with respect to the atmospheric O 2 concentration, as surface ocean O 2 is usually not at saturation. In most of the low-latitude surface ocean, O 2 is slightly supersaturated due to photosynthesis and sluggish ventilation. In the high latitudes surface waters, in particular in regions where convection occurs, O 2 can be undersaturated due to rapid heat loss and mixing of O 2 -depleted deep waters to the surface (Ito et al., 2004). In this study, O 2 diseq is included in ΔO 2 circ/rem/diseq : ΔO 2 circ/rem/diseq can therefore be calculated as the difference between ΔO 2 and ΔO 2 sol . We also estimate the ventilation time and the organic matter respiration. For the ventilation time, we use an ideal age tracer, which is set to zero in the surface ocean and ages at a rate 1 year year −1 below that (Thiele & Sarmiento, 1990). We use the changes in the export production of particulate organic matter at 100 m depth as a proxy for organic matter remineralization at subsurface. A quantitative analysis of all O 2 transport terms, as has been done previously (Gnanadesikan et al., 2012), is not possible, because of a lack of diagnostics in these multimillennial simulations.
We calculate the Atlantic meridional overturning circulation (AMOC) and the Southern Ocean meridional overturning circulation (SMOC) as the vertically and zonally integrated meridional transport in the Atlantic basin and in the Southern Ocean, respectively: where ν is the meridional velocity expressed as the sum of resolved and parameterized mesoscale and submesoscale processes (Griffies, 2012), x is the zonal coordinate (increasing eastward), z is the vertical coordinate (increasing upward), η is the ocean free surface, and x w and x e are the westward and eastward boundaries of the North Atlantic basin or the Southern Ocean region. We report the strength of the AMOC and SMOC as the annual mean below 500 m of the vertical maximum north of 20°N for the AMOC and the minimum south of 60°S for the SMOC.

Global Biogeochemical Cycles
FRÖLICHER ET AL.
We assess the effect of ocean-atmosphere interactions on surface density changes by calculating the sum of the thermal (F ρ,thermal ) and haline (F ρ , haline ) density fluxes according to Liu et al. (2017): where α and β are the thermal expansion and haline contraction coefficients, respectively. SST is the sea surface temperature and SSS the sea surface salinity. Q denotes the net heat flux from the atmosphere into the ocean, c p the specific heat capacity of seawater, ρ(0, SST) is the density of freshwater with salinity of 0 and temperature of SST, and E, P, R, and I the freshwater flux from evaporation, precipitation, river runoff, and sea ice melt, respectively.
To investigate changes in the stability of the water column (i.e., stratification), we calculate the squared Brunt-Väisälä frequency (N 2 ; Talley et al. (2011)): where g it the local acceleration of gravity, ρ the in situ density, and dρ/dz the vertical gradient of the in situ density. To assess whether stability changes are caused by salinity or temperature changes, we also determine N 2 based on the density with either the temperature (N 2 (ρ T0,S )) or the salinity (N 2 (ρ T,S0 )) prescribed to preindustrial (0) values.
For all analyses, we use the model output as is, that is, we do not correct the output with the results from the preindustrial control simulation. However, we generally illustrate the corresponding preindustrial control simulation to indicate potential drift.

Global Ocean Oxygen Changes
The multimillennial response of the global oceanic O 2 content to changes in atmospheric CO 2 and associated global warming is characterized by three phases (Figure 2a). Phase I covers Years 1-70 and is characterized by steadily increasing atmospheric CO 2 and globally decreasing O 2 levels; phase II encompasses Years 71-720 and is characterized by constant atmospheric CO 2 levels, yet global O 2 levels continue to decline; phase III covers Years 721-4500 and is marked by constant atmospheric CO 2 and a gradually recovering global O 2 inventory.
During phase I, under steadily increasing atmospheric CO 2 forcing, the global oceanic O 2 inventory decreases by −3.6 Pmol (−1.6%; black line in Figure 2a). After the CO 2 forcing has stabilized during phase II, the global O 2 inventory further decreases such that the O 2 content is 17.5 Pmol (−7.8%) smaller by Year 720 than at the beginning of the simulation. The warming-induced reduction in solubility, O 2 sol , accounts for 40% of the total decrease by the end of phase I, amounting 68% of the total O 2 loss by the end of phase II (blue line in Figure 2a). The additional O 2 loss results from changes in ocean circulation and organic matter respiration as indicated by the decrease in O 2 circ/rem/diseq (orange line in Figure 2a). The relative contribution of O 2 sol to the total change during phase I is consistent with earlier multimodel studies, which reported that 18-50% of the O 2 depletion could be attributed to warming over the 21st century (Bopp et al., 2013). It is also consistent with detailed analysis of observations for the relatively well observed North Atlantic for the past 50 years (Stendardo & Gruber, 2012), whereas the global analyses by

Upper Ocean Deoxygenation in the Tropics
In the (sub)tropical upper (30°N-30°S, 200-1,000 m) ocean, where O 2 is generally low, the O 2 content decreases by 0.8 Pmol during the first two phases, that is, up to 720 years into simulation. It recovers only marginally during phase III (black line in Figure 2c), such that at the end of the simulation, the O 2 content is still 0.7 Pmol lower than at the beginning. As is the case for the global ocean O 2 content, most (77%) of the total O 2 decrease occurs during phase II, when atmospheric CO 2 has already stabilized. During phases I and II, the decrease in O 2 sol (blue line in Figure 2c) dominates the decrease in O 2 . O 2 circ/rem/diseq increases (orange line in Figure 2c), but the increase is smaller in magnitude than the decrease

10.1029/2020GB006601
Global Biogeochemical Cycles increases strongly during phase III around Antarctica (see section 3.3) and also, yet to a lesser extent, in the Indian Ocean.
In the midlatitudes of the South Pacific, the increase in O 2 circ/rem/diseq is larger than the decrease in O 2 sol leading to the overall increase in O 2 (Figure 4). This is not the case in the tropical Atlantic, where changes in O 2 sol are larger than changes in O 2 circ/rem/diseq , and therefore, O 2 decreases. In general, changes in O 2 circ/rem/diseq largely explain the pattern of the O 2 changes.
The changes in the volume of the OMZs ( Figure 5) mirror to some extent the changes in the upper ocean O 2 concentrations described above, irrespective of whether these regions are defined by a threshold of 50 or 80 mmol m −3 . OMZs expand substantially during phase I and especially during phase II. The near-anoxic regions (O 2 < 5 mmol m −3 ) nearly double in volume (+16·10 15 m 3 or +44%). During phase III, the OMZs start to shrink but remain larger by the end of the simulation compared to preindustrial. Note that the volume of the OMZs increases over the course of the preindustrial control simulation (dashed lines in Figure 5), but this increase alone cannot explain the overall expansion of the OMZs in the 4,500 year-long simulation.

Deep Ocean Deoxygenation and Reoxygenation
In the deep ocean below 2,000 m (black line in Figure 2d), the dissolved O 2 inventory evolves fundamentally differently. The O 2 inventory decreases by 13 Pmol during phases I and II. However, the deep O 2 content recovers during phase III and overshoots preindustrial levels by the end of the simulation (+5 Pmol O 2 ) indicating that the recovery and overshoot of the global O 2 inventory (black line in Figure 2a (Figures 6h and 6i). During phase III, O 2 circ/rem/diseq also strongly increases in the Ross Sea, along the Chilean Margin and the North Pacific, in part opposite to the changes simulated during phase II. The large O 2 decrease in the equatorial Pacific during phase II is also mainly driven by a decrease in O 2 circ/rem/diseq . The simulated regional changes in O 2 sol are generally negative (Figures 6d-6f), yet more uniform and subdued than the changes in O 2 circ/rem/diseq .  Figure 7a) and net primary and export production (Figure 7b). We first discuss the changes in the upper ocean of the (sub)tropics and then in the deep ocean.

Processes Controlling Changes in O
The increase in O 2 circ/rem/diseq in the upper ocean of the (sub)tropics during phase I and in the early part of phase II is associated with a large decrease (−70 years) in ventilation age-the average time required for

10.1029/2020GB006601
Global Biogeochemical Cycles waters to travel from the surface to 200-1,000 m depth interval (orange line in Figure 7a), particularly in the equatorial Pacific (Figures 8a and 8b). Concurrently, net primary production (blue dashed line in Figure 7b), export production (orange dashed line in Figure 7b) and by inference the respiratory O 2 demand in the upper ocean decrease during phase I (Figure 8g and 8h). Thus, both increased ventilation and the reduced O 2 consumption from remineralization contribute to the increase in O 2 circ/rem/diseq . The latter process might be even underestimated in our model, since compared to other CMIP5-type Earth system models, its simulated reduction in NPP and export production in response to transient global warming tend to be on the low end of the range (Laufkötter et al., 2015(Laufkötter et al., , 2016. After about 450 years, ventilation ages in the upper ocean of the (sub)tropics increase again and overshoot the preindustrial water mass age by Year 2500 (Figure 7a). Afterward, ventilation ages decrease slowly but remain higher than under preindustrial conditions. In particular, the upper waters of the tropical and midlatitude Pacific Ocean are bathed by older waters by the end of the simulation (Figure 8c). During phases II and III, export production in the (sub)tropics continues to decrease slightly (Figure 7b), leading to total reduction of −4% or −0.3 Pg C/year integrated by the end of phase III. Interestingly, this decrease in export production occurs while NPP in the sub (tropics) temporarily increases by 2% or 0.8 Pg C/year during phase II (blue dashed line in Figure 7b). This decoupling between export production and NPP could be a consequence of the strong temperature dependence of the phytoplankton growth rates responding to the higher temperatures and a decreasing export efficiency due to a shift in the phytoplankton production toward small phytoplankton (Laufkötter et al., 2016).
Next, we turn our discussion to drivers of O 2 circ/rem/diseq changes in the deep ocean.  (Figures 4g-4i). The ventilation age strongly decreases in the Southern Ocean and in the Ross Sea, as well as in the Arabian Sea and along the Chilean Margin during phase III (Figure 8f). The North Pacific becomes first more sluggishly ventilated during phase II (Figure 8e) leading to a decrease in O 2 circ/rem/diseq ( Figure 6), but then during phase III, it becomes better ventilated ( Figure 8f) and oxygenated (Figure 6i). In contrast, the deep North Atlantic remains less ventilated until the end of the simulation. The simulated overall decrease in export production (orange solid line in Figure 7b) and associated reduction in O 2 consumption may also contribute to the O 2 increase in the abyssal ocean. However, the strong correlation between changes in ventilation age and O 2 changes indicate that changes in export production are not the main drivers of deep ocean O 2 changes. An exception relates to the eastern equatorial Pacific, where O 2 decreases

Long-Term Recovery of the Southern Ocean Ventilation
Next, we consider the processes leading to the ventilation increase in the Southern Ocean at the onset of phase III. The AMOC decreases from 25.0 Sv (1 Sv ¼ 10 6 m 3 s −1 ) to 15.6 Sv during the first 80 years of the simulation but then recovers and overshoots its initial value after 1,000 years by nearly 10% (red line in Figure 9a). The AMOC decrease and "overshoot" is qualitatively in line with previous coupled climate model simulations (e.g., Battaglia & Joos, 2018;Krasting et al., 2018;Rugenstein et al., 2016;Schmittner et al., 2008;Zickfeld et al., 2013). The AMOC becomes shallower during phase I and remains relatively shallow even Similarly, the SMOC reaches its weakest state (−0.2 Sv) soon after atmospheric CO 2 has stabilized (blue line in Figure 9a). In contrast to the AMOC, the SMOC does not steadily converge to a new equilibrium thereafter. Within the first few decades after the minimum strength, SMOC recovers to its preindustrial strength but then weakens again until the end of phase II. However, the small weakening is superimposed by large interannual variability. At the beginning of phase III, the SMOC increases again and reaches a new equilibrium strength of −18.5 Sv by the end of phase III (38% stronger than at preindustrial). This strengthening leads to a stronger ventilation of the deep Southern Ocean and subsequently parts of the Indo-Pacific Ocean (Figure 8f) and, therefore, to an increase in O 2 circ/rem/diseq (Figure 6i) and O 2 levels ( Figure 6c).
The abrupt reinvigoration of the SMOC is somewhat unexpected, as the thermal and especially the haline surface density flux actually decrease during phases I and II (Figure 10a)-circumstances that are expected to lead to enhanced water-column stratification and a reduction of the formation of Antarctic Bottom Water (de Lavergne et al., 2014). Changes in atmospheric wind patterns over the Southern Ocean also cannot explain the sudden onset of deepwater formation in the Southern Ocean. During phase I, the model simulates a poleward shift of the Southern Hemisphere westerlies by about 1°(blue line in Figure 10b) and a small strengthening (red line in Figure 10b). Both changes are common features simulated in Earth system models under transient global warming (Russell et al., 2006). However, during phases II and III, the westerlies gradually weaken and move equatorward. Thus, they do not show any sudden changes in intensity or position that could explain the abrupt onset of deepwater formation.
The abrupt increase in SMOC during the transition from phases II to III is primarily driven by the onset of deepwater formation in the Ross Sea. The September mixed layer in the Ross Sea abruptly deepens from about 170 m to more than 600 m within a few years (orange line in Figure 10c) with a concomitant abrupt sea ice retreat (orange line in Figures 10d and 11). In contrast, the mixed layer in the Weddell Sea steadily deepens during phase II after an initial shoaling (green line in Figures 10c and 10d). Also, the sea ice extent in the Weddell Sea changes only gradually. Thus, the changes in the Weddell Sea cannot be responsible for the abrupt increase in SMOC. But the Weddell Sea is likely responsible for the partial SMOC recovery at the beginning of phase II and the partial recovery and the slowdown in O 2 circ/rem/diseq decrease toward the end of phase II (Figure 2d).
Which processes lead to the abrupt change in the mixed layer depth in the Ross Sea? During phases I and II, the uppermost 75 m of the Ross Sea waters freshen (decrease in salinity; Figure 12b) due to an increase in precipitation and in melting rates (green line in Figure 10a). This freshening together with a very small warming ( Figure 12) decreases surface density ( Figure 12d) and increases the vertical stratification of the upper water column as expressed by the squared Brunt-Väisälä frequency (Figure 12c). This increase in the stratification reduces the transport and mixing of subsurface waters, which usually carry saltier and warmer waters toward the surface, where they are cooled by the atmosphere and freshened by an excess of freshwater input from the atmosphere and the melting of sea ice (Haumann et al., 2020). Through the reduction in this upward transport and mixing, the increased stratification leads to a buildup of positive anomalies in temperature (up to 1.5°C; Figure 12a) and in salinity (up to 0.1 psu; Figure 12b) below~75 m. These anomalies tend to compensate each other in terms of their effect on subsurface density (Figures 12d-12f). However, as outlined by Martin et al. (2013), these subsurface changes prime the system toward a potential instability. This can happen once these anomalies breach the stratification and make it to the surface. There the warm waters would lose their signal quickly through cooling by the atmosphere, while the excess salinity would stay behind, as the freshening occurs more slowly. This leads then to especially dense surface waters which can rapidly initiate deep convection. In addition to this priming effect by the accumulation of heat and excess (a) (b) (c) (d) Figure 11. September sea ice concentration in the Southern Ocean for the (a) first 5 years of the simulation, (b) end of phase II (average over period 701 to 720), (c) start of phase III (average over period 721 to 740), and (d) end of phase III (average over period 4481-4500). The simulated 80% sea ice concentration contour is shown in red, and the September mixed layer depth contours of 2,000 and 4,000 m are shown as magenta dashed and solid lines, respectively. In (a), the green line shows the observation-based 80% sea ice concentration contour averaged over the 1986-2005 period (Rayner et al., 2003) and the gray solid lines depict the boundaries of the Ross Sea (RS) and Weddell Sea (WS) region, respectively.

10.1029/2020GB006601
Global Biogeochemical Cycles salinity at depth, sea ice is continually retreating during phases I and II (Figures 11a and 11b) permitting a rapid cooling through heat release into the atmosphere when these warm and salty water mass eventually reach the surface. Therefore, at the end of phase II, the Ross Sea is optimally primed by meeting several conditions that favor open ocean deep convection.
During the transition from phases II to III, the deep water formation initiates in the Ross Sea, which results in an abrupt change in density, temperature, and salinity structure within only a few years ( Figure 12) as well as to intensified sea ice melting (orange line in Figures 10d and 11c) and enhanced air-sea heat flux (orange line in Figure 10a). The onset of deepwater formation also leads to a rapid increase in atmospheric temperature south of 55°S (orange line in Figure 1b). The open ocean deep convection in the Ross Sea persists throughout phase III, with the exception of the 1150-1200 period. In the new equilibrium, the water column is well mixed in the upper 500 m (Figure 12c). The exact timing of the abrupt increase in deep water formation in the Ross Sea is stochastic in nature, in agreement with previous modeling studies (Martin et al., 2013). We demonstrate this by using an additional simulation that was run over the period 673-822 with identical forcing but started from slightly different initial conditions (see section 2.2 for description of the simulation). The onset of the deepwater formation in the Ross Sea is shifted by 22 years between the two simulations (not shown) pointing to the stochastic nature of the deepwater formation onset. However, certain preconditioning such as high subsurface salinity must be met so that deep water formation can occur.

Discussion and Conclusion
Our study suggests that a millennial-scale recovery and strengthening of the deep Southern Ocean ventilation leads to a global ocean reoxygenation under long-term global warming. Even though the substantial warming of the whole ocean (

10.1029/2020GB006601
Global Biogeochemical Cycles overcompensates this solubility-driven reduction. This long-term equilibrium response to global warming is opposite to that occurring during the transient warming on centennial timescales, when both circulation and solubility tend to reinforce each other. This difference between the transient and the long-term equilibrium response of O 2 to global warming is largely driven by the deep ocean, reflecting its dominating contribution to the O 2 inventory of the ocean. Oxygen in the upper ocean of the low latitudes and especially in the OMZ responds in a distinctly different manner. There, O 2 decreases on centennial timescales and remains depressed on millennial timescales. This is a consequence of the effect of warming never being compensated by ventilation/biology on all timescales.
Our result of a continued loss of O 2 in the (sub)tropical upper ocean and an associated expansion of OMZs under protracted global warming appears to be in contrast to the results from Fu et al. (2018). Using simulations with the NCAR CESM1 up to the Year 2300 under the RCP8.5 scenario and its extension, they show that after nearly two centuries of declining O 2 trends, the O 2 content of the tropical upper ocean undergoes a strong trend reversal around 2150, resulting in a near-full recovery of the O 2 concentration and the OMZ volume back to the preindustrial state. Unfortunately, the results from the two studies cannot be compared directly, since the two simulations employed two very different forcing scenarios, that is, 1% increase in atmospheric CO 2 for 70 years up to nearly 600 ppm, followed by a fixed forcing in our case, versus a slower increase over 400 years all the way up to about 2,000 ppm in the Fu et al. case. The achieved global ocean warming in the latter simulation was nearly 6°C, while it is less than 2°C in our simulation. In addition, the two studies compared two rather different time periods, that is, 200 years during which the forcing continued to increase, versus a multicentury to millennial response after the forcing had been stabilized in our case. Still, the two very different responses to ocean warming hold two important messages.
First, it points to the fundamental challenge associated with the modeling of the O 2 content of the tropical upper ocean (cf. Cabré et al., 2015;Cocco et al., 2013). Changes in the O 2 concentration in these regions are often the net result of strongly opposing tendencies stemming from warming/solubility and circulation/biology, with subtle changes in either of the two leading to fundamentally different net results. Similar to our study, Fu et al. (2018) also find that the export production decreases and that the ventilation slightly increases, overcompensating the solubility-driven loss of O 2 . In our simulation, however, the ventilation decreases toward the end of the simulation and the solubility-driven decrease in O 2 overcompensates any O 2 change due to the decrease in biological production and the changes in ventilation. The global export production decreases only by a few percent in the GFDL-ESM2M (Figure 7b), whereas the global export production declined by about 30% in the model simulations analyzed in the Fu et al. study .
The second message emerging from the comparison is that both studies emphasize that when it comes to O 2 changes in the tropical upper ocean under global warming, trends during one phase of a transient (e.g., 21st century) cannot be extrapolated to another phase of the transient (22nd century). And even further, trends during any phase of a transient cannot be used as an indication for what the final equilibrium response to warming will be. We conclude that there is high uncertainty in how the ocean and its biogeochemistry will change under long-term global warming. So far, the climate modeling community has strongly focused on centennial timescales and developed the model components needed for that relatively short timescale. However, to understand the millennial-scale equilibrium response of the Earth system to global warming, we need to focus on the processes that impact the Earth system on these long timescales, such as ice sheet-ocean interactions, and the representation of the deepwater formation and the biological pump.
Our simulation indicates a strong recovery of the deep ocean ventilation under a warming climate. Similar deep ocean circulation changes have been simulated in previous multicentennial to millennial global warming simulations with reduced-complexity models (Battaglia & Joos, 2018;Bi et al., 2001;Ridgwell & Schmidt, 2010;Rugenstein et al., 2016;Schmittner et al., 2008). Using different greenhouse gas scenarios, these studies show a decrease in the Antarctic Bottom Water formation with transient global warming but a partial or full recovery under a (nearly) equilibrated warmer climate. The only available study to date using a comprehensive ESM to investigate changes in deep ocean ventilation and oxygenation on millennial timescales confirms our findings (Yamamoto et al., 2015). They also show that the SMOC recovers after approximately 500 to 700 years followed by a subsequent oxygenation of the deep ocean. However, Yamamoto et al. (2015) suggest that the recovery of the deep ocean ventilation is entirely caused by the recovery of the Weddell Sea convection, in contrast to our simulation suggesting strong Ross Sea convection under 10.1029/2020GB006601

Global Biogeochemical Cycles
warming climate in addition to Weddell Sea convection. Yamamoto et al. (2015) also suggests that the O 2 decrease in the deep North Atlantic may be caused by a weaker AMOC. However, in our simulation, the AMOC recovers, but O 2 remains low. The reduced ventilation and associated O 2 decrease is rather caused by the shallowing of the AMOC cell, which leads to less ventilation of the deep North Atlantic even under a stronger AMOC. As neither Yamamoto's ocean model nor our model includes an interactive ice sheet and both models have a relatively coarse resolution, the results of both studies should be viewed with caution.
Our results of a better ventilated and oxygenated deep ocean but a less oxygenated upper ocean under a warming climate is also consistent with marine proxy records that indicate generally lower O 2 levels in the deep ocean yet higher O 2 concentrations in the upper ocean during the last glacial period (Galbraith & Jaccard, 2015;Jaccard et al., 2014;Jaccard & Galbraith, 2012). Interestingly, the deep ocean became better oxygenated as the Earth climate emerged from the last ice age as deep ocean ventilation increased in the Southern Ocean, possibly associated with sea-ice retreat and/or wind shifts (Jaccard et al., 2016;Skinner et al., 2010). In any case, additional long-term simulations with Earth system models are needed to systematically assess changes in the ocean beyond 2100.
Our simulation suggests that in both, the upper and deep ocean, the largest changes in O 2 occur many centuries or even millennia after atmospheric CO 2 has stabilized. In fact, about 80% of the deep ocean and subsurface ocean O 2 decrease occurs after atmospheric CO 2 has stabilized. Therefore, the largest impacts on marine life may occur long after the stabilization of the radiative forcing. In the geological past, reductions in intermediate-water O 2 had significant impacts on marine biodiversity and functioning (Moffitt et al., 2015;Myhre et al., 2017). In the upper ocean at low latitudes, where O 2 levels are already low, a decrease in O 2 may have a wide range of biological, ecological, and biogeochemical consequences, including shifts in body size, taxonomic composition, diversity, lifestyles, life stage, and trophic patterns of marine organisms (Breitburg et al., 2018;Cheung et al., 2013;Levin, 2018) as well as changes in ocean biogeochemistry, primarily the marine nitrogen cycle with consequences for marine productivity (Gruber, 2008). A large decrease in O 2 may therefore have deleterious consequences on marine life, in addition to stress arising from other adverse ocean ecosystem drivers such as heatwaves, acidification, and nutrient stress (Bopp et al., 2013;Frölicher et al., 2016Frölicher et al., , 2018Gruber, 2011).
Even though our simulated O 2 changes, in particular the deep ocean O 2 changes, compare qualitatively well with the past O 2 changes covering the last deglaciation and with recent simplified ESM simulations, a number of caveats need to be discussed. First and most importantly, the model lacks a dynamic ice sheet component, and meltwater input (and albedo changes) from Antarctic and Greenland ice melting is not included. Any additional freshwater might lead to a larger reduction in the SMOC and AMOC than simulated in our study and might strongly delay or prohibit the recovery of the SMOC and AMOC and consequently also deep ocean O 2 . However, the mechanism (surface freshening and subsurface warming) responsible for the SMOC recovery seen in our simulation may also be at play under an even stronger stratification with stronger surface freshening but also stronger subsurface warming.
Second, the GFDL-ESM2M predominantly forms deep water in the Southern Ocean through open ocean convection in the Weddell Sea, illustrated by the overly deep mixed layer and low bias in sea ice extent in this region (Figure 11a). The combination of sea ice and iceberg export of freshwater and winter heat loss leads to open ocean convection in the Weddell Sea in the GFDL-ESM2M. In the "real" Southern Ocean, deep waters are predominantly formed on the continental shelves of the Weddell and Ross Sea through brine rejection and interactions with grounded ice sheets (Drucker et al., 2011;Orsi, 2010) and only rarely in open ocean polynyas (Killworth, 1983). Therefore, the resupply of O 2 to the deep ocean through Southern Ocean deep water formation in the Ross and Weddell Sea in the GFDL-ESM2M may not be realistic. However, no Earth system model currently represents the complex physical processes in the Southern Ocean correctly (Heuzé, 2020), because in part the resolution is too coarse, and in part the models are missing some crucial processes, such as cooling underneath the glacial ice.
Third, current ESMs including GFDL-ESM2M show large discrepancies in simulating present-day and future trends in OMZs and generally overestimate the observed OMZ volume (Cabré et al., 2015;Cocco et al., 2013;Fu et al., 2018). Higher resolution ESMs will hopefully improve the representation of OMZs 10.1029/2020GB006601

Global Biogeochemical Cycles
FRÖLICHER ET AL. and deep water formation in the Southern Ocean and offer an improved way for investigating the response of OMZs and the Southern Ocean circulation to global warming (e.g., Busecke et al., 2019;Morrison et al., 2016;Newsom et al., 2016). However, the length of such high-resolution model simulations will be restricted to a few hundred years due to the enormous computational costs. To extend such high-resolution simulations to the millennial timescale will require a significant increase in community-level computational resources.
Fourth, the magnitude and timescale of response of the meridional overturning circulation is sensitive to the selected forcing scenario (Battaglia & Joos, 2018;Rugenstein et al., 2016;Stouffer & Manabe, 2003) and model formulation, and the SMOC as well as the AMOC may not recover under stronger greenhouse gas forcing scenarios. This limitation could be overcome by analyzing a suite of different multimillennial global warming simulations. The new LongRunMIP offers such a unique opportunity (Rugenstein et al., 2019).
And fifth, the TOPAZv2 lacks certain processes that may alter future O 2 changes in addition to the processes that we describe above. For example, recent studies have indicated that the carbon-to-nitrogen utilization ratio of organic matter production and therefore O 2 respiration may increase under higher ocean carbon conditions (Riebesell et al., 2007). While these perturbations to the nutrient cycling may alter the ocean's O 2 content , these changes would mainly be restricted to the upper ocean. In addition, the TOPAZv2 module does not include a temperature dependent remineralization (Bendtsen et al., 2015;Laufkötter et al., 2017) and particulate sinking speed, and therefore, a possible decreasing or increasing transfer efficiency of organic matter to the deep sea under global warming cannot be ruled out as an additional driver of upper and deep O 2 changes in the real world. Also, the potential changes in atmospheric nutrient deposition and sediment interactions are not considered here.
In summary, our study suggests that the equilibrium response in oceanic O 2 is not simply a scaled version of the transient response, especially in the deep ocean, and that future changes in O 2 may occur centuries after atmospheric CO 2 has stabilized. Given the potential for substantial long-term multimillennial-scale changes in ocean O 2 under a doubling of atmospheric CO 2 , this study provides further evidence that future global warming should be kept at a minimum to prevent substantial and sustained negative impacts on marine life.

Data Availability Statement
Code and data can be found under: https://doi.org/10.5281/zenodo.3982620. respectively). T. L. F. also acknowledges support from the Swiss National Supercomputing Centre (CSCS). This publication has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement 820989 (project COMFORT, Our common future ocean in the Earth system-quantifying coupled cycles of carbon, oxygen, and nutrients for determining and achieving safe operating spaces with respect to tipping points).