Hazards of decreasing marine oxygen: the near-term and millennial-scale beneﬁts of meeting the Paris climate targets

. Ocean deoxygenation is recognized as key ecosystem stressor of the future ocean and associated climate-related ocean risks are relevant for current policy decisions. In particular, beneﬁts of reaching the ambitious 1.5 ◦ C warming target mentioned by the Paris Agreement compared to higher temperature targets are of high interest. Here, we model oceanic oxygen, warming and their compound hazard in terms of metabolic conditions on multi-millennial timescales for a range of equilibrium temperature targets. Scenarios where radiative forcing is stabilized by 2300 are used in ensemble simulations with the Bern3D Earth System Model of Intermediate Complexity. Transiently, the global mean ocean oxygen concentration decreases by a few percent under low forcing and by 40 % under high forcing. Deoxygenation peaks about a thousand years after stabilization of radiative forcing and new steady-state conditions are established after AD 8000 in our model. Hypoxic waters expand over the next millennium and recovery is slow and remains incomplete under high forcing. Largest transient decreases in oxygen are projected for the deep sea. Distinct and near-linear relationships between the equilibrium temperature response and marine O 2 loss emerge. These point to the effectiveness of the Paris climate target in reducing marine hazards and risks. Mitigation measures are projected to reduce peak decreases in oceanic oxygen inventory by 4.4 % ◦ C − 1 of avoided equilibrium warming. In the upper ocean, the decline of a metabolic index, quantiﬁed by the ratio of O 2 supply to an organism’s O 2 demand, is reduced by 6.2 % ◦ C − 1 of avoided equilibrium warming. Deﬁnitions of peak hypoxia demonstrate strong sensitivity to additional warming. Volumes of water with less than 50 mmol O 2 m − 3 , for instance, increase between 36 % and 76 % ◦ C − 1 of equilibrium temperature response. Our results show that millennial-scale responses should be considered in assessments of ocean deoxygenation and associated climate-related ocean risks. Peak hazards occur long after stabilization of radiative forcing and new steady-state conditions establish after AD 8000.


