Enhanced Organic Carbon Burial in Sediments of Oxygen Minimum Zones Upon Ocean Deoxygenation

Oxygen minimum zones (OMZs) in the ocean are expanding. This expansion is attributed to global warming and may continue over the next 10 to 100 kyrs due to multiple climate CO2-driven factors. The expansion of oxygen-deficient waters has the potential to enhance organic carbon burial in marine sediments, thereby providing a negative feedback on global warming. Here, we study the response of dissolved oxygen in the ocean to increased phosphorus and iron inputs due to CO2-driven enhanced weathering and increased dust emissions, respectively. We use an ocean biogeochemical model coupled to a general ocean circulation model (the Hamburg Oceanic Carbon Cycle model, HAMOCC 2.0) to assess the impact of such regional deoxygenation on organic carbon burial in the modern ocean on time scales of up to 200 kyrs. We find that an increase in input of phosphorus and iron leads to an expansion of the area of the OMZ impinging on continental margin sediments and a significant decline in bottom water oxygen in the open ocean relative to pre-industrial conditions. The associated increase in organic carbon burial could contribute to the drawdown of ~1,600 Gt of carbon, which is equivalent to the total amount of CO2 in the atmosphere predicted for the year 2100 in a business as usual scenario, on time scales of up to 50 kyrs. The corresponding areal extent of sediments overlain by bottom waters with little or no oxygen as estimated by the model is not very different from the minimum area estimated for two major oceanic anoxic events in Earth's past. Such events were associated with major perturbations of the oceanic carbon cycle, including high rates of organic carbon burial. We conclude that organic carbon burial in low oxygen areas in the ocean could contribute to removal of anthropogenic CO2 from the atmosphere on long time scales.

Oxygen minimum zones (OMZs) in the ocean are expanding. This expansion is attributed to global warming and may continue over the next 10 to 100 kyrs due to multiple climate CO 2 -driven factors. The expansion of oxygen-deficient waters has the potential to enhance organic carbon burial in marine sediments, thereby providing a negative feedback on global warming. Here, we study the response of dissolved oxygen in the ocean to increased phosphorus and iron inputs due to CO 2 -driven enhanced weathering and increased dust emissions, respectively. We use an ocean biogeochemical model coupled to a general ocean circulation model (the Hamburg Oceanic Carbon Cycle model, HAMOCC 2.0) to assess the impact of such regional deoxygenation on organic carbon burial in the modern ocean on time scales of up to 200 kyrs. We find that an increase in input of phosphorus and iron leads to an expansion of the area of the OMZ impinging on continental margin sediments and a significant decline in bottom water oxygen in the open ocean relative to pre-industrial conditions. The associated increase in organic carbon burial could contribute to the drawdown of ∼1,600 Gt of carbon, which is equivalent to the total amount of CO 2 in the atmosphere predicted for the year 2100 in a business as usual scenario, on time scales of up to 50 kyrs. The corresponding areal extent of sediments overlain by bottom waters with little or no oxygen as estimated by the model is not very different from the minimum area estimated for two major oceanic anoxic events in Earth's past. Such events were associated with major perturbations of the oceanic carbon cycle, including high rates of organic carbon burial. We conclude that organic carbon burial in low oxygen areas in the ocean could contribute to removal of anthropogenic CO 2 from the atmosphere on long time scales.