Introduction
Oxygen (O 2 ) is a sparingly soluble gas and its abundance in the ocean is decreasing under ongoing global warming (IPCC, 2013). Decreasing O 2 concentrations, warming and changes in other environmental parameters forced by anthropogenic greenhouse gas (GHG) emissions pose high risks for marine ecosystems (Gattuso et al., 2015). The parties to the United Nations Framework Convention on Climate Change (UNFCCC, 2018) and to the Paris Agreement note "the importance of ensuring the integrity of all ecosystems, including Oceans" (Paris Agreement, 2018) and "the threats of irreversible damage" (UNFCCC, 2018, Article 3). Marine changes are projected to evolve over multi-century and millennial timescales with peak impacts occurring potentially long after stabilization of atmospheric GHG concentrations and peak temperatures. However, only a few studies have assessed millennial-scale impacts of anthropogenic GHG emissions on the ocean and the reversibility of marine changes in oxygen. Explicit quantification of the benefits of meet-Published by Copernicus Publications on behalf of the European Geosciences Union. source: https://doi.org/10.7892/boris.118488 | downloaded: 5.3.2020 G. Battaglia and F. Joos: O 2 projections for multi-millennial timescales ing the 2 or 1.5 • C climate targets mentioned by the Paris Agreement with respect to the reversibility and avoidance of implied impacts on marine oxygen and related environmental parameters, including ocean circulation, ocean warming, metabolic viability and biological productivity on multimillennial timescales is missing.
Typical thresholds leading to O 2 stress for many macroorganisms (hypoxia) are around 50 mmol O 2 m −3 . Water with lower O 2 concentrations are effectively dead zones for many higher animals (reviewed in Keeling et al., 2010;Storch et al., 2014;Breitburg et al., 2018). Species are also sensitive to thermal stress (Gattuso et al., 2015) and their sensitivity to hypoxia increases with higher temperatures (Pörtner, 2010). In the modern ocean, oxygen-poor zones with O 2 < 50 mmol m −3 occupy about 5 % of its volume (Garcia et al., 2014;Bianchi et al., 2012). The expansion of oxygenpoor waters leads to habitat compression, mortality and major changes in community structure where energy preferentially flows into microbial pathways to the detriment of higher trophic levels. Suboxic (< 5 mmol O 2 m −3 ) or anaerobic conditions can also lead to production of poisonous H 2 S within sediments (reviewed in Diaz and Rosenberg, 2008;Breitburg et al., 2018), and decreasing O 2 concentrations potentially lead to higher production and emissions of the greenhouse gas nitrous oxide.
Observational Ito et al., 2017) and modeling studies  indicate an overall decline in the oceanic oxygen content over past decades. Systematic discrepancies exist for the typically low-oxygen tropical thermocline, where observations suggest O 2 has decreased and most models simulate increased O 2 levels over the past decades. Model projections to the end of the 21st century consistently project the global ocean oxygen inventory to further decline with anthropogenic climate change (Matear, 2000;Plattner et al., 2001;Bopp et al., 2002Bopp et al., , 2013Frölicher et al., 2009;Cocco et al., 2013). The most recent generation of Earth system models simulate global deoxygenation by the end of the 21st century of around −1.81 % (RCP2.6) to −3.45 % (RCP8.5; IPCC, 2013). Impact studies have highlighted potential habitat compression (Deutsch et al., 2015;Mislan et al., 2017) and reduced catch potential (Cheung et al., 2016) associated with climate change at the end of the century. Large model-model differences remain in projections of oxygen minimum zones (OMZs; Cocco et al., 2013;Bopp et al., 2013).
Given the long residence time of anthropogenic CO 2 in the atmosphere and long equilibration timescales of the ocean overturning circulation, anthropogenic climate change will grow and persist beyond the end of the 21st century, the typical near-term assessment timescale of climate change (Clark et al., 2016). Only a few studies have assessed ocean biogeochemistry and the oceanic oxygen content beyond this nearterm timescale. Available studies employ a range of physical and biogeochemical complexity levels from box models to general circulation models (GCMs). Oxygen concentra-tions are simulated to decline beyond the 21st century on multi-centennial timescales (Matear and Hirst, 2003;Hofmann and Schellnhuber, 2009;Mathesius et al., 2015). Simulations covering 2 millennia show a recovery phase thereafter (Schmittner et al., 2008;Yamamoto et al., 2015). In most studies, simulated oxygen concentrations have not reached new steady-state conditions at the end of the simulation. Low-order Earth system models and those of intermediate complexity (EMIC) integrated by up to 100 000 years have demonstrated the potential for long-term ocean oxygen depletion in response to carbon dioxide emissions and the long equilibration timescales of ocean biogeochemical variables in response to carbon emissions (Shaffer et al., 2009;Ridgwell and Schmidt, 2010). Multi-millennial simulations are therefore required to assess the full amplitude of ocean biogeochemical changes and new steady-state conditions due to anthropogenic climate change.
The distribution of O 2 in the ocean results from the sum of its solubility component set through air-sea exchange, the effect of O 2 production by phytoplankton in the euphotic zone and O 2 consumption during organic matter remineralization at depth. In modeling studies, it is possible to identify the drivers of O 2 changes by considering changes due to solubility and changes due to oxygen consumption. When assessing the near-term timescale at the end of 21st century, studies have shown that different depths in the water column tend to be associated with different dominant mechanisms of change. In the surface ocean, O 2 changes tend to be determined by changes in O 2 solubility. In the subsurface, both changes in solubility and utilization may reinforce (mid-and high latitudes) or compensate for each other (tropics; e.g., Cabre et al., 2015;Bopp et al., 2017). In the deep ocean, simulated O 2 changes are dominated by changes in O 2 utilization, which is in turn controlled by ocean ventilation (see also Matear and Hirst, 2003;Yamamoto et al., 2015, for longer timescales). Changes in the oceanic heat content and in ocean circulation are therefore crucial for O 2 changes.
Well-defined metrics that summarize the Earth system response are useful in many aspects and may facilitate communication in the mitigation policy context of the Paris Agreement. The Transient Climate Response to Cumulative Carbon Emissions (TCRE, Allen et al., 2009) or the Transient Earth System Response to Cumulative Carbon Emissions (TREX, Steinacher and Joos, 2016) are such metrics. These link changes in global surface air temperature and environmental parameters to cumulative CO 2 emissions relying on near-linear relationships. Similar metrics may be developed for oceanic oxygen.
Deoxygenation is one of several marine ecosystem stressors including warming, acidification, hypocapnia, changes in food supply and sea-level rise (Gruber, 2011;Cocco et al., 2013;Bopp et al., 2013;Gattuso et al., 2015;Sweetman et al., 2017). Several key marine and coastal ecosystems may face severe risks due to climate change even if the suggestions for low emission pathways are followed through the Earth Syst. Dynam., 9, 797-816, 2018 www.earth-syst-dynam.net/9/797/2018/ end of the century (Gattuso et al., 2015;Magnan et al., 2016;Breitburg et al., 2018). This growing body of concern helped drive the 21st Conference of the Parties (COP21) to reach the Paris Agreement. A goal of the Paris Agreement is to "hold the increase in the global average temperature to well below 2 • C above pre-industrial levels and pursuing efforts to limit the temperature increase to 1.5 • C above pre-industrial levels" (Article 2a, Paris Agreement, 2018). Article 4 of the Paris Agreement further states that "[in] order to achieve the long-term temperature goal set out in Article 2, Parties aim to reach global peaking of greenhouse gas emissions as soon as possible . . . and to undertake rapid reductions thereafter in accordance with best available science, so as to achieve a balance between anthropogenic emissions by sources and removals by sinks of greenhouse gases in the second half of this century" (Article 4, Paris Agreement, 2018).
In this study, we assess the effectiveness of the Paris climate targets in reducing hazards of decreasing oceanic oxygen, ocean warming and marine export productivity as simulated by the Bern3D Earth System Model of Intermediate Complexity. We prescribe in the model four different, idealized scenarios where anthropogenic GHG forcing is stabilized by AD 2300. The four scenarios are designed to reach an equilibrium warming of 1.5, 1.9, 3.3 and 9.2 • C above preindustrial. Simulations are run to year AD 10 000 by which time the ocean has reached new steady-state conditions. This allows us to assess reversibility and the full amplitude of changes, acknowledging the long equilibration timescale of biogeochemical variables with peak hazards potentially occurring long after stabilization of radiative forcing in the atmosphere. We summarize the outcomes developing global metrics which quantify avoided marine hazards per avoided global warming on three different time horizons. The first time horizon is the end of the 21st century, the typical assessment timescale of climate change hazards. Here, changes in a variable are related to changes in surface air temperature (SAT) in the year 2100. Those are contrasted to the millennial-scale perspective where peak changes in the variable in the course of the simulation and equilibrium changes at the end of the simulation are related to the corresponding equilibrium warming.
In Sect. 2, we briefly describe the Bern3D model and the experimental setup. Four different radiative forcing stabilization scenarios to meet four temperature targets (1.5, 1.9, 3.3 and 9.2 • C above preindustrial) are considered. The observation-constrained 100-member ensembles used to explore parameter uncertainties for each scenario is introduced. In Sect. 3, physical changes, including changes in overturning, water mass age, sea ice, temperature, salinity and density as well as biogeochemical changes, including changes in global oxygen inventory, the extent of oxygen minimum zones and productivity, are presented. The compound effects of warming and oxygen changes are assessed in the form of a metabolic index (Deutsch et al., 2015). Underlying physical and biogeochemical processes and mechanisms are dis-cussed. Following earlier studies, we attribute the contributions of O 2 changes from changes in solubility, and the interplay of ocean biology and ventilation by carrying four explicit O 2 tracers and an ideal age tracer. The graphical illustration of spatial changes is focused on the 1.5 • C equilibrium warming target at the point of peak O 2 decline. Additional supporting figures are given in the appendix. In Sect. 4, the relationship between change in global mean surface air temperature ( SAT) and selected impact-relevant parameters is quantified. Relationships are established for the near-term (AD 2100), the time of the peak decline in oxygen around AD 3000 to 4000 and at year AD 10 000 when a new equilibrium has been reached in the model. Often a near-linear relationship is found between the change in a variable of interest and the change in SAT as simulated across the range of scenarios and ensemble members at a distinct time. This allows us to develop new metrics to quantify avoided marine hazards per unit change in SAT at different points in time. These quantitatively illustrate the benefits of meeting the Paris target in terms of marine hazards. Each modeling exercise is associated with uncertainties, and in Sect. 5, we discuss relevant uncertainties, mention neglected processes and compare our findings to other studies. Finally, in Sect. 6 we present implications and conclusions, and we summarize our findings graphically for a "1.5 • C world" and contrast peak changes across the range of temperature targets.

Bern3D
The Bern3D Earth System Model of Intermediate Complexity is a three-dimensional frictional geostrophic balance ocean model (Müller et al., 2006), which includes a sea-ice component coupled to a single-layer energy and moisture balance model of the atmosphere (Ritz et al., 2011) and a prognostic marine biogeochemistry module (Tschumi et al., 2011;Parekh et al., 2008). A version with 41 × 40 horizontal grid cells and 32 vertical layers is used (see also Roth et al., 2014;Battaglia et al., 2016 for model evaluation). The NCEP/NCAR monthly wind stress climatology (Kalnay et al., 1996) is prescribed at the surface. Air-sea gas exchange, carbonate chemistry and natural 14 C of dissolved inorganic carbon is modeled according to OCMIP-2 protocols Orr and Najjar, 1999;Orr and Epitalon, 2015). The global mean air-sea transfer rate is reduced by 19 % compared to OCMIP-2 to match observation-based estimates of natural and bomb-produced radiocarbon .
The biogeochemical module is based on phosphorus and simulates production and remineralization/dissolution of organic matter, calcium carbonate and opal. Production of particulate organic phosphorus (POP) within the euphotic zone (top 75 m) depends on temperature, light availability, phosphate and iron following Doney et al. (2006). POP reminer-alization within the water column follows a power law profile (Martin et al., 1987). Organic matter falling on to the sea floor is remineralized in the deepest box. Two-thirds of organic phosphorus production form dissolved organic phosphorus (DOP), which decays with an e-folding lifetime of 1.5 years. An updated remineralization scheme assigns remineralization of POP and DOP to aerobic and anaerobic pathways depending on the mean grid cell dissolved O 2 concentration (see Battaglia and Joos, 2018). We introduce two power law profiles with two distinct remineralization length scales for aerobic and anaerobic remineralization (α aerob and α denit ). Constant stoichiometric ratios are used for both aerobic and anaerobic remineralization to convert biological P fluxes into carbon and alkalinity fluxes (P : Alk : C = 1 : 17 : 117). The O 2 demand for complete aerobic remineralization is 170 mol O 2 mol PO 4 and no oxygen is consumed for anaerobic remineralization. Accordingly, aerobic remineralization in the ocean is smaller than O 2 production in the euphotic zone leading to an O 2 outgassing for steadystate conditions. The atmospheric oxygen inventory is constant. This is justified as 99.5 % of the ocean-atmosphere inventory is in the atmosphere and potential net fluxes of O 2 from the ocean and land to the atmosphere and fossil fuel burning have a small impact on atmospheric O 2 . O 2 components from O 2 production, consumption and solubility are carried as explicit model tracers to attribute changes. Tracers add up to within 10 −14 Pmol with mean inventories of 23.2, −239, 430 yielding a total of 214 Pmol, respectively (median values given). O 2 components inferred from O 2 saturation can result in systematic errors from surface disequilibrium (Ito et al., 2004). The use of explicit tracers avoids such systematic errors in the O 2 components. As changes in the O 2 production term are small, we combine the O 2 production and consumption tracers to a O 2 biology tracer when displaying sections.
We include evaluation of a metabolic index, , which was proposed by Deutsch et al. (2015). It combines temperature and pO 2 as indicators of metabolically viable environments and is defined as the ratio of O 2 supply to an organism's resting O 2 demand. We consider only relative changes in (t) relative to a reference time, t 0 (average over 1870-1899): where T is the absolute temperature, k B is Boltzmann's constant and the exponential function and the parameter E 0 characterize the temperature dependence of the baseline metabolic rate. E 0 only weakly affects the relative influence of temperature and O 2 gradients, and relative changes in are therefore independent of species (Deutsch et al., 2015).
Here, we consider E 0 = 0.87 eV (representative of Atlantic cod). For the calculation of pO 2 we pressure-correct the equilibrium constant following Eq. (5) in Weiss (1974). The metabolic index , as proposed by Deutsch et al. (2015), is linear in pO 2 (representing the rate of O 2 supply) and de-creases non-linearly with temperature (indicative of the resting metabolic demand). One may note that the exponential curve varies approximately linearly for typical global warming associated temperature changes as E 0 /k b (≈ 10 000 K) is large. The current set up does not include sediment interactions, temperature-dependent remineralization, variable stoichiometry, nitrogen cycle feedbacks, atmospheric nutrient deposition, dynamic wind or freshwater input/albedo changes from melting of continental ice sheets.