INTRODUCTION
Oxygen minimum zones (OMZs) in the ocean are currently expanding, likely due to global warming (e.g., Schmittner et al., 2008;Stramma et al., 2010;Cocco et al., 2013;Schmidtko et al., 2017;Breitburg et al., 2018). The long-term climate change initiated by anthropogenic CO 2 emissions (Archer, 2005;Eby et al., 2009) is expected to continue in the future (Hansen, 2005;Friedlingstein et al., 2014a). Rising temperatures decrease the solubility of oxygen in surface waters and enhance stratification (Bopp et al., 2002;Stramma et al., 2009) thereby reducing the ventilation of the interior of the ocean. A 1 • C warming throughout the upper ocean alone could result in a 3-fold increase in the volume of ocean waters containing <5 µM oxygen (Deutsch et al., 2011). Changes in dissolved oxygen in the ocean have major implications for biogeochemical cycles. Low oxygen levels alter the availability of key nutrients, such as phosphorus, iron and nitrogen (Duce, 1986). This can directly impact the carbon cycle by affecting rates of organic matter production, respiration and burial in the ocean. On time scales of a hundred thousand to a million years, increased burial of organic carbon acts as a sink for CO 2 (e.g., Freeman and Hayes, 1992;Kuypers et al., 1999;Barclay et al., 2010) and provides a negative feedback on global warming. Most studies on ocean CO 2 sinks on short to long-time scales (e.g., Cox et al., 2000;Archer, 2005;Friedlingstein et al., 2014b) focus on ocean uptake of CO 2 through its solubility in seawater, the biological pump, reactions involving CaCO 3 and silicate weathering. The role of organic carbon burial as a sink for CO 2 in the future ocean has received little attention, despite its potential for long-term carbon storage.
Given current trends in fossil fuel emissions and the slow transition to decarbonization, atmospheric CO 2 levels are likely to continue to increase in the near future (Friedlingstein et al., 2014a). After reaching a maximum, CO 2 will decrease within a few centuries, but enough CO 2 will likely persist on timescales of 100 kyrs to affect climate dynamics (Archer, 2005). The impact on dissolved oxygen is therefore expected to continue far into the future. Many studies have assessed the response of OMZs to climate change over the next 100 to 1,000 years (e.g., Oschlies et al., 2008;Hofmann and Schellnhuber, 2009;Stramma et al., 2010Stramma et al., , 2012Cabré et al., 2015), but only a few have looked into the longer-term consequences. In one such study, Watson et al. (2017), used a biogeochemical box model to assess the response of the ocean to anthropogenically driven phosphorus inputs and concluded that this could lead to an entirely oxygen depleted ocean within 2000 years. Results of several studies with general circulation models including marine biogeochemical processes suggest a significant expansion of OMZs in response to changes in either CO 2 -driven changes in phosphorus inputs from weathering, ocean circulation or oceanic temperatures on time scales of 10 to 100 kyrs. However, these models do not show evidence for full scale oceanic anoxia (Palastanga et al., 2011(Palastanga et al., , 2013Beaty et al., 2017;Niemeyer et al., 2017;Kemena et al., 2019).
Besides phosphorus, other nutrients may limit primary productivity, such as iron and nitrogen (Duce, 1986). In fact, most current marine productivity is limited by nitrogen. Whether this will still be the case in the future ocean is uncertain because the nitrogen inventory is strongly affected by ocean deoxygenation (e.g., Deutsch et al., 2007;Gruber, 2011;Oschlies et al., 2019). Moreover, at least 30% of modern marine primary productivity is iron-limited, particularly in the high nitrate, low chlorophyll regions (e.g., Martin et al., 1989Martin et al., , 1990Martin, 1991;Price et al., 1991Price et al., , 1994Moore et al., 2001Moore et al., , 2013. Aeolian dust is a primary source for iron. Dust inputs are largely driven by atmospheric moisture, wind patterns, rainfall, soil moisture and vegetation cover (Duce, 1995), all of which are expected to change in the future. Despite large modeling efforts, predicting future changes in dust inputs remains a challenge Evan et al., 2014;Kok et al., 2018) and global warming may lead to a decrease or an increase in the input of aeolian dust (e.g., Harrison et al., 2001;Mahowald and Luo, 2003;Tegen et al., 2004;Mahowald et al., 2005). The most recent hypothesis is that rates of modern global dust emission may double in at least the next 100 years due to an increase in wind speed and a decrease in soil moisture in key source regions (Kok et al., 2018). Iron inputs from dust lead to important shifts in phytoplankton species and sizes, and alter nitrogen to iron ratios in the ocean. At least 30% of modern marine primary productivity is iron-limited, particularly in the high nitrate, low chlorophyll regions (e.g., Martin et al., 1989Martin et al., , 1990Martin, 1991;Price et al., 1991Price et al., , 1994Moore et al., 2001Moore et al., , 2013. Consequently, changes in the oceanic iron cycle may lead to variability in global rates of ocean productivity, dissolved oxygen concentrations and organic carbon export (Moore et al., 2001(Moore et al., , 2013Altabet et al., 2002). Besides dust, continental margin sediments may also act as a key source of iron to ocean waters, and enhanced future supply to surface waters is expected as OMZs expand (Severmann et al., 2008;Homoky et al., 2012;Henderson et al., 2018).
At present, OMZs have a size of 100 × 10 6 km 3 , of which a tenth is considered the 'core' of the OMZs, characterized by waters with an oxygen concentration of < ∼20 µM (Paulmier and Ruiz-Pino, 2009). Permanent OMZs with suboxic cores are located in the Eastern South Pacific, Eastern (Sub-) Tropical North Pacific, Arabian Sea and Bay of Bengal, and account for 0.8% of the global ocean's volume. In all areas, except the Eastern Tropical North Pacific, these OMZs also include large areas with permanently anoxic (or nearly anoxic) bottom waters (e.g., Stramma et al., 2008;Paulmier and Ruiz-Pino, 2009). The total area where suboxic or anoxic waters of permanent OMZs impinge on the seafloor is estimated at ∼1.15 × 10 6 km 2 (∼0.3% of the total ocean floor). Anoxic waters account for 70% of this area (Helly and Levin, 2004). Note that the low oxygen area in the Eastern Atlantic is not considered an OMZ here because its core does not contain <20 µM of oxygen. Since the early 1960s, anoxia has also increased on continental shelves (e.g., Tropical Pacific and Atlantic, Northern Gulf of Mexico, and in coastal waters of Eastern USA, Western and Northern Europe, Korea, China, Taiwan and Japan) and in semi-enclosed basins (e.g., the Baltic Sea, Black Sea, and Mediterranean Sea) as summarized in Stramma et al. (2010) and Breitburg et al. (2018). Due to limitations in spatial resolution, such near coastal areas are not considered in this study.
Oxygen minimum zones play an essential role in the oceanic nitrogen cycle, by contributing to removal of nitrogen through denitrification and anammox and to the release of nitrous oxide, a powerful greenhouse gas (Codispoti et al., 2001). Denitrification requires low oxygen concentrations and is observed in sediments and the water column of OMZs where oxygen concentrations are <20 µM (Smethie, 1987). With the exception of the OMZ in the Eastern Sub-Tropical North Pacific, all permanent OMZs are associated with a denitrification zone. Together, these zones account for a third of the total volume of all OMZ cores (∼3.4 million km 3 ; Paulmier and Ruiz-Pino, 2009). The area where benthic denitrification occurs could be larger than that of water column denitrification (DeVries et al., 2013). Moreover, OMZs favor benthic recycling of nutrients (Fisher et al., 1982). In the sediments of OMZs, rates of organic carbon burial are typically enhanced due to a combination of a high primary productivity and increased preservation of organic matter under low oxygen conditions (Emerson, 1985;Arthur et al., 1998;Hedges et al., 2001).
In the geological past, episodes of widespread ocean anoxia that lasted for more than 100 kyrs also led to changes in organic carbon burial (Schlanger and Jenkyns, 1976). Key examples of such episodes are the mid-Cretaceous Oceanic Anoxic Event 2 (OAE 2; 94 Myr ago) and the Toarcian Oceanic Anoxic Event (T-OAE, 183 Myr). Both events occurred during a period of frequent volcanic activity, high concentrations of greenhouse gases in the atmosphere and global warming (e.g., Schlanger et al., 1987;Berner, 1992;Hesselbo et al., 2000;Huber et al., 2002;Wilson et al., 2002;Kemp et al., 2005). During OAE 2, ∼5% of the global sea floor was covered by anoxic and sulfidic bottom waters (Owens et al., 2013(Owens et al., , 2016Dickson et al., 2016). This area is thought to have been mainly limited to the proto-North Atlantic, including various of its continental margins (e.g., Schlanger and Jenkyns, 1976;Ruvalcaba Baroni et al., 2014;van Helmond et al., 2014) whereas the degree of oxygenation of the Pacific Ocean is uncertain. OAE 2 was attributed to a reduced solubility of oxygen due to high sea surface temperatures (Barron et al., 1993) and changes in ocean circulation, enhanced nutrient supply and recycling (e.g., Parrish and Curtis, 1982,?;Föllmi, 1999;Mort et al., 2007;Wagner et al., 2007;Kraal et al., 2010;Trabucho Alexandre et al., 2010;Topper et al., 2011;Poulton et al., 2015). Geochemical analysis revealed high rates of organic carbon burial in both the open and especially the coastal ocean (Owens et al., 2018) with an estimated removal of 26% of atmospheric CO 2 during OAE 2 (Barclay et al., 2010).
The exact spatial extent of the anoxia during the T-OAE is less well-known than that for OAE 2. However, geological evidence and model results suggest that during T-OAE protracted anoxia/euxinia was geographically widespread (Gill et al., 2011;Them et al., 2018), but likely mainly confined to relatively shallow areas (e.g., van de Schootbrugge et al., 2005;Dickson et al., 2017). These areas account for at least ∼0.3% of the seafloor (Ruvalcaba Baroni et al., 2018) and provide evidence for high rates of organic carbon burial. Trace metal and thallium isotope analysis further suggest that bottom water deoxygenation did not remain constant over the event (e.g., Montero-Serrano et al., 2015;Them et al., 2018). Organic carbon burial during T-OAE is thought to have contributed to at least 8 to 14% of the drawdown of atmospheric CO 2 (Xu et al., 2017;Ruvalcaba Baroni et al., 2018). Both OAE 2 and T-OAE share similar biogeochemical characteristics, such as prolonged periods of regional photic zone euxinia, major losses in marine biodiversity, large perturbations in most elemental cycles and high rates of organic carbon burial (e.g., Schlanger and Jenkyns, 1976;Jenkyns, 1985Jenkyns, , 2010McArthur et al., 2008;Dickson et al., 2017).
In this study, we use a biogeochemical model coupled to a general ocean circulation model to assess changes in oxygen and organic carbon dynamics in the modern ocean in response to increased CO 2 -driven weathering inputs of phosphorus and increased inputs of iron from dust on time scales of up to 200 kyrs. We discuss the effects on ocean deoxygenation, export production and burial of organic carbon and the potential impact on CO 2 drawdown from the atmosphere. Our study indicates that long term changes in oxygen concentrations and organic carbon burial in the modern ocean upon continued CO 2emissions could be comparable to those inferred for periods of widespread oceanic anoxia in the geological past.

Model Description
We use the Hamburg Oceanic Carbon Cycle model (Maier-Reimer et al., 1993) in its annually averaged version (HAMOCC 2.0; Heinze et al., 1999Heinze et al., , 2003Heinze et al., , 2006, which is specifically designed for long-term integrations. The physical model has a horizontal resolution of 3.5 • × 3.5 • and includes 11 layers for the oceanic water column with a total ocean surface of 359 301 107 km 2 , a zonally mixed atmosphere compartment and a bioturbated sediment layer of 10 cm, separated into 11 sublayers. The model assumes a fixed modern ocean circulation and is based on the velocity and thermohaline fields of the Hamburg Large-Scale Geostrophic ocean general circulation model (Maier-Reimer et al., 1993). The biogeochemistry of the model describes the marine carbon, silicon, oxygen and phosphorus cycles (Heinze et al., 1999(Heinze et al., , 2003(Heinze et al., , 2006. Sedimentary phosphorus cycling and the oceanic iron cycle were implemented by Palastanga et al. (2011) (P-model) and Palastanga et al. (2013) (PFe-model), respectively. In the later study, the PFe-model was used to study glacial-interglacial variability. Here, we use their model version (that includes both phosphorus and iron) and their preindustrial configuration and apply it to study the modern ocean. For a full description of the model, we refer to these earlier papers.
The model includes both aerobic and anaerobic degradation of organic matter, assuming first-order kinetics, with the anaerobic pathway becoming active when oxygen is nearly depleted (i.e., below 5 µM). Rate constants for aerobic and anaerobic degradation of organic matter vary greatly between different depositional environments (Tromp et al., 1995). Here, rate constants for the aerobic and anaerobic degradation of organic matter were selected based on the best fit of particulate organic carbon (POC) in the surface sediments (Seiter et al., 2004). Differences in values of the first-order rate constants for continental margins (aerobic: 0.01 yr −1 ; anaerobic: 0.008 yr −1 ) and the deep sea (aerobic: 0.005 yr −1 ; anaerobic: 0.002 yr −1 ) reflect the lower quality of organic matter reaching the seafloor in the latter environment. The model pattern of POC in sediments is in good agreement with global observations of organic carbon contents in surface sediments as demonstrated by Palastanga et al. (2011). Key biogeochemical processes related to sediment phosphorus cycling in the model are degradation of organic phosphorus, precipitation of authigenic carbonate fluorapatite and binding of dissolved phosphate to and release from iron-oxides. The sedimentary iron cycle includes both formation and dissolution of iron oxides. When oxygen falls below 5 µM, anaerobic diagenesis becomes active and both iron and phosphorus are released to the porewater upon dissolution of the iron-oxides. The model simulates mixing and burial of all solid phases and diffusion of solutes. While the release of dissolved phosphate from the sediment to the overlying water is explicitly modeled, that of iron is prescribed using the parameterization of Moore and Braucher (2008). In the model, preferential regeneration of phosphorus in the open ocean accelerates the expansion of anoxia mainly in the eastern and western tropical Pacific Ocean, and the southeastern Pacific sector (Palastanga et al., 2011). Results of sediment phosphorus profiles for the pre-industrial ocean are in good agreement with observations in the open ocean and along continental margins (e.g., Wheat et al., 1996;Harrison et al., 2005;Parekh et al., 2005). For further details on the implementation of the sedimentary processes we refer to the work of Palastanga et al. (2011) and Palastanga et al. (2013) and key references therein (e.g., Ruttenberg and Berner, 1993;Filippelli and Delaney, 1996;Slomp et al., 1996). The model is forced with the preindustrial fields of annual mean ocean circulation from Winguth et al. (1999) and annual mean dust deposition from Mahowald et al. (2006), with a total input of 0.35 Tmol Fe yr −1 (Table 1) as described in detail by Palastanga et al. (2013).

Experimental Design
We performed 5 experiments with the model (Table 1): a preindustrial simulation (Palastanga et al., 2013), an experiment with a 2-fold increase of the input of phosphorus (2xP), an from dust increased by a factor of 2.5 and 4, respectively experiment with a 2-fold increase of the input of iron from dust (2xFe), an experiment with a 2-fold increase of the input of phosphorus from rivers and of iron from dust (2xPFe), the same but then with a 2.5-fold increase of both riverine phosphorus and dust iron (2.5xPFe), and finally, an experiment with 2.5-fold increase in phosphorus and 4-fold increase in iron (2.5xP4xFe). The pre-industrial simulation was run over 200 kyrs to ensure near-equilibrium state considering that the residence time of phosphorus in the ocean is about 120 kyrs. The other 4 experiments started from the pre-industrial scenario but with a stepwise increase in the nutrient supply and were further run over 200 kyrs. To assess the role of limitation of oceanic productivity by phosphorus alone vs. limitation by phosphorus and iron upon increased nutrient inputs, we compared the oxygen concentrations and POC export from our 2xPFe experiment to that presented for a doubling of the phosphorus input by Palastanga et al. (2011) (Supplementary Material S1). Due to different model characteristics, such as its horizontal and vertical low resolution, its yearly integrated forcing, its simplified topography and its assumed fixed rate constants for organic carbon degradation, complete loss of oxygen near the seafloor is strongly underestimated in HAMOCC 2.0. Therefore, we consider anoxia to occur when the model oxygen concentrations are close to zero (i.e., below 12 µM) and hereafter use this approximation. Consequently, the terms hypoxia and suboxia here refer to oxygen concentrations that are between 60 and 20 µM and between 20 and 12 µM, respectively. The range of phosphorus inputs from rivers tested here is in line with estimates of the change in rates of weathering for OAE 2 (Blättler et al., 2011). The iron concentrations in dust follow the hypothesis of a future increase in global dust emission as simulated by various Earth System models when forced by a new emission scheme that best captures dust emission measurements (Kok et al., 2018). Note that our 2.5xP4xFe scenario simulates an ocean where primary productivity is less limited by iron than in all other scenarios.

The Pre-industrial Ocean
In our pre-industrial simulation (Figures 1A,B), the total volume of the suboxic and anoxic waters of the OMZs are 25.6 × 10 6 and 3.4 × 10 6 km 3 ( Table 2). These volumes account for ∼2% and ∼0.3% of the total ocean volume, respectively. When compared to recent volume estimates for modern open ocean waters containing <20 µM of O 2 (∼10.3 × 10 6 km 3 with an average thickness of 340 m-i.e., accounting for 0.8% of the ocean volume; Paulmier and Ruiz-Pino, 2009), our preindustrial simulation overestimates the volume of suboxic waters by nearly a factor of 3. Anoxic waters are present between water depths of 400 and 900 m over an area of ca. 9.6 × 10 6 km 2 , representing about 3% of the total ocean surface (Figures 1A,B and Table 2). Hypoxic bottom waters containing <60 µM of O 2 are present over an area of 3.1 × 10 6 km 2 (∼0.6% of the ocean floor), of which 0.1 × 10 6 km 2 contains <∼20 µM of O 2 in the model (∼0.03% of the ocean floor). The areal extent The percentage of the ocean's surface and volume that is anoxic in each scenario is calculated assuming a total ocean area and volume of 359 × 10 6 km 2 and 1,349 × 10 6 km 3 , respectively. The factor by which the anoxic volume is increased when compared to the pre-industrial simulation. The results of the 2xP experiment were similar to those of the Preindustrial experiment and therefore, are not shown. of naturally occurring suboxic bottom waters on continental margins is on the order of ∼1.15 × 10 6 km 2 (Helly and Levin, 2004), which is an order of magnitude higher than what we find in the model. Hence, the model overestimates the volume of the OMZs but underestimates the area where it impinges on the continental margins. Note however, that the bottom area considered in Helly and Levin (2004) may include some shallow areas that cannot be captured by our model resolution.

Model Response to a Doubling of Phosphorus and Iron Inputs
Upon a 2-fold increase in both phosphorus and iron inputs for 200 kyrs, a significant expansion of low oxygenation occurs when compared to pre-industrial conditions (Figures 1A-D). This expansion is even more pronounced when doubling the riverine phosphorus input in the P-model of Palastanga et al. (2011), as iron in all our scenarios acts as the limiting nutrient for POC export production in most of the ocean (Supplementary Material S1; Figures S1.1, S1.2). Otherwise well-oxygenated intermediate waters in the Southern and Western Pacific, and the Atlantic Ocean show a decrease in oxygen concentrations of ∼40 to 70 µM ( Figure 1E). Bottom water oxygen concentrations also decrease, with the largest change (of up to 80 µM of O 2 ) observed for tropical and equatorial continental margins (Figure 1F). The most expanded OMZ is that of the Pacific Ocean, where a large anoxic core develops throughout most of the tropical and subtropical regions in the Eastern Pacific ( Figure 1C). The volume of suboxic waters with <20 µM increases by a factor of 4.4 relative to the pre-industrial scenario (with a total volume of 128 x 10 6 km 3 ), whereas the volume of the anoxic core increases by a factor of nearly 10 ( Table 2). However, the anoxic volume remains confined to a relatively thin intermediate layer spanning from 450 to 1,000 m ( Figure 1D and Figure S2). The sum of all anoxic waters in the 2xPFe experiment covers about 13% of the ocean surface ( Figure S2a) with an area of 47 × 10 6 km 2 ( Table 2). Low-oxygen conditions in bottom waters increase and anoxic bottom waters now develop covering about 0.3% of the ocean floor ( Table 2).
The spatial distribution in organic carbon burial in our 2xPFe scenario is similar to that in the pre-industrial simulation and is characterized by high burial rates on continental margins (Figures 2A,B). In the 2xPFe scenario, on average ca. 33% more organic carbon is buried than during the pre-industrial at a global rate of 0.045 Gt of carbon per year. Maximum carbon contents in sediments of the 2xPFe scenario are 4 wt% (vs. <3 wt% in the pre-industrial run), and are found in sediments of equatorial margins, along the northern margins of the Indian Ocean and eastern margins of the tropical Pacific Ocean, in the Gulf of Mexico and along the margins of northern Europe. In the 2xPFe scenario, the average organic carbon content of sediments at the global scale is ∼0.8 wt%, while that for continental margins is 1.5 wt% and that for the open ocean is 0.7 wt% ( Figure 2C). In all our scenarios with high nutrient inputs, the increase in organic carbon content in sediments of the open ocean is larger than that of the continental margins when comparing to pre-industrial conditions (e.g., Figure 2D).

Model Sensitivity to Increased Iron and Phosphorus Input
In an ocean where iron is limiting in 85% of the surface ocean, as is the case in the pre-industrial ocean in the model used here, increased iron input leads to an enhanced primary productivity. When, for example, the iron input from dust alone is doubled (2xFe), the area where primary productivity is limited by iron decreases to 65% of the surface ocean. Consequently, POC export production increases by 30% within 10 kyrs. However, a rise in iron alone cannot sustain this production and the iron consumption by primary producers and scavenging by sinking particles (see Palastanga et al. (2013)) leads to a gradual decrease of the iron concentrations by 10% at 200 kyrs ( Figure 3A). Changes in oxygen demand follow those in POC export production and global oxygen concentrations rise slightly after an initial decline ( Figure 3B). Results of the 2xP experiment are not shown because the lack of iron then limits primary productivity, resulting in rates similar to pre-industrial values. The response to an increase in the input of one single nutrient clearly contrasts with the rise in POC export and decline in oxygen observed when both iron and phosphorus inputs to the ocean are doubled, as in our 2xPFe scenario (Figures 3A,B). The rise in POC export and decline in oxygen is even more pronounced in the scenarios where nutrients are increased by a factor >2 (2.5xPFe and 2.5xP4xFe; Figures 3A,B). The associated increase in oxygen loss below the surface waters in the ocean leads to a major expansion of the OMZs (Figure S2b: depending on the exact scenario, between 1% to 8% of the ocean volume is anoxic at a depth of 1 000 m) and anoxic bottom waters cover between 0.3 to 1.3% of the seafloor (Figure 3C and Table 2). Increased POC export production and oceanic anoxia promote the burial of organic carbon. The cumulative organic carbon burial over 200 kyrs in the scenarios with both increased phosphorus and iron input ranges from 8,547 to 11,085 Gt of carbon ( Figure 3D). In all our model experiments with increased phosphorus input, the average phosphorus concentration in the ocean increased during the experiment. This is because iron is limiting in large parts of the ocean and phosphorus is only partially consumed. To illustrate the consequences for the biogeochemistry of the ocean, we compared the results of the 2xPFe scenario (i.e., performed with the PFe-model) to those of the scenario in Palastanga et al. (2011) with a doubling of riverine phosphorus input (i.e., performed with the P-model). We find that in the latter model, the same increase in phosphorus input (i.e., a doubling), when compared to the pre-industrial ocean, leads to lower oxygen concentrations, higher POC export and higher rates of organic carbon burial after 200 kyrs (Figure S1.1).

Increased Nutrient Availability as a Driver of Ocean Deoxygenation
In our model simulations, increased inputs of iron and phosphorus led to enhanced export production and ocean deoxygenation. Importantly, phosphorus in surface waters was not fully consumed by primary producers and, over time, enhanced riverine phosphorus inputs led to a significant increase in the oceanic phosphorus reservoir. Other biogeochemical ocean models also show such an increase in the global phosphorus inventory in response to increased riverine supply (e.g., Palastanga et al., 2011;Niemeyer et al., 2017;Watson et al., 2017;Kemena et al., 2019). However, the time-scale, magnitude, causes and consequences of the increase in oceanic phosphorus differs between models. Typically, the increase in phosphorus is largest in models where other nutrients, such as nitrogen and iron are also limiting, because this prevents full phosphorus uptake in surface waters. This is illustrated by our comparison of the model results of the P-model of Palastanga et al. (2011) to those of the PFe-model in this study (Figure S1.2).
In models that consider phosphorus as the only limiting nutrient, primary productivity mainly follows phosphorus inputs. This may lead to an overestimation of ocean deoxygenation and organic carbon burial.
In our simulations, the productivity in most of the surface ocean is limited by iron, and any additional iron supply is rapidly consumed by primary producers. In the real ocean, spatial differences in iron to carbon ratios and, therefore, in iron limitation depend on the phytoplankton species and size (e.g., Moore et al., 2001Moore et al., , 2013. These effects are not included in our model, therefore, the corresponding detailed spatial variability in iron limitation and primary productivity is not captured in our results. However, future changes in the dust field are highly uncertain and our results show how sensitive the ocean is to dust input in the model: a 2-fold increase leads to a rise in global primary productivity and at least a 3-fold increase in the anoxic water volume in the ocean ( Table 2).
Nitrogen also controls primary productivity (Duce, 1986;Tyrrell, 1999;Moore et al., 2013). A direct consequence of an OMZ expansion, such as that observed in our scenarios, is the spreading of its associated denitrification zone (Paulmier and Ruiz-Pino, 2009). The large increase in the anoxic core of the OMZs in our simulations (e.g., of a factor 10 in the 2xPFe scenario; Figure 1 and Table 2) suggests that denitrification in the water column indeed would be enhanced. Benthic denitrification could also increase due to the decline in bottom water oxygen (Middelburg et al., 1996;DeVries et al., 2013). Enhanced N-loss from the ocean and iron supply to the surface waters favor N 2 -fixation which could more than compensate for the N-loss, especially on long time scales (Codispoti, 1989;Tyrrell, 1999). However, whether such a compensation will actually occur in the future ocean is uncertain. If we would assume that a doubling of phosphorus and iron inputs from weathering and dust, respectively, as implemented in our model, is a realistic scenario for the future ocean and that phosphorus and iron are the key limiting nutrients, this would imply a major expansion of low oxygen areas in the ocean on time scales of up to 200 kyrs. In our model simulations, low-oxygenated waters in the open ocean expand and spread from existing OMZs. Suboxic and anoxic bottom waters may cover 2 and 0.3% of the ocean's sea floor, respectively. Given that our model underestimates bottom water anoxia, these percentages are rather conservative. Such an oxygen depleted ocean unavoidably leads to major changes in ocean biogeochemistry and to a reduced biodiversity (e.g., Chameides and Perdue, 1997;Vaquer-Sunyer and Duarte, 2008;Wright et al., 2012).

Role of Oxygen Demand vs. Oxygen Supply
Oxygen concentrations in the ocean depend on the balance between oxygen demand and supply. In our model scenarios, we specifically focused on changes in oxygen demand in intermediate and bottom waters due to nutrient-driven increased export production. Recently, Beaty et al. (2017) used the same model (PFe-model) to assess changes in oxygen supply over 10 kyrs due to increased radiative CO 2 forcing (2-8x pre-industrial pCO 2 levels), in scenarios where the associated temperature rise affected oxygen solubility, while keeping POC export at pre-industrial values. Comparison of our results for the first 10 kyrs of the scenario with a doubling of nutrients (2xPFe) with these scenarios from Beaty et al. (2017) shows that only their most extreme scenario (8xCO2; i.e., an atmospheric pCO 2 of 2238.24 ppmv and change in sea water temperature of 11.5 C • ) results in a less-well oxygenated ocean than in our scenario (Figures 4A-D). This illustrates that while a reduced solubility due to warming of surface waters is a key driver of ocean deoxygenation (Bopp et al., 2002;Plattner et al., 2002;Keeling et al., 2009), changes in nutrient supply have the potential to become equally important on long time scales.
A striking difference between the results of the 2xPFe and 8xCO2 scenarios is that in the latter scenario bottom water oxygen concentrations in the Pacific Ocean are up to ∼30 µM lower (Figures 4D,E,G). This suggests that bottom water concentrations in the ocean are particularly vulnerable to changes in oxygen supply. In a set of additional scenarios, FIGURE 4 | POC export production (A-C) and dissolved oxygen concentrations (D-F) after 10 kyrs for our scenario with a doubling of nutrients (2xPFe) and 2 scenarios of Beaty et al. (2017): one with a 8x increase in CO 2 radiative forcing (8xCO2) and one with a reduction of 25% in the ocean ventilation (25%vent). Difference in oxygen concentrations (G,H) between our scenario and those of Beaty et al. (2017). Positive values indicate that changes in nutrient-driven oxygen demand are smaller than those in temperature-or ventilation-driven oxygen supply, hence allowing higher concentrations of oxygen to remain. Negative values indicate the opposite. Beaty et al. (2017) also assessed the response of ocean oxygen to a reduced ventilation, without changing the circulation pattern or the mixing due to diffusion. The results of these scenarios reveal that reducing ocean ventilation by up to 75% has a comparatively much smaller effect on oxygen concentrations than changes in oxygen solubility or POC export (Beaty et al., 2017). To illustrate this, we show the results of the scenario for a 25% reduction in ventilation from Beaty et al. (2017) and the difference with our 2xPFe scenario (Figures 4C,F,H). Overall, our results suggest that, on long time scales, increased nutrient input with or without a reduced oxygen solubility could lead to expanded OMZs and low oxygen concentrations in bottom waters, particularly in the Pacific Ocean.
Recent results from an earth-system model simulating global warming on millennial time scales show a gain of ∼6% of dissolved oxygen concentrations in the ocean with respect to modern conditions . This gain of oxygen is largely explained by an increase in oxygen in deep waters due to a reduction in the typical aerobic remineralization pathway (which consumes oxygen) which is replaced by the anaerobic denitrification pathway. Additionally, earth-system models under radiative forcing typically show an oxygen gain in the ocean after several thousand years [e.g., about 16% in Ridgwell and Schmidt (2010)]. This has been attributed to warming-induced instabilities of the stratification and its final break down in the Southern Ocean surface. While earthsystem models can account for changes in marine oxygen due to both warming and biogeochemistry, they do not yet account for enhanced nutrient supply due to CO 2 -induced weathering. When comparing our equilibrium results of oxygen concentrations in bottom waters for the 2xPFe and the preindustrial scenarios (Figure 3B), we obtain a loss of about 25% of oxygen in the global ocean when a doubling of nutrient is applied. This highlights that the loss of oxygen caused by enhanced nutrient supply due to CO 2 -induced weathering is of the same order of magnitude as the possible gain of oxygen shown by earth-system models. Because all these drivers occur simultaneously but on different time scales, it is possible that on a given time the oxygen loss from enhanced weathering exceeds the oxygen gain. In such case, enhanced nutrient supply due to CO 2induced weathering could lead to widespread deoxygenation, explaining the occurrence of ocean anoxic events.

Impact of Ocean Deoxygenation on Long Term Organic Carbon Burial in Sediments
The modern atmosphere currently contains about ∼870 Gt of carbon in the form of CO 2 (Pachauri et al., 2014) and this is expected to increase to ca. 1600 Gt by the year 2100 (in a business-as-usual model scenario; Cox et al., 2000). Atmospheric CO 2 is removed from the atmosphere by terrestrial and oceanic processes that can act on a range of time-scales and are sensitive to a range of environmental conditions. For example, the solubility of CO 2 in seawater decreases with rising sea surface temperatures and is expected to turn the ocean into a less efficient sink for CO 2 on decadal time scales (Sabine et al., 2004). In contrast, the rate of silicate weathering is expected to increase upon rising pCO 2 and this sets the ultimate maximum duration of the anthropogenic carbon cycle perturbation, acting on timescales of ∼100 kyrs (Sundquist, 1991).
Our scenarios show that increased organic carbon burial on continental margins and in the deep ocean could potentially bury 1600 Gt of carbon within 50 ± 10 kyrs (Figure 2 and Figure 3D). This implies that organic carbon burial alone could "permanently" remove all CO 2 present in the atmosphere in the year 2100, when assuming no other sources and sinks. In reality, the atmosphere is continuously affected by CO 2 fluxes and associated climate feedbacks. The model results suggest that 20 ± 10 kyrs are required to sequester the carbon from anthropogenic emissions (about half of the 1600 Gt) through organic carbon burial alone. The model results also suggest that, besides continental margins, sediments in the deep sea also contribute to organic carbon burial (Figure 2). Since our model does not include near-coastal areas, organic carbon burial on continental margins is likely underestimated.
In the model, the burial of organic carbon in the sediment depends on the input of POC from the overlying water and rates of organic carbon degradation in the sediment. The POC input depends on its export from the surface ocean and hence is critically dependent on nutrient availability, i.e., the assumed increase in river input of phosphorus and dust input of iron. This nutrient input must be sustained, otherwise productivity, ocean oxygen and rates of organic carbon burial will return to pre-industrial values. For example, if nutrient inputs are abruptly returned to pre-industrial inputs, after reaching equilibrium in the 2xPFe scenario, waters at 700 m and 1000 m depth reoxygenate within ∼2000 and 500 years, respectively ( Figure S3). While a 2-fold sustained increase of phosphorus input due to weathering is in line with the suggested response to high pCO 2 based on the geological record (Blättler et al., 2011), the long term changes in input of iron are uncertain.
In the model, anaerobic degradation of organic matter is assumed to be a factor 0.8 and 0.4 slower than aerobic degradation in continental margin and deep sea sediments, respectively. This results in a global increase in POC burial of at least ∼30% due to POC preservation under low-oxygen conditions alone. These are relatively conservative estimates for the change in degradation upon oxygen depletion when compared to data for modern systems and other model studies (Canfield, 1994;Hartnett et al., 1998;Moodley et al., 2005;Reed et al., 2011). This implies that our estimates of POC burial could be minimum values. We conclude that enhanced organic carbon burial upon ocean deoxygenation could contribute to long term removal of atmospheric CO 2 .

Comparison to Periods of Ocean Anoxia in Earth's Past
In all our scenarios with increased phosphorus and iron inputs, existing OMZs expanded and between 0.3 and 1.3% of the seafloor was covered by anoxic bottom waters after 50 kyrs. This compares well to the known areas with anoxic bottom waters during the Toarcian-OAE and the Cretaceous OAE 2, during which at least ∼0.3% (Dickson et al., 2017;Ruvalcaba Baroni et al., 2018) and ∼5% (Owens et al., 2013(Owens et al., , 2016Dickson et al., 2016) of the global ocean sea floor was overlain by waters depleted in oxygen, respectively. High atmospheric pCO 2 and high rates of nutrient input, sustained over a longer period of time, similar to our scenario for the modern ocean, have been deduced for these past OAEs (e.g., Schlanger and Jenkyns, 1976;Poulsen et al., 2001;Cohen et al., 2004Cohen et al., , 2007Kemp et al., 2005;Barclay et al., 2010;Goddéris et al., 2012;van Bentum et al., 2012;Poulton et al., 2015).
An important difference between our model projections and conditions inferred for the past ocean relates to the location and intensity of the anoxia. In our model results, anoxic bottom waters are observed mostly along continental margins, in particular in the eastern Equatorial Pacific Ocean (Figure 1). In our simulations (and those of Beaty et al., 2017), the OMZ expands and bottom water anoxia is most pronounced at intermediate water depths. In such a setting, widespread sulphidic conditions in the water column are unlikely to develop. In contrast, during the Toarcian and Cretaceous OAEs, most of the water column, i.e., from the seafloor to the bottom of the photic zone, was anoxic and sulphidic (e.g., Sinninghe Damsté and Köster, 1998;Dickson et al., 2017). This difference in vertical extent and intensity of the anoxia is likely due to the differences in continental configuration and ocean circulation, allowing parts of the ocean to be more prone to anoxia in the past than today (Donnadieu et al., 2006;Goddéris et al., 2012). During the Toarcian, this holds specifically for the northern European Epicontinental Shelf, which was subject to a flow pattern that promoted the development of anoxia (Ruvalcaba Baroni et al., 2018). During OAE 2, this was the case in the proto-North Atlantic (e.g., Sinninghe Damsté and Köster, 1998;Kuypers et al., 2002;Owens et al., 2013). Sedimentary rock records indicate that large amounts of organic-rich sediments were deposited in low oxygen areas in the ocean during the T-OAE and OAE 2 (Jenkyns, 1985(Jenkyns, , 2010. This is thought to have led to sequestration of atmospheric CO 2 and may have contributed to the termination of these past events (Sinninghe Damsté and Köster, 1998;van Bentum et al., 2012;Xu et al., 2017;Ruvalcaba Baroni et al., 2018).
In our model simulations, we find that 1600 Gt of carbon can be sequestered in the modern ocean during a period of 50 kyrs (Figure 2). To assess whether this is realistic we can compare these values to the estimates of organic carbon burial for the T-OAE and OAE 2. During the T-OAE, the suboxic and anoxic areas on the European shelf alone are thought to have buried 2400 Gt of carbon during a period of 300-500 kyrs (Ruvalcaba Baroni et al., 2018). We note that nothing is known about the burial in other coastal settings and in open ocean environments at this time, hence, the total burial during the T-OAE will likely have been larger. More is known about the burial of organic carbon during OAE 2 and a recent calculation suggests it amounts to at least 41,000 Gt of carbon, as estimated by Owens et al. (2018) for an approximated event duration of 500 kyrs. In our scenarios, the average rate of organic carbon burial in the modern ocean (here, approximated as the final cumulative amount of organic carbon divided by the time of the run) range from ∼0.032 to 0.055 Gt C yr −1 (Figure 3D). Our results compare well to the average rate of organic carbon burial estimated during, for example, OAE 2 (41,000 Gt C/500 krsy ≈0.082 Gt C yr −1 ), and therefore, are not unrealistic. The difference in the average rate of organic carbon burial in the modern ocean vs. that during OAE 2 may be explained by the widespread presence of anoxic and sulphidic bottom waters during the event, which will have promoted preservation of organic matter through sulfidization (e.g., Meyers, 2007;Raven et al., 2018) and, possibly, a higher productivity and hence higher input of organic matter to the sediment Recent findings on OAEs reveal a time lag between marine deoxygenation and significant organic carbon accumulation in sediments (Owens et al., 2016;Ostrander et al., 2017;Them et al., 2018). This suggests that the contribution of organic carbon burial to atmospheric CO 2 sequestration became effective only upon widespread deoxygenation and the establishment of anoxic bottom waters. Our scenarios do not show such a time lag. This difference is related to differences in nutrient availability: while in our scenarios, nutrient inputs are elevated from the start of the simulation, in reality, nutrient input are expected to gradually increase and the timing of that increase can determine this time lag. However, we conclude that the rates of organic carbon burial in our scenario are realistic when compared to those for the geological record.
In conclusion, anthropogenic CO 2 emissions as considered in the "business as usual" projection (i.e., no reduction of fossil fuel emissions for the next 100 years) could drive the modern ocean toward a state of long-term widespread anoxia and increased organic carbon burial. However, the time scale of recovery from anoxia may be different from that during an OAE. This is because anthropogenic emissions represent a single pulse of CO 2 on geological time scales, whereas during the T-OAE and OAE 2, volcanic inputs replenished the CO 2 in the atmosphere. Thus, if no new CO 2 sources are activated, organic carbon burial could be more efficient as a sink for CO 2 in the future ocean than it was in the past. However, this is highly dependent on the response of other natural carbon sources to global warming, such as, for example, the evasion of methane from permafrost, freshwater systems or clathrates, since this will affect the total carbon emission to the atmosphere.

CONCLUSIONS
Using a biogeochemical model coupled to a general circulation model for the modern ocean, we show that a long-term increase in nutrient inputs could lead to a major increase in the area and volume of suboxic and anoxic waters in the ocean. Existing oxygen minimum zones are projected to expand, with anoxia potentially impacting more than 20% of the area of the ocean. In all our scenarios, the expanding oxygen minimum zones impinge on the seafloor of continental margins. When accounting for a minimum of a 2-fold increase of nutrient inputs, at least 0.3 to 1.3% of the seafloor becomes covered by anoxic bottom waters.
Ocean deoxygenation is known to lead to major changes in ocean biogeochemistry. One major consequence of the nutrientdriven decline in ocean oxygen in our model is the increase in organic carbon burial on continental margins and in the deep sea. This increase is the combined result of increased production of organic matter in the surface ocean and enhanced preservation of organic matter. Our results suggest that, on a time scale of up to 50 kyrs, this burial of organic matter can significantly contribute to the removal of atmospheric carbon from fossil fuel burning and other human activities. Our results critically depend on increased nutrient inputs from CO 2 -driven enhanced weathering on land and from dust over longer periods of time. Such longterm changes in nutrient input, anoxic area and organic carbon burial are not unlike those that are thought to have occurred during periods of oceanic anoxia in Earth's past. This implies that modern human activities will likely impact marine life and ocean biogeochemistry beyond timescales of several hundred years and even create a new episode of long-term widespread ocean deoxygenation.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
IR, VP, and CS have made substantial intellectual contributions to this work and approved it for publication.

ACKNOWLEDGMENTS
We thank Christoph Heinze for substantial modeling assistance and Teresa Beaty for providing her modeling results of scenarios with reduced ventilation and increased CO 2 . The code used in this study is available upon request.