Ensemble and scenarios
To explore potential oxygen changes we set up four 100member ensembles each targeting a different equilibrium temperature response (∼ 1.5, 1.9, 3.3 and 9.2 • C above preindustrial). A feedback parameter λ (W m −2 K −1 ; Ritz et al., 2011), accounting for climate feedbacks that are not explicitly treated in the Bern3D model, is chosen in combination with radiative forcing from the Representative Concentration Pathways (RCPs; Meinshausen et al., 2011) to achieve these stabilization targets. RCP2.6, stabilizing by 2300, is run with λ values of −0.71 and −1 W m −2 K −1 , reaching the 1.5 and 1.9 • C targets, respectively. RCP4.5, stabilizing after 2100, is run with −1 W m −2 K −1 yielding a 3.3 • C temperature response and RCP8.5, stabilizing in the 23rd century, with −0.71 W m −2 K −1 yielding a 9.2 • C response (median values given for temperature targets). Each member is spun up over 5000 years to AD 1765 boundary conditions. The radiative forcing follows RCP scenarios (RCP2.6, 4.5 and 8.5; Meinshausen et al., 2011). The RCP scenarios are extended to year AD 10 000 by which time the ocean has equilibrated to new steady-state conditions. Radiative forcing includes an 11-year solar cycle up to year AD 3000. After that, all forcings are kept constant. We employ a single-model setup, and assess uncertainties arising from organic matter remineralization (α aerob and α denit ) and vertical mixing (k diff−dia ). The three parameters are sampled using the Latin hypercube sampling technique (McKay et al., 1979). The parameter ranges are chosen such that all members achieve similar skill scores with respect to observation-derived fields of natural radiocarbon (Key et al., 2004) and dissolved O 2 (Garcia et al., 2014;Bianchi et al., 2012) and correspond to the values chosen in Battaglia and Joos (2018 , Table 1). A normal distribution is used to sample α aerob with a standard value of −0.83 and a standard deviation of −0.0625. α denit is sampled uniformly between −0.1 and −0.01. And a lognormal distribution is used to sample k diff−dia (standard value = 2.25 × 10 −5 m 2 s −1 , shape parameter = 0.2, location parameter = 0). We choose a single ensemble member with parameter values close to the standard values as representative ensemble member to illustrate spatial anomalies (α aerob = −0.85, α denit = −0.037, k diff−dia = 2.05 × 10 −5 m 2 s −1 ).

Pre-industrial characteristics
The ensemble produces a range in overturning strengths, remineralization fluxes and O 2 distributions. The following numbers represent the 90 % confidence ranges of important model characteristics across the ensemble. The maximum of the Atlantic meridional overturning circulation (AMOC) stream function below 400 m depth ranges from 16.5 to 19.7 Sv. The minimum of the Indo-Pacific meridional overturning stream function below 400 m depth (Indo-Pacific MOC) ranges between −13.6 and −15.6 Sv. Export of particulate organic matter at 75 m ranges from 9.0 to 11.4 Gt C yr −1 . The simulated oxygen inventory ranges between 195 and 230 Pmol given the three parameters and the simulated oxygen distribution covers the observational range and spatial pattern well (see Figs. 3 and 7a and Table D1 of Battaglia and Joos, 2018). Biases exist in the simulated extent of OMZs. The volume of suboxic conditions (O 2 < 5 mmol m −3 ) is overestimated by a factor of 5 but water column denitrification fluxes are well within current estimates (Table D.1., Fig. 2c in Battaglia and Joos, 2018). This is a common model bias in EMICs and GCMs (Cocco et al., 2013;Bopp et al., 2013;Cabre et al., 2015). Vastly enhanced spatial resolution may be required to simulate equa-torial physics and ecosystems in better agreement with observations .

Marine changes in temperature, circulation and biogeochemistry
We first describe the evolution of important physical quantities that impact O 2 concentrations. Figure 1 displays the temporal changes in global mean surface air and ocean temperature, the evolution of annual mean sea-ice area in the Northern and Southern hemispheres and the Atlantic and Indo-Pacific meridional overturning circulations. In response to the RCP scenarios, atmospheric temperatures rise and stabilize after ∼ 1000 years (Fig. 1a). The four ensembles reaching 1.5, 1.9, 3.3 and 9.2 • C above preindustrial surface air temperature show an equilibrium ocean warming of 1.1, 1.3, 2.0 and 5.5 • C, respectively (median values given). Sea ice retreats in both hemispheres ( Fig. 1c and d). The retreat is more pronounced for higher forcing. In the Southern Hemisphere, even the lower forcing levels show strong decline in the annual mean sea-ice area and sea ice vanishes for higher forcing. The warming perturbation causes the AMOC and Indo-Pacific MOC to decline transiently (Figs. 1e, f and A1). The larger the forcing and im- plied changes in stratification, the larger the peak decline in overturning ( Fig. 1e and f). The decline is likely driven by upper ocean warming, leading to increasing surface-to-deep density gradients as further modulated by salinity changes (Fig. A2). The deep ocean water mass age increases in response to the slowed overturning (Figs. 2d and 3d). As retreating sea-ice increases wind stress over these newly exposed areas, younger water masses form in the upper ocean of the Southern Ocean (Fig. 3d). Reduced convection may contribute to a younger upper ocean through decreased entrainment of old deep water (not quantified within the scope of this paper). As the model tends to equilibrate under the sustained radiative forcing, the surface-to-deep gradients in the density anomalies diminish (Fig. A2), the meridional overturning circulation recovers ( Fig. 1e and f) and anomalies in water mass age become again smaller (Fig. 3d versus Fig. A3). The final circulation state is close but not identical to the preindustrial steady-state circulation. Maximum overturning strength in AMOC and the Indo-Pacific MOC varies by less than ±1 Sv around the initial value. At the new steady state, the maximum in AMOC below 400 m tends to be lower under higher forcing, whereas the maximum in the Indo-Pacific MOC below 400 m tends to be higher under higher forcing ( Fig. 1e and f). It is difficult and beyond the scope of this paper to conclusively explain such subtle changes in ocean dynamics and overturning (Fig. A1), likely linked to the complex changes in density (Fig. A2) and sea-ice retreat ( Fig. 1c and d). However, these differences have direct consequences for the projected global water mass age and by that for oceanic oxygen (Fig. 2) at the new equilibrium as further discussed below.
The response in oceanic oxygen is complex and characterized by an initial decline followed by a recovery phase (Fig. 2a). In line with earlier studies (Matear and Hirst, 2003;Schmittner et al., 2008;Shaffer et al., 2009;Ridgwell and Schmidt, 2010;Mathesius et al., 2015), our results demonstrate the potential for large changes in marine oxygen under anthropogenic forcing, a large inertia in the response and a slow and partially incomplete recovery of the perturbation. Transiently, the whole ocean oxygen inventory decreases by a few percent (6 %) under low forcing and by as much as 40 % under high forcing (median values given). The minimum in oxygen occurs about a thousand years after stabilization of radiative forcing, and it takes several millennia to approach a new equilibrium. Then, the global ocean O 2 inventory is a few percent higher than at preindustrial conditions under low and intermediate forcing and remains depleted by around 8 % in the high forcing case. global oxygen inventory (Fig. 2a) strongly correlate with water mass age (Fig. 2d) and are also impacted by gradual oxygen loss due to warming as evidenced by the evolution of the O 2 solubility tracer (Fig. 2b). Inventory changes based on the O 2 production tracer (Fig. 2c) are negligible; changes equilibrate with the atmosphere and only a small fraction remains in the ocean. The O 2 consumption tracer (Fig. 2e) determines the shape of the global O 2 signal (Fig. 2a). Its decline and recovery phase is strongly correlated with the evolution of ideal age ( Fig. 2e and d, see also Fig. 3b and d). As has been shown previously (e.g., Yamamoto et al., 2015;Bopp et al., 2017), the high correlation between changes in O 2 and ideal age and the absence of a direct relationship between changes in remineralization fluxes and O 2 , indicate that circulation changes are the major contributors to changes in O 2 . Changes from remineralization fluxes include both changes in absolute aerobic remineralization fluxes and changes in the relative share of denitrification (Fig. 2i). An increased share of denitrification at organic matter remineralization, for instance, effectively constitutes an implicit O 2 gain. Denitrification fluxes correlate with the volumetric expansion of OMZs and are also impacted by changes in remineralization fluxes within them (Fig. 2i). The recovery level of the O 2 consumption tracer (Fig. 2e) reflects the global recovery level of ideal age (Fig. 2d), where younger water masses are associated with less O 2 consumption and therefore higher O 2 concentrations. The total O 2 recovery level (Fig. 2a), on the other hand, is diminished due to O 2 loss from solubility (Fig. 2b). As such, 1.5 to 3.3 • C warming targets reach similar global O 2 equilibrium levels for different reasons. The 1.9 and 3.3 • C warming targets tend to result in younger water masses, which would increase O 2 due to less O 2 consumption compared to 1.5 • C warming targets. As those scenarios are also associated with higher warming, they lose more O 2 due to less solubility compared to 1.5 • C warming targets and yield similar global anomalies despite more pronounced spatial patterns. The 9.2 • C warming target reaches a lower equilibrium O 2 inventory compared to preindustrial due to high O 2 loss from solubility (−44.1 Pmol). We illustrate spatial changes in critical variables for a single representative ensemble member (see Sect. 2.2) at its peak O 2 decline, which occurs at year AD 3150 and amounts to 5 % (Fig. 3). The member eventually reaches the 1.5 • C warming target. Figure 3 displays anomalies in total O 2 (Fig. 3a), and the contributions from biologically mediated changes (termed "biology" below, Fig. 3b) combining the changes in the O 2 production and consumption tracer and from changes in solubility (Fig. 3c). In the upper ocean O 2 concentrations tend to increase due to biology and decrease due to solubility. Such compensating mechanisms have been documented elsewhere (e.g., Cabre et al., 2015;Bopp et al., 2017). The resulting changes in O 2 are less pronounced than  the changes in each component. The increase in O 2 due to biology stems from younger water masses and less export in the low and mid-latitudes (see next paragraph and Fig. 4c). O 2 changes show strong spatial correlation with changes in water mass age ( Fig. 3a and d). Largest decreases in O 2 are simulated in bottom waters in line with older water mass age. The equilibrium response in O 2 for this 1.5 • C warming case is characterized by slight O 2 decreases in the Atlantic, caused mainly by less solubility, and increases in the Southern Ocean and deep Pacific, caused by higher overturning and less sea-ice coverage in the Southern Ocean compared to preindustrial (Fig. A3). Global export production is simulated to decline over the first few centuries, and reach higher values under new steadystate conditions (Fig. 2g). The decline is stronger for higher forcing, while the recovery level of global export production is similar across the scenarios. Bern3D transiently simulates decreased export in the mid-and low latitudes ( Fig. 4c; see also Steinacher et al., 2009;Battaglia and Joos, 2018) as a result of increased stratification (Fig. A2c, f and i) and reduced nutrient concentrations in the surface ocean (Fig. 4b).
In the high latitudes, the model simulates increased export production, as a result of less temperature and light limitation as surface waters warm and sea ice retreats. This pattern of decreased export in mid-and low latitudes and increased export in high latitudes is similar across the scenarios. Export production in the low latitudes fully recovers for lower forcing and partially recovers for higher forcing. The lower recovery level in the low latitudes is compensated by higher increases in the high latitudes for high forcing. The magni-tude of positive and negative changes increases with forcing, but the global anomalies remain comparable at the end of the simulation.
Next to changes in export, we consider the evolution of a metabolic index in the upper ocean which integrates effects of changes in O 2 and temperature at the organism level (Figs. 2f and 4e). The globally averaged, upper ocean (depth < 400 m) metabolic index declines throughout the simulation dominated by increased temperatures (Fig. 2f). The metabolic index, (Deutsch et al., 2015), decreases in most places in line with warming and lower pO 2 (Fig. 4a, d and e). The O 2 gain in upper ocean waters is able to counteract the adverse effect of warming in some high-latitude environments. In other places with higher pO 2 , the temperature increase dominates the response in . Near-bottom waters in the Pacific are prone to largest reductions in , driven by large decreases in pO 2 (Fig. 4e).
Oxygen-poor waters (O 2 < 50 mmol m −3 , Fig. 2h) are simulated to transiently increase across all scenarios. The response is characterized by high uncertainty as introduced by the sampled parameters. Under new equilibrium conditions, the volume of low O 2 waters is reduced for low and intermediate forcing and remains higher than pre-industrial levels in the high forcing case.
Turning to uncertainties in our perturbed parameter ensemble, we find that variations in the vertical diffusion parameter (k diff−dia ) dominate the uncertainty in the globally averaged evolution of ideal age, sea-ice cover, temperature ter. Whether a threshold in O 2 concentration is met depends on the pre-industrial tracer distribution. Longer remineralization length scales bring more remineralization to depth, leading to higher O 2 consumption.

Metrics linking global warming to marine hazards
The purpose of this section is to quantify the relationship between changes in global mean surface air temperature (SAT), the target variable of the Paris Agreement, with selected aggregated metrics for marine ecosystem stressors. In this way, we link marine hazards to the temperature target of the Paris Agreement and quantify avoided marine change per unit of avoided global warming. Specifically, we investigate the relationship of SAT with changes in the marine O 2 inventory, ocean temperature and the metabolic index of Deutsch et al. (2015), the volume occupied by hypoxic water and in low-latitude export production (30 • S-30 • N) across the range of warming scenarios in our ensemble (Fig. 5). Distinct and often near-linear relationships emerge. Near-linearity allows us to characterize the benefits of avoided warming by single sensitivities, corresponding to the slopes of the relationships displayed in Fig. 5. The relationships between SAT and marine hazard metrics critically depend on the time horizon considered (Fig. 5). We compare and contrast changes at the end of the 21st century, the typical assessment timescale of climate change, to peak and equilibrium changes at the millennial timescale. Peak and equilibrium changes are related to the corresponding equilibrium temperature response, while changes at the end of the 21st century are related to the transiently realized warming at the end of the 21st century. Larger magnitudes are simulated on millennial timescales compared to the nearterm end of the 21st century. Assessment of ocean deoxy-genation by the end of the 21st century, therefore, underestimates the full amplitude of change.
Transient (end of 21st century), peak (AD ∼ 3000) and equilibrium (AD ∼ 8000) oxygen changes exhibit distinct relationships with their corresponding warming (Fig. 5a). At the end of the 21st century, simulated oxygen decreases by 0.68 % • C −1 of realized warming (median values). At peak oxygen decline, this sensitivity increases and oxygen decreases by 4.4 % • C −1 of equilibrium temperature response. In other words, an avoided warming of 1 • C, avoids a peak decline in marine O 2 inventory of 4.4 %. The linear relationship breaks down for the equilibrium response. While 1.5 to 3.3 • C warming targets lead to similar and higher oxygen levels, the 9.2 • C warming target results in lower oxygen levels compared to preindustrial as discussed in the previous section. The relationships generally hold across the sampled parameter space.
The volume of low-oxygen waters is particularly sensitive to warming and parameter uncertainty (Fig. 5b). We illustrate the sensitivities at the example of the volume of waters with O 2 < 50 mmol m −3 . At the end of the 21st century, there is a 1.7 % increase in this volume per • C of realized warming. Peak increases scale with 63 % • C −1 of equilibrium temperature response. Uncertainties in remineralization cause a spread in this response ranging from 36 to 76 % • C −1 of equilibrium temperature response (90 % confidence range): the longer the remineralization length scale, the higher this sensitivity. Pre-existing low O 2 waters expand and new low O 2 waters may develop in near-bottom environments for higher forcing levels. While the lower temperature targets yield lower volumes of low-oxygen waters, the 9.2 • C target yields higher low O 2 volumes under new steady-state conditions. In brief, hypoxic waters expand over the next millennium across the scenario range and recovery towards modern conditions is slow and, in the case of high forcing, incomplete. Acknowledging millennial timescales, the hazard of expanding low O 2 waters is much larger than when assessed on the near-term timescale.
The sensitivity of changes in low-latitude export production (30 • S-30 • N) is similar at the end of the 21st century and at time of peak decline. Changes scale with 2.2 % • C −1 . Export production in the low latitudes recovers for lower forcing, but remains reduced in the high forcing case.
Decreases in the metabolic index of the upper ocean scale linearly with forcing: 5.1 % • C −1 of realized warming at the end of the 21st century and 6.2 % • C −1 of equilibrium temperature response. Likewise, global mean oceanic temperatures increase by 0.099 • C per degree of realized warming at the end of the 21st century and by 0.56 • C per degree of warming at equilibrium. In conclusion, the compound hazards related to deoxygenation and warming, as indicated by the metabolic stress index, evolve over millennia and increase with increasing anthropogenic forcing and with time.
The magnitude and intensity, as well as the duration of oxygen-related transient hazards, and thus the severity of the hazards, increase with increasing temperature targets. The severity combines magnitude and duration of a hazard into one quantity. It may be defined as the time integral of a hazard. The severity of the hazard of expanding hypoxic waters, for example, corresponds to the area under the scenario curve shown in Fig. 2h (the area enclosed by the null line and the modeled evolution, here until the end of the simulation). Figure 2a, h, and g illustrate that the severity of the three hazards of decreasing mean oxygen concentration, expanding hypoxic waters and reduced export of particulate organic matter providing food for deep sea organisms increases strongly from low to high temperature targets.

Uncertainties in O 2 projections
The pattern and magnitude of simulated global O 2 changes are determined by the response of the overturning circulation. O 2 loss due to less O 2 solubility at higher temperatures also gradually decreases oceanic O 2 . Only a few multimillennial simulations with GCMs currently exist. The response of the overturning circulation on long timescales differs among available model simulations (including EMICs and GCMs). Uncertainties in the equilibrium climate sensitivity additionally impact projections of O 2 loss due to solubility. These uncertainties directly impact projections of oceanic oxygen.
Similar circulation dynamics as simulated here ( Fig. 1e and f) were found by Rugenstein et al. (2016) based on EMIC simulations over 10 000 years with ECBILT-CLIO, which features a dynamic, quasi-geostrophic atmosphere. Schmittner et al. (2008) also found similar AMOC and Indo-Pacific MOC characteristics for their EMIC UVic 2.7, which includes an atmospheric energy balance model with fixed wind fields similar to the Bern3D model, over a 2000-year simulation. Yamamoto et al. (2015), on the other hand, found different overturning characteristics in a simulation with a stateof-the art Earth system model (MIROC 3.2 for a 4 × CO 2 ) over 2000 years. There AMOC slowed down with no recovery, while Antarctic Bottom Water formation decreased only slightly and gradually increased thereafter. Predictions of AMOC have received more attention so far, and AMOC slowdown and partial or full recovery emerges in other multimillennial simulations (Zickfeld et al., 2013;Li et al., 2013;Weaver et al., 2012). AMOC and Southern Ocean overturning in CMIP5 Earth system models were analyzed by Heuzé et al. (2015). They found AMOC and Southern Ocean overturning is positively correlated in most CMIP5 models by the end of the 21st century. Generally, preindustrial circulation states, magnitudes and timing of changes are highly model and scenario dependent such that the long-term evolution of meridional overturning is uncertain. As oxygen changes are dominated by circulation changes, this makes the oxygen prediction highly model and scenario dependent as well. The simulated timing and strength of peak O 2 decrease in Bern3D Earth Syst. Dynam., 9, 797-816, 2018 www.earth-syst-dynam.net/9/797/2018/ is similar to what Schmittner et al. (2008) found (AD 3000, 30 % for SRES A2 high emission scenario/SAT ∼ 10 • C in UVic 2.7). Other comparable simulations show earlier peaks and smaller magnitudes such as AD 2600, 16 % decrease for RCP8.5/ SAT ∼ 7 • C in CLIMBER-3α in Mathesius et al. (2015), and after 800 model years, 10 % for 4 × CO 2 / SAT ∼ 8.5 • C in MIROC 3.2 in Yamamoto et al. (2015). Major physical limitations of our simulations concern prescribed winds and ice sheets. Future model studies may include sensitivity simulations with prescribed changes in the wind stress over the ocean (e.g., Tschumi et al., 2008) and prescribed meltwater fluxes or apply Earth system models with interactive atmospheric dynamics and ice sheets. Our study, as is the case for most climate change simulations, does not include melting of continental ice sheets, which would tend to further (transiently) reduce circulation (Bakker et al., 2016) and increase the equilibrium climate sensitivity.
Current generation GCMs, such as is the case for Bern3D, have difficulty simulating the current distribution of OMZs due to missing physical processes operating at small spatial scales, such as eddies and zonal jets (Cocco et al., 2013;Bopp et al., 2013) or missing biogeochemical characteristics. Large model-data and model-model discrepancies exist . Laufkötter et al. (2017) recently achieved improved representation of OMZs introducing temperature and oxygen dependence of the remineralization profile within a GCM (GFDL ESM2M). In our ensemble, the magnitude of peak increases in low O 2 waters depend strongly on the rate of organic matter remineralization. Temperature-dependent feedback mechanisms, neglected here, may be addressed in future studies. Both particulate sinking speed and local remineralization rates, which control the remineralization profile, have been shown to be sensitive to temperature. While higher temperatures increase bacterial activity and therefore remineralization (Bendtsen et al., 2014) they decrease viscosity and therefore increase sinking speed (Taucher et al., 2014). The net effect on the remineralization profile is correspondingly uncertain. In addition, ecosystem structure influences the size and density of organic particles available for export (Armstrong et al., 2001(Armstrong et al., , 2009. Given these existing uncertainties and the coarseresolution physical models, the projections of OMZs has to be viewed with caution. Despite simulated lower background concentrations of O 2 in the subsurface ocean, the volumes of low O 2 waters decrease for steady-state conditions in the model. We neglect a number of biogeochemical feedback mechanisms that could alter biological productivity in the surface ocean and thus remineralization fluxes in the water column. Any mechanisms that would increase remineralization would tend to decrease oceanic oxygen, and mechanisms that decrease remineralization would increase the oceanic oxygen content. Future studies may address feedbacks from sediment interactions and imbalances from riverine input and burial (such as Roth et al., 2014;Niemeyer et al., 2017), temperature-dependent remineralization and variable stoichiometry. Further investigations may also address nitrogen cycle dynamics and assess the interplay of denitrification and N fixation and of external atmospheric and terrestrial nitrogen sources.

Implications and conclusion
In Bern3D, strong deoxygenation in all basins is projected to peak long after the end of the 21st century, and new steady-state conditions establish after AD 8000 in scenarios where radiative forcing is stabilized in the next century. The equilibration timescale of oceanic oxygen is therefore longer than the thermal equilibration timescale of both the atmosphere (∼ 1000 years) and the ocean (∼ 4000 years). Based on CMIP5 models, Sweetman et al. (2017) discuss the deep sea ecosystem implications of climate change by 2100. Deep sea ecosystems provide a range of services including habitat provision, nursery grounds, trophic support and refugia to biodiversity (reviewed in Sweetman et al., 2017). Biogeochemical changes such as deoxygenation, warming, acidification and less food availability will likely be accompanied by exploitation of mineral resources, overfishing and dumping of pollutants and microplastics. We project the largest biogeochemical changes to occur after 2100 and to aggravate over millennia. How these changes will affect deep sea ecosystems is poorly understood, but the adaptation to stress may be limited by their slow growth rates and long generation times (Sweetman et al., 2017). Figure 6a contrasts near-term (AD 2100) and peak changes (relative to 1870-1899) in definitions of metabolically viable habitats in the upper ocean, hypoxia and food availability as projected by Bern3D for a 1.5 • C warmer world. Export in low latitudes (30 • S-30 • N) as an indicator of food availability is reduced by maximally 4 % over the course of the simulation in this scenario. Median decreases in the metabolic index, representing viable habitat reductions of the upper ocean, amount to 11 % for a 1.5 • warmer world. The volume of low-oxygen waters is particularly sensitive to anthropogenic warming and peak changes occur after the end of the 21st century. The volume of water with O 2 < 50 mmol m −3 changes by 6.6 % by the end of the century and by 14 % at its peak. Meeting the 1.5 • C climate target of the Paris Agreement requires very fast and very stringent emission reductions Sanderson et al., 2016;Millar et al., 2017). Estimates by Steinacher et al. (2013) for a range of scenarios show that post-2017 allowable carbon emissions from fossil fuel need to be lower than 320 Gt C to meet the 1.5 • C target with 66 % probability (320 Gt C is derived by adjusting the emission limit as displayed in Fig. 4 of Steinacher et al., 2013, for  decades of current fossil fuel use. The nationally determined contributions, outlining emission mitigation actions by the parties of the Paris Agreement, need to be more ambitious and broader in scope to meet the 1.5 or 2 • C targets (Joeri et al., 2016). Such efforts would lead not only to lower warming compared to the current emission trajectory, but also have the benefit of reduced marine hazards as investigated here (Fig. 6b). Higher temperature targets increase the hazard of ecosystem impacts as expressed in the investigated variables. In particular, definitions of peak hypoxia exhibit strong sensitivity to additional warming (Fig. 6b). Definitions of deoxygenation, marine food scarcity and marine aerobic habitat reduction are aggravated for the 2 • C compared to the 1.5 • C temperature target and investigated hazards are strongly amplified in a world where surface air temperature is stabilized at 3.3 • C (Fig. 6b). Unbounded use of carbon emissions from existing fossil resources is projected not only to lead to a global warming on the order of 10 • C (Figs. 1a and 5; Randerson et al., 2015;Zickfeld et al., 2013), but also to a peak reduction in global mean oxygen inventory by almost a factor of 2 (Figs. 2a and 5a).
We find a near-linear relationship between impact-relevant marine hazards and global mean surface air temperature. This allows us to quantify avoided hazards per unit of avoided global warming. For example, emission mitigation measures would help to reduce peak O 2 loss by 4.4 % • C −1 of avoided equilibrium warming.
The Earth system response timescale to climate change spans several millennia such that anthropogenic perturbations to greenhouse gas concentrations commit the Earth system to long-term, irreversible climate change (Clark et al., 2016). Our simulations show that the long-term fate of oceanic oxygen is characterized by an initial decline followed by a recovery phase. Peak decline and associated potential adverse ecosystem impacts are projected long after stabilization of radiative forcing in the atmosphere. This adds to the list of long-term Earth system commitments including warming, acidification and sea-level rise assessed elsewhere (Eby et al., 2009;Ridgwell and Schmidt, 2010;Lord et al., 2016;Pfister and Stocker, 2016;Steinacher and Joos, 2016;Clark et al., 2016). Long-term, multi-millennial perspectives are thus required for a full account of climate-related ocean risks.
Data availability. Model output is available upon request to the corresponding author (battaglia@climate.unibe.ch).