Ocean Phosphorus Inventory and Ocean Deoxygenation: Large Uncertainties in Future Projections on Millennial Timescales

Previous studies have suggested that weathering and benthic phosphorus (P) fluxes, triggered by 10 climate warming, can increase the oceanic P inventory on millennial time scales, promoting ocean productivity and deoxygenation. In this study, we assessed the major uncertainties in projected P inventories and their imprint on ocean deoxygenation using an Earth system model of intermediate complexity for a business-as-usual carbon dioxide (CO2) emission scenario until year 2300 and subsequent linear decline to zero emissions until year 3000. 15 Model results suggest a large spread in the simulated oceanic P inventory due to uncertainties in (1) assumptions for weathering parameters, (2) the representation of bathymetry on slopes and shelves in the model bathymetry, (3) the parametrization of benthic P fluxes and (4) the representation of sediment P inventories. Our best estimate for changes in the global ocean P inventory by the year 5000 caused by global warming amounts to +30% compared to pre-industrial levels. Weathering, benthic 20 and anthropogenic fluxes of P contributed +25%, +3% and +2% respectively. The total range of oceanic P inventory changes across all model simulations varied between +2% and +60%. Suboxic volumes were up to 5 times larger than in a model simulation with a constant oceanic P inventory. Considerably large amounts of the additional P left the ocean surface unused by phytoplankton via physical transport processes as preformed P. Nitrogen fixation was not able to adjust the oceanic 25 nitrogen inventory to the increasing P levels or to compensate for the nitrogen loss due to increased denitrification. This is in contrast to palaeo reconstructions of large-scale deoxygenation events. We suggest that uncertainties in P weathering, nitrogen fixation and benthic P feedbacks need to be reduced to achieve more reliable projections of oceanic deoxygenation on millennial timescales.

Phosphorus is considered the ultimate limiting nutrient for ocean productivity at the global scale (Tyrrell, 1999).Elevated supply of P to the ocean stimulates production and export of organic matter and deoxygenation, which possibly drives more intense oxygen depletion in the oxygen deficient zones 35 and along the continental margins, with release of additional P from sediments turning anoxic (Van Earth Syst.Dynam.Discuss., https://doi.org/10.5194/esd-2018-58Manuscript under review for journal Earth Syst.Dynam.Discussion started: 2 October 2018 c Author(s) 2018.CC BY 4.0 License.
Cappellen and Ingall, 1994;Palastanga et al., 2011).Such a positive feedback was discussed for a global warming scenario under present-day conditions (Niemeyer et al., 2017) as well as for large-scale deoxygenation events in the Cretaceous era, the so-called oceanic anoxic events (OAEs) (Tsandev and Slomp, 2009;Monteiro et al., 2012;Ruvalcaba Baroni et al., 2014).For the Cretaceous, it has been 40 suggested that atmospheric carbon dioxide (CO 2 ) concentrations as high as 1000 to 3000 ppmv, driven by enhanced CO 2 outgassing from volcanic activity (Jones and Jenkyns, 2001;Kidder and Worsley, 2012), have triggered OAEs (Damsté et al., 2008;Méhay et al., 2009;Bauer et al., 2016).The warmer climate during past OAEs increased weathering on land (Blättler et al., 2011;Pogge von Strandmann et al., 2013), leading to an enhanced supply of nutrients, in particular P, increasing the oceanic nutrient 45 inventory and driving the positive feedback mentioned above.Furthermore, the enhanced release of P from sediments were suggested to maintain high levels of productivity in the Cretaceous ocean (Mort et al. 2007;Kraal et al. 2010), which would contribute to the development of OAEs.Evidence in the palaeo record indicates that the Earth had experienced several climate OAE-like states, with large-scale anoxia, euxinia and mass extinctions (Kidder and Worsley, 2010).50 Could such OAEs also appear in the near future under contemporary global warming?High CO 2 concentrations in the atmosphere seem to be one driver for initiating OAEs and ocean deoxygenation.
Projected anthropogenic CO 2 emissions may lead to atmospheric CO 2 concentrations exceeding 1000 ppmv at the beginning of the 22 nd century if emissions continue to increase in a business-as-usual scenario (Meinshausen et al., 2011).Although anthropogenic CO 2 emissions occur over a short period 55 compared to the long-term and relatively constant volcanic CO 2 emissions during OAEs (Kidder and Worsley, 2012), elevated atmospheric CO 2 concentrations could persist for many millennia (Clark et al., 2016).This may provide the conditions for long-term climate change and large-scale deoxygenation.There is thus some concern that anthropogenic CO 2 emissions could potentially trigger another OAE (Watson et al., 2017).Yet, Kidder and Worsley (2012) argue that emissions of global 60 fossil fuel reserves are insufficient to drive a modern OAE, but may instead lead to widespread suboxia.
During climate warming, ocean productivity could switch from P to nitrogen (N) limitation (Saltzman, 2005).N limitation could arise from enhanced denitrification in a more anoxic ocean, but at the same time low N to P ratios would be expected to stimulate N 2 -fixation by diazotrophs (Kuypers et al., 65 2004).During warmer and wetter climates, however, a reduced supply of aeolian Fe may limit marine N 2 -fixation (Kidder and Worsley, 2010).N 2 -fixation in regional proximity with OMZs can lead to net N losses due to mass balance constraints (Landolfi et al., 2013), which may even reverse the net effect of N 2 -fixation on the nitrogen inventory.
Recently, Niemeyer et al. (2017) showed in a model study that P weathering and sedimentary P release 70 in a business-as-usual CO 2 -emission (RCP8.5)scenario could strongly enlarge the marine P inventory and lead to a 4 to 5-fold increase in the suboxic water volume (dissolved oxygen (O 2 ) concentrations less than 5 mmol m -3 ) on millennial timescales.Here, we build on this study and test the sensitivity of the marine P and O 2 inventories in a climate change scenario on millennial timescales to different model formulations of P weathering and benthic fluxes.We aim to provide better constraints on future 75 ocean deoxygenation and assess the biogeochemical feedbacks triggered by P addition.In Sect.

Model
We applied the University of Victoria (UVic) Earth System Model (ESM) version 2.9 (Weaver et al., 2001), which has been used in several studies to investigate ocean oxygen dynamics (Schmittner et al., 85 2007;Oschlies et al., 2008;Getzlaff et al., 2016;Keller et al., 2016;Landolfi et al., 2017).The UVic model consists of a terrestrial model based on TRIFFID and MOSES (Meissner et al., 2003), an atmospheric energy-moisture-balance model (Fanning and Weaver, 1996), a sea-ice model (Bitz and Lipscomb, 1999) and the general ocean circulation model MOM2 (Pacanowski, 1996).Horizontal resolution of all model components is 1.8° latitude x 3.6° longitude.The ocean model has 19 layers 90 with layer thicknesses ranging from 50 m for the surface layer to 500 m in the deep ocean.The marine ecosystem was represented by a NPZD model (Keller et al., 2012).Organic matter transformations (production, grazing, degradation) were parameterized using fixed stoichiometric molar ratios (C:N:P, 106:16:1) and directly related to the production and, in oxygenated waters, utilization of O 2 (O:P, 160).
When O 2 is depleted in the model, organic matter is respired using nitrate (NO 3 -) (i.e.microbial 95 denitrification).An O 2 concentration of 5 mmol m -3 was used as switching point from aerobic respiration to denitrification.Sedimentary denitrification was not considered in this model configuration, such that water column denitrification and N 2 -fixation dictate the oceanic N balance.No explicit iron cycle was simulated and iron limitation was approximated with prescribed seasonally varying dissolved iron concentrations (Keller et al., 2012).Parameterizations of benthic and weathering 100 fluxes of P were extended from the study of Niemeyer et al. (2017).Implementations of a calcium carbonate sediment model (Archer, 1996) and a parameterization for silicate and carbonate weathering (Meissner et al., 2012) were applied in all simulations.When P weathering and anthropogenic P fluxes were applied (see Sect. 2.2), the global P flux was distributed over all river basins, in every grid box, weighted by river discharge rates.105

Experimental Design
Twelve different model simulations were performed to explore the range of uncertainties for the longterm development of the oceanic P inventory (Table 1).Each simulation started from an Earth system state close to equilibrium under preindustrial atmospheric CO 2 concentrations, prescribed wind fields and present-day orbital forcing.Spin-up runs lasting 20,000 simulation years or longer were made for 110 each simulation to reach equilibrium.In the spin-up runs for simulations with benthic P fluxes (purple and red in Table 1), the marine P inventory was kept constant by instantaneously compensating oceanic Earth Syst.Dynam.Discuss., https://doi.org/10.5194/esd-2018-58Manuscript under review for journal Earth Syst.Dynam.Discussion started: 2 October 2018 c Author(s) 2018.CC BY 4.0 License.
P loss (burial) by P weathering fluxes to the ocean.For model simulations without benthic P fluxes (black and blue in Table 1), one common spin-up run was performed without P weathering fluxes.
All transient simulations started in the year 1765 and ended in year 5000.Simulations were forced with 115 anthropogenic CO 2 emissions (fossil fuel and land use change) according to the extended RCP 8.5 scenario until the year 2300 (Meinshausen et al., 2011), followed by a linear decline to zero CO 2 emissions by the year 3000.Warming from non-CO 2 greenhouse gases and the effect of sulphate aerosols were prescribed as radiative forcing (Eby et al., 2013).Non CO 2 -emission effects from landuse change were not considered.The reference simulation (Ref) was performed without weathering and 120 without burial fluxes of P, meaning that the P inventory of the ocean remained unchanged.The remaining transient simulations applied either variable climate-sensitive weathering anomalies (without burial) or time-variable burial fluxes (with constant weathering) to the ocean (Table 1).

Burial experiments
The water column model is not coupled to a prognostic and vertically resolved sediment model.125 Instead, P burial in the sediment (BUR P ) was determined in every grid box from the difference between the simulated detritus P rain rate to the sediment (RR P ) and the benthic release of dissolved inorganic P from the sediment (BEN P ): where RR P is the detritus flux from the ocean (in P units).BEN P was calculated locally by a "transfer function", which parameterizes sediment/water exchange of P as a function of the rain rate of organic 130 matter and the bottom water O 2 concentration.Preferential P release, relative to carbon (C), is observed in sediments overlain by O 2 -depleted bottom waters (Ingall and Jahnke, 1994).Benthic P release was dependent on the dissolved inorganic carbon release (BEN C ) from organic matter degradation in the sediment and the C:P regeneration ratio r C:P (Wallmann, 2010; equation 2): BEN C was computed (Eq.3a) as the difference of the carbon rain rate to the sediment (RR C ) and a 135 'virtual' organic carbon burial flux (BUR C ).There is no explicit treatment of organic C burial in the model, and instead all organic C is remineralized in the deepest ocean layer.BUR C is dependent on the simulated organic C rain rate and bathymetry (Flögel et al., 2011).Burial of organic C is more efficient on the shelf and continental margins (Eq.3b) than for the deep sea (Eq.3c, sediment below 1000m water depth): 140 where RR C is in mmol C m -2 a -1 .r C:P (in Eq. 4) depends on the bottom water oxygen concentration and was calculated according to (Wallmann, 2010;equation 4).
where O 2 is in mmol m -3 and the coefficients and their uncertainties are Y F =123±24; A=112±24; Under low O 2 conditions, r C:P is lower than 106, which leads to a preferential P release from organic 145 matter and, eventually, a net release of P from the sediment (BEN P > RR P , in Eq. 1).
The UVic model has a coarse standard model bathymetry with a horizontal resolution of 1.8° latitude x 3.6° longitude.This does not adequately represent continental shelves and slopes.To overcome this limitation, sinking organic matter interacts with the sediment on a detailed subgrid bathymetry (Somes et al., 2013).The subgrid bathymetry was inferred from ETOPO2v21 (National Geophysical Data 150 Center, 2006).Fractional coverage of every ocean grid box by seafloor was calculated on each model depth level.Sinking organic matter is partially intercepted at the bottom of each grid box by a sediment layer and the intercepted amount depends linearly on the fractional coverage of the grid box by seafloor.At the seafloor, organic matter is remineralized on the subgrid bathymetry for P in accordance with Eq. (1) and Eq. ( 2), whereby organic C and N are completely remineralized under oxygen or 155 nitrate utilization without any burial.The subgrid-scale parameterization leads to a better vertical representation of benthic fluxes of P and simulated sediment rain rates of detritus (see Sect. 3.2).The subgrid bathymetry does not affect other processes like circulation, advection or mixing.
Burial fluxes of P were applied in the simulations Bur, Bur_noSG, Bur_Dun, Bur_low, Bur_high and Bur_res.The default Bur model configuration uses Eq. (3) (Flögel et al., 2011) and the subgrid-scale 160 bathymetry.Uncertainties in benthic P fluxes were examined by modifying this default model configuration.For the simulation Bur_noSG (i.e.without subgrid-scale parameterization), P fluxes at the sediment-ocean interface were calculated using the coarser standard model bathymetry, which barely reproduce the global coverage of shelf areas (compare hypsometries in suppl.Fig. S1).
In the Bur_Dun simulation BUR C was calculated using Eq. ( 5) with RR C in mmol C m -2 d -1 ; Dunne et 165 al. (2007): Where c = 7 mmol C m -2 d -1 .This parameterization leads to high (low) organic C burial rates for high (low) organic C rain rates.This formulation is different to the standard formulation of burial in Eq. (3b, c) where burial depends on the C rain rates and in addition on the water depth.In the standard formulation, C burial is an magnitude larger in slope and shelf regions compared to the deep ocean.170 We examined the sensitivity of P burial to the uncertainty of the parameters in Eq. ( 4) describing the carbon to phosphorus regeneration ratio r C:P .Given means and standard deviations for the parameters Y F =123±24; A=112±24; r=32±19 and assuming a Gaussian distribution, 100,000 independent coefficient combinations were assembled to calculate offline a range of global P burial estimates.For the offline calculation, preindustrial fields of O 2 and RR C were extracted from the simulation Bur with 175 a temporal resolution fine enough to resolve seasonal variations in the data.Global P burial varied between 0.21 TmolP a -1 (Bur_low) and 0.60 TmolP a -1 (Bur_high) for a confidence interval of 90% (coefficients are shown in Table 1).Individual spin-ups were performed for the Bur_low and Bur_high simulation to check that the offline calculated P burial corresponded to the online values from the spinup.Only minor differences between the O 2 fields of the Bur spin-up and the spin-ups for Bur_low and 180 The implemented transfer functions (Eq. 2 and 4) assume unlimited local reservoirs of sedimentary P, meaning that the cumulative release of P may exceed the local inventory of P in the sediment if the benthic fluxes are sustained over a longer period of time.In the simulation Bur_res we tested the 185 impact of this simplification by applying sediment inventory restrictions to sediment P release.An upper limit sediment P inventory was calculated based on the following assumptions for the continental shelf and slope.We assume that the top 10 cm of the sediment column are mixed by organisms and are hence regarded as the active surface layer that is in contact with the overlying bottom water.
Considering a mean porosity of 0.8 and a mean density of dry particles of 2.5 g cm -3 , the mass of solids 190 in this layer is 5 g cm -2 (Burwicz et al., 2011).The mean concentration of total P in continental shelf and slope sediments is 0.07 wt-% equal to 22.6 µmol/g (Baturin, 2007).Together, these assumptions convert to a local inventory of total solid P in the active surface layer of 113 µmol cm -2 .We assume that shelf and slope sediments can release up to 100 % of this total inventory (RES P ) under low oxygen conditions.The reservoir can be fully replenished by P supply from the water column.Any excess P is 195 assumed to be permanently buried: Local values of RES P adjust during the spin-up according to the environmental conditions.Our pragmatic sediment inventory approach most likely overestimates the upper limit of P that can be released from the sediments.For example, under low O 2 conditions, part of the releasable or reactive P is transformed into authigenic P and permanently buried (Filippelli, 2001).200 All Bur experiments applied a constant global weathering flux (W P,const ) as established during the respective spin-up run (see Table 1 for values of W P,const for the different Bur experiments).

Weathering Experiments
Uncertainties in the ocean P inventory due to weathering processes and anthropogenic fluxes of P were examined with the model simulations Anthr, Weath0.05,Weath0.10,Weath0.15 and Weath0.38. 205 In simulations Weath0.05,Weath0.10,Weath0.15,Weath0.38 the global weathering flux of P to the ocean (W P ) was parameterized in terms of an anomaly relative to a preindustrial P weathering flux (W P,0 ) according to Eq. ( 8).
The weathering function f is given in Eq. ( 9).Values of W P,0 are given in Table 1 and derived below.
The chosen anomaly approach assumes that, at steady state, W P,0 is balanced by a respective global 210 burial flux and hence can be neglected during the spin-up.In these simulations no benthic P fluxes were applied and for preindustrial conditions the weathering function f(NPP,SAT) equals 1 and hence W P equals 0 TmolP a -1 .The dynamic weathering function f (Eq.9) was adopted from Niemeyer et al. (2017) and is originally based on an equation from Lenton and Britton (2006) for carbonate and silicate weathering.Following Niemeyer et al. (2017), we assumed that the release of P is proportional to the 215 chemical weathering of silicates and carbonates on a global scale.Equation ( 9 with NPP 0 and SAT 0 being the respective preindustrial values.Increasing SAT and NPP lead to enhanced weathering.The upper estimate of W P,0 in Weath0.38 was inferred from the P burial reference 220 simulation Bur, assuming that the global integral of burial is compensated by the preindustrial global weathering flux (i.e. the global marine P inventory is in steady state).With the simulations Weath0.05,Weath0.10,Weath0.15,Weath0.38 we explored the range of W P,0 estimates as derived from observational studies, which range from 0.05 to 0.30 TmolP a -1 (see Fig. 1, Benitez-Nelson, 2000;Compton et al., 2000;Ruttenberg, 2003).These studies indicate that total P fluxes to the oceans are 225 higher than interfered from dissolved inorganic P fluxes.A small amount of fluvial P is delivered to the ocean as dissolved inorganic P, but the majority (90%) is particulate (inorganic and organic) P (Compton et al., 2000).The fast transformations between dissolved and particulate P in rivers (seconds to hours) (Withers and Jarvie, 2008) suggest a much higher amount of P that is available for marine organism than derived from dissolved inorganic P concentrations.A large amount of bioavailable P in 230 rivers is present as loosely sorbed and iron-bound P. Estimates of bioavailable P are given in Fig. 1 (Benitez-Nelson, 2000; Compton et al., 2000;Ruttenberg, 2003), which are much higher than the estimates for dissolved inorganic P (0.018 TmolP a -1 from Seitzinger et al. (2005) or 0.03 TmolP a -1 from Filippelli (2002)).Taking into account only fluxes of dissolved inorganic P would strongly underestimate the effect of weathering fluxes as a P source to the ocean.The weathering 235 parametrization (Eq.9) was used to scale preindustrial fluvial fluxes of bioavailable P that is delivered in UVic to the ocean as dissolved inorganic P. In the model, no distinction was made between particular and dissolved fluvial fluxes of P.
Uncertainties to other weathering parameterizations were not investigated in this study.Our parameterization predicts similar weathering rates to other weathering formulations (Meissner et al.,240 2012, their Fig.6a).Since weathering is calculated on a global scale, we cannot study the effects of regional lithology and soil shielding on weathered P (Hartmann et al., 2014).
Finally, global anthropogenic P fluxes from fertilization, soil loss due to deforestation and sewage as projected by Filippelli (2008) were prescribed in the simulation Anthr.

Uncertainties in Phosphorus Inventory 245
The large range of projected global phosphorus (P) fluxes to the ocean from sediments or weathering (Fig. 2a) leads to uncertainties in future P inventories by up to 60% of the present-day value until year 5000 (Fig. 2b).All simulations show negligible differences in atmospheric CO 2 concentrations and hence undergo a similar climate development.Maximum CO 2 concentrations of 2200 ppmv were reached in year 2250 and then declined to 1100 ppmv by year 5000, comparable to results from Clark 250 et al. (2016).We found that the large range in P fluxes was not related to differences in the climate or atmospheric CO 2 forcing, but rather to differences in parameterizations of P land-ocean (Sect.3.1) and sediment-ocean (Sect.3.2) interactions.

Fluvial P Fluxes: Weathering and Anthropogenic
Largest uncertainties in the P inventory are related to the large range of P weathering fluxes (Fig. 2, 255 blue curves).Upper and lower estimates of P weathering fluxes differ by a factor of 6 (Fig. 2a, blue lines).In our weathering simulations, weathering anomalies depend linearly on the preindustrial weathering flux, W P,0 , estimate (see Eq. 8) because the climate development is essential equal across the simulations.Therefore, the choice of W P,0 (Fig. 1a) is a major source of uncertainty for projected future land-ocean P fluxes.260 Weathering fluxes increased from the pre-industrial value by a factor of 2.5 until year 5000 for atmospheric CO 2 concentrations of 1100 ppmv.This is comparable with the two-to four-fold increase in weathering fluxes estimated during OAE 2 approximately 91 Ma ago (Pogge von Strandmann et al., 2013) when atmospheric CO 2 concentrations increased to about 1000 ppmv (Damsté et al., 2008).
In contrast to weathering-induced P input, anthropogenic P fluxes (Filippelli, 2008) influence the 265 global marine P inventory only in the near future (Fig. 2a, black dashed line).A decline in anthropogenic P fluxes after year 2100 is expected due to the depletion of the easily reachable phosphorite mining reserves (Filippelli, 2008).

Sediment Fluxes: Parameterizations, Subgrid Bathymetry, Sediment Reservoir
The release of P from the sediment is strongly dependent on the O 2 concentration in the water above 270 the sediments (Wallmann 2003;Flögel et al. 2011).Climate warming reduces O 2 solubility and ventilation of the ocean, which decreases the global O 2 content (more details in Sect.4).The general decrease in ocean O 2 content may therefore cause preferential release of P from marine sediments.
Differences in sediment P fluxes in our simulations are related to uncertainties in the parameterization of the transfer function (Fig. 2, red lines, -0.01 to 0.22 TmolP a -1 by year 5000), to different 275 representations of the bathymetry (Fig. 2, purple dashed line, 0.06 (without subgrid) and 0.12 (Bur) TmolP a -1 ) and to the way sediment P reservoirs in the sediment are represented (Fig. 2, purple solid line, -0.01 (limited reservoir) and 0.12 (unlimited reservoir, Bur) TmolP a -1 ).
The global P burial of approximately 0.2 TmolP a -1 (Fig. 3) (Filippelli and Delaney, 1996;Benitez-Nelson, 2000;Ruttenberg, 2003) is relatively well reproduced by simulations Bur_low and Bur_Dun.280 The simulation with the standard UVic bathymetry (Bur_noSG) underestimates P burial by 60% while the simulations Bur_high, Bur and Bur_res overestimate P burial by 180%, 90% and 80% with respect to estimates based on observations.The transient response of the P release to O 2 was stronger for simulations with low burial and vice versa (Fig. 2), except for simulation Bur_res.In Bur_res, a significant reduction in the transient P release occurred due to the implementation of a finite P 285 reservoir, with net global P loss due to enhanced burial at the end of the simulation.In year 5000, global P concentrations increased in Bur_res by only 0.06 mmolP m -3 compared to the global mean preindustrial concentration of 2.17 mmolP m -3 .This is six-fold smaller than the increase of 0.36 mmolP m - 3 in simulation Bur with an assumed unlimited P reservoir.The small increase in the oceanic P inventory in Bur_res can be explained by the reduction in P sediment inventory rather than by changes 290 in the rain rate of particulate organic matter to the sediment (RR C ).In Bur, a rapid increase in the benthic P flux appeared in areas where the water turned suboxic and thus drove a positive benthic  2) and to other field data studies reporting a range from 900 to 2300 TgC a - 1 (Fig. 4) (Muller-Karger et al., 2005;Burdige, 2007;Dunne et al., 2007;Bohlen et al., 2012).
In summary, subgrid bathymetry leads to a substantial improvement of the representation of RR C to the 300 sediment.More realistic benthic fluxes of P could be also attained by adjusting parameters for r C:P (Eq.4) or by using the function of Dunne et al. (2007) to calculate BUR C (Eq. 5).The implementation of a finite P reservoir in the sediment has a substantial impact on the transient development of the global P inventory on millennial time scales.

Ocean Deoxygenation and Suboxia 305
Climate change influences ocean oxygen content by changes in circulation, ocean temperature and the degradation of organic matter.In warming surface waters, the solubility of O 2 decreases along with an increase in stratification, which together cause the deeper ocean to becomes less ventilated (Bopp et al., 2002;Oschlies et al., 2018).Changes in export production and the degradation of organic matter in the ocean interior also affects O 2 content.In the following, we analyze the impact of different ocean P 310 inventories on ocean deoxygenation and suboxia (Fig. 5).For a more detailed analysis we compare Weath0.15 to the Ref simulation.In the Weath0.15simulation, the assumed preindustrial weathering flux compares well to estimates from observations (Fig. 1).
In the Ref simulation, global suboxic volume increased due to climate change from 0.3 to 1% until year 5000 and the suboxic sediment area increased from 0.06 to 0.23% (Fig. 5, black line).In the Weath0.15315 simulation, the increase in suboxic volume (suboxic sediment area) was more than 2 (3) times higher than for the Ref simulation.The expansion of suboxic sediment areas was also enhanced for simulations with benthic fluxes, which could be related to regional feedbacks between increasing marine productivity, decreasing oxygen and enhanced sedimentary P release ( (Tsandev and Slomp, 2009).The explicitly simulated finite sedimentary P reservoir in simulation Bur_res places an upper 320 limit to the benthic release of P and dampens these regional feedbacks, resulting in a weaker spreading of suboxic waters by only 17% compared to the Ref simulation.
In the following sections, we show how the expansion of suboxia is related to net primary production in the ocean (NPP), the export of organic matter (Sect.4.1) and to nitrogen limitation (Sect.4.2).

Enhanced Biological Pump
The biological carbon pump can be summarized as the supply of biologically sequestered CO 2 to the 330 deep ocean.In the euphotic zone phytoplankton and diazotrophs take up CO 2 , a process that is intensified by elevated PO 4 concentrations in the surface ocean (Fig. 6a).Part of the organic matter sinks out of the euphotic zone (Fig. 6b) to the ocean interior, where it is respired using O 2 .It is therefore P supply to the surface waters that explains the differences in deoxygenation between the simulations.Circulation changes could also affect the supply of O 2 to the ocean interior.However, no 335 significant differences in climate and circulation appeared among the simulations and therefore the global-warming induced circulation changes affected all simulations in the same way.
In the Ref simulation, net primary production (NPP, Fig. 6a black line) increased from 45 to 70 TmolP a -1 (57 to 89 GtC a -1 ) by the end of the simulation.In Weath0.15,enhanced P supply to the ocean led to a doubling of NPP compared to the Ref simulation.The P inventory increased continuously, but NPP 340 did not follow this trend and instead peaked in year 4000.In year 5000, all simulations, excluding Weath0.38,showed a similar response of NPP to the P addition with an increase in NPP of 19 TmolP a - 1 (relative to the Ref simulation) per 10% increase in P inventory.In Weath0.38 the response was weaker and NPP increased by 8 TmolP a -1 per 10% rise in the P concentration.P is less effectively utilized in simulations with large oceanic P inventories.Higher ocean temperatures enhanced 345 remineralization of organic matter in the shallower ocean so that the overall export to NPP ratio decreased from its preindustrial value of 0.12 to an average value among all simulations of 0.08 by year 5000.To summarize, NPP and export of organic matter is sensitive to P addition.However, the proposed positive feedback between P, NPP, export of organic matter, and deoxygenation was limited in our simulations due to a negative feedback related to nitrate availability.This is showed and 350 explored in the following section.

Nitrogen Limitation
At the end of the spin-up the N sink by denitrification and the N source by N 2 -fixation were balanced.
In the Ref simulation, climate warming enlarged the oxygen minimum zones, which enhanced denitrification in the tropics (not shown).In all simulations, N 2 -fixation was stimulated by the addition 355 of P to the ocean and was sensitive to rapid changes in the supply of P (compare Fig. 7a and Fig 2a).However, N 2 -fixation by diazotrophs (Fig. 7a) was not able to balance the loss by denitrification because low temperatures in polar regions and iron limitation at lower latitudes inhibit growth of diazotrophs (Fig. S2), leading to a substantial amount of excess phosphate in the surface waters of these regions (Fig. S3).As a consequence, nitrate decreased globally by 4 mmolN m -3 until year 5000 360 (Fig. 7b).In low-N and high-P environments, diazotrophs have a competitive advantage over other phytoplankton.In the simulations with P supply, N 2 -fixation by diazotrophs was stimulated (Fig. 7a) and partly counteracted the nitrate loss by denitrification.The loss in nitrate led to a decrease in globally averaged N to P ratios.In the Ref simulation, N:P decreased from 14 to 12 and for the Weath0.15simulation it decreased to 10, which contributed further to a N limiting ocean.The nitrogen 365 cycle was not able to recover from the decrease in N:P ratio with respect to pre-industrial values.In the phytoplankton to adapt to changing N:P ratios was not considered.Phytoplankton organisms may utilize the excess in P if they are able to adapt their stoichiometry to the decreasing N:P ratios, which could then lead to an increase marine biological productivity and greater deoxygenation.370

Temporal Variations of Deoxygenation
Anomalies in circulation, ocean temperature and remineralisation of organic matter affect oceanic O 2 levels in a climate-warming scenario.In the Ref simulation, the O 2 inventory (Fig. 8a) decreased by 60 Pmol O 2 by the year 3000 and then recovered by year 5000.In Weath0.15,weathered P enhanced deoxygenation and led to a greater decrease in O 2 than in the Ref simulation.The O 2 decrease was up to 375 70 Pmol in year 3300 and O 2 still showed a negative anomaly of 24 Pmol O 2 by the year 5000.Global anomalies in oxygen were due to changes of the Apparent Oxygen Utilization (AOU, Fig. 8b) and the O 2 saturation level (Fig. 8c).Changes in O 2 saturation were similar across the model simulations and followed with a delay surface ocean temperature.The circulation and ventilation of the ocean were similar in the model simulations because differences in surface temperatures were negligible and the 380 atmospheric forcing of the ocean circulation was identical.So that differences in AOU depended almost only on biological O 2 consumption and AOU anomalies were directly yet inversely related to the changes in O 2 levels.Hence, biological consumption explained variations in O 2 content among the different model simulations (compare Fig. 8a and 8b).Increasing O 2 utilization contributed to the decrease of the O 2 until year 3000.Thereafter, a distinct negative trend in AOU with a similar slope 385 was observed among all simulations and contributed to a re-oxygenation of the ocean.For simulations with larger P inventories, the AOU had a larger positive offset to the Ref simulation.
In a model with constant stoichiometry for elemental exchange by biological processes, anomalies in AOU (Fig. 9, blue lines) can be explained by the difference between total integrated nutrients (Fig. 9, red and black solid lines as anomalies) and preformed nutrients (Fig. 9, red and black dashed lines as 390 anomalies).Preformed nutrients correspond to the fraction that leaves the surface ocean unutilized by phytoplankton.For example in the Southern Ocean, a large fraction of nutrients leaves the surface as preformed nutrients.The fraction of utilized and preformed nutrients can change during a transient simulation and could affect the oxygen state of the ocean.
In the Ref simulation (Fig. 9a), the anomaly of preformed dissolved inorganic P was directly inverse to 395 the anomaly of AOU because the oceanic P Inventory was conserved in this simulation.Until year 2200, changes in circulation and climate are likely the main cause for the reduction in preformed N and 11 Sv in year 2200.The continuous warming and stratification of the ocean reduces the supply of nutrients to the surface layer from the deep ocean.This is consistent with a reduction of the export of organic matter until year 2200 (Fig. 6b).The balance between exported P out of the surface ocean and supplied P controls changes in AOU.We speculate that a weaker overturning increased the residence time of water and nutrients in the surface ocean.Nutrients staying longer in the euphotic zone are with Earth Syst.Dynam.Discuss., https://doi.org/10.5194/esd-2018-58Manuscript under review for journal Earth Syst.Dynam.Discussion started: 2 October 2018 c Author(s) 2018.CC BY 4.0 License.a higher probability biologically consumed.This implies more efficient utilization of nutrients and, hence, the reduction in preformed nutrients and an increase in AOU.
Enhanced suboxia after year 2200 drove excess denitrification and a decline in nitrate (Fig. 9a red solid line) in the Ref simulation.The decline in nitrate could explain the negative trend in AOU anomalies (Fig. 9a blue solid line) and therefore a negative feedback on the global deoxygenation.In year 2200, 410 overturning had started to recover quickly and increased to 21 Sv in year 3000 (+24% relative to preindustrial values), which drove a faster overturning of organic matter in the surface ocean and a decrease in global AOU.We assume that the slight increase in export by 5% (relative to preindustrial values) was not strong enough to compensate for the by +24% faster overturning, which reduced the residence time of nutrients in the surface ocean.415 P addition in the Weath0.15simulation stimulated N 2 -fixation by diazotrophs and counteracted N-loss by denitrification (Fig. 9b, red solid line).This led to an increase in N inventory by 17 Pmol O 2equivalents compared to the Ref simulation.Furthermore, the high availability of P seems to reduce preformed N by 6 Pmol O 2 equivalents.Both explain the difference in AOU between Weath0.15 and Ref of 24 Pmol O 2 at the end of the simulation (Fig. 8b).However, denitrification still exceeded N 2 -420 fixation, which led to low levels of nitrate.From year 5000 approximately all of the added P in the Weath0.15simulation remained unused by phytoplankton, left at the surface ocean as preformed P and was afterwards stored in the deep ocean.Phytoplankton was not able to utilize the excess P because it was limited in nitrate.Diazotrophs were not able to compensate for the lack in N due to iron limitation and low surface temperatures in the polar oceans.The denitrification feedback driven by the spread of 425 suboxic conditions in the tropics had reduced further N availability for the phytoplankton and reduced the effect of P addition on the global oxygen level.

Discussion and Conclusions
The P inventory is very sensitive to the weathering and benthic flux parameterizations tested in our model.Large uncertainties (Fig. 2, blue lines) derive from poorly constrained estimate for the 430 preindustrial P weathering flux that ranges from 0.05 to 0.30 Tmol P a -1 (Benitez-Nelson, 2000; Compton et al., 2000;Ruttenberg, 2003).The preindustrial weathering flux in simulation Weath0.15(0.15 Tmol P a -1 ) is well in this range.In this simulation, enhanced weathering leads to an increase in the global ocean P inventory by 25% until year 5000 (Fig. 2, blue dotted line).Benthic fluxes of P were simulated using transfer functions on a subgrid bathymetry.Applying the transfer functions without 435 taking into account the local sedimentary P inventory can greatly overestimate the release of benthic P on long time scales.In the UVic model, the application of finite benthic P inventories limited the benthic release significantly.Under low-oxygen conditions, sediments were P depleted already after a few years to decades.In our simulation, this resulted in an increase in the global oceanic P inventory by just 3% (Fig. 2, magenta solid line).This implies that benthic release of P is actually negligible in 440 comparison to the weathering fluxes of P, but the UVic model does not resolve coastal processes such as the deposition of reactive particulate P from rivers on the continental shelves and its dissolution and release to the water column.A more detailed representation of coastal processes would be necessary to simulate burial and release of fluvial P. Further, we find that a more realistic bathymetry substantially improves the simulated rain rate of particular organic carbon to the sediment (Table 2), particularly on 445 the shelf, which most models do not resolve.Anthropogenic P fluxes increased the global P inventory by just 2% (Fig. 2, black dashed line).In summary, our best estimate for changes in the total global ocean P inventory by the year 5000 amounts to +30%, which was dominated by weathering.This seems to be surprisingly high, but several studies indicate that changes in past climate could also have been accompanied with substantial changes in the P inventory but at a much lower pace (Planavsky et 450 al., 2010;Monteiro et al., 2012;Wallmann, 2014).
The increased P inventory (Fig. 2b) promotes deoxygenation (Fig. 5) and expansion of suboxia, but it also causes a net loss of nitrate, which appears to further limit the full utilization of P by phytoplankton in our simulations.Wallmann (2003), using a box model, already recognized that for a eutrophic ocean, nitrate might ultimately limit marine productivity.As a consequence, large amounts of P leave the 455 surface ocean as preformed P (Fig. 9b) with no further impact on O 2 levels in the ocean interior.Low N/P ratios are thought to give N 2 -fixers a competitive advantage over ordinary phytoplankton and lead to an increase in N 2 -fixation (Fig. 7a).A substantial increase in N 2 -fixation was also inferred from a palaeo study of sedimentary nitrogen isotopes by Kuypers et al. (2004).However, high denitrification rates remove nitrate from the global ocean and in the UVic model N 2 -fixers are not able to compensate 460 for this loss (Fig. 7b) because low temperatures in polar regions and iron limitation at lower latitudes inhibit growth of diazotrophs (Fig. S2) and a substantial amount of excess phosphate remains in the surface waters in these regions (Fig. S3).General circulation models without a N cycle, or box models without realistic representation of habitats suitable for N 2 -fixers, would miss this important negative feedback limiting global deoxygenation.465 Some additional model limitations are a cause for uncertainty in our results.We considered a fixed Redfield-ratio stoichiometry.In future deoxygenation studies, an optimality-based model for nutrient uptake with variable nutrient ratios (Pahlow et al., 2013) could be applied to investigate how well marine organisms adapt to a changing nutrient availability in the global ocean.Sea level change and the implied bathymetry change were not simulated in the UVic model.In future projections, higher surface 470 air temperatures would lead to a rise in sea level, which increase global coverage of shelf areas.Burial of P is more effective on the shelf (Flögel et al., 2011), which would remove P from the ocean and lead to a lower marine P residence time (Bjerrum et al., 2006).
To conclude, climate warming leads to a larger oceanic P inventory mainly due to addition of P by weathering, but also due to the release of P from the sediment and due to anthropogenic fluxes.A 475 realistic representation of shelf bathymetry improves the predicted benthic P fluxes.Transfer functions for benthic P release should consider the sedimentary P inventory.However, the largest uncertainties in the projection of oceanic P inventory are due to poorly constrained weathering fluxes of P.Although additional deoxygenation is driven by P addition to the ocean, the degree of deoxygenation -and hence the positive redox-related feedback on benthic P fluxes is eventually limited by the availability of N 480 and the apparent inability of the modelled N 2 fixation to respond to the larger P inventory.

Tables: 700
Table 1: Overview of simulations.P fluxes are given in TmolP a -1 .In the P weathering simulations, only weathering anomalies were applied.In the P burial simulations, a constant P weathering flux (W P,0 ) balances P burial (BUR P ) during the spin-up simulations.The preindustrial P inventory is identical in all simulations.More detailed information can be found in the text.Time in yrs 2 we Earth Syst.Dynam.Discuss., https://doi.org/10.5194/esd-2018-58Manuscript under review for journal Earth Syst.Dynam.Discussion started: 2 October 2018 c Author(s) 2018.CC BY 4.0 License.present the experimental design and the model parameterizations of continental P weathering and of benthic P release.In Sect. 3 we assess uncertainties in P fluxes due to different assumptions about the P weathering fluxes, different model formulations of benthic P fluxes, improved representation of bathymetry and anthropogenic P fluxes.Consequences for deoxygenation and for the biogeochemical 80 cycling of nutrients are discussed.
) describes the sensitivity Earth Syst.Dynam.Discuss., https://doi.org/10.5194/esd-2018-58Manuscript under review for journal Earth Syst.Dynam.Discussion started: 2 October 2018 c Author(s) 2018.CC BY 4.0 License. of terrestrial weathering to the change of global terrestrial net primary production (NPP) and global mean surface air temperature (SAT): Earth Syst.Dynam.Discuss., https://doi.org/10.5194/esd-2018-58Manuscript under review for journal Earth Syst.Dynam.Discussion started: 2 October 2018 c Author(s) 2018.CC BY 4.0 License.feedback between P release, productivity and deoxygenation.A limited supply of P from the sediment (Bur_Res) dampens this feedback.With a more realistic subgrid bathymetry, simulated pre-industrial RR C increased significantly from 295 180 to 1040 TgC a -1 on the shelf and globally from 900 to 1500 TgC a -1 compared to simulations without subgrid bathymetry.Pre-industrial RR C with subgrid bathymetry agrees better to estimates by Bohlen et al. (2012) (Table

Finally
, we show how changes in O 2 solubility and utilization vary over time and affect the global O 2 325 inventory (Sect.4.3).The latter approach gives another perspective because changes in O 2 inventories are a global integrated signal in comparison to the extent of suboxia, which are consequence of more local processes.Earth Syst.Dynam.Discuss., https://doi.org/10.5194/esd-2018-58Manuscript under review for journal Earth Syst.Dynam.Discussion started: 2 October 2018 c Author(s) 2018.CC BY 4.0 License.

P
in the Ref simulation since global N and P inventories were constant (Fig 9a, solid red and black line).During continuous and intense ocean warming, a weakening of the meridional overturning (not shown) reduced ocean ventilation.The meridional overturning decreased from 17 Sv (pre-industrial) to 400 Earth Syst.Dynam.Discuss., https://doi.org/10.5194/esd-2018-58Manuscript under review for journal Earth Syst.Dynam.Discussion started: 2 October 2018 c Author(s) 2018.CC BY 4.0 License.

Table 2 : Rain rate of particulate organic carbon (RR C ) to the seafloor for the shelf, slope and deep-sea areas from the observational estimate by Bohlen et al. (2012) and for the UVic model simulation Bur with and 710 without subgrid bathymetry. Preindustrial RR C shows no significant differences among all model simulations (expect for simulation Bur_noSG).
Earth Syst.Dynam.Discuss., https://doi.org/10.5194/esd-2018-58Manuscript under review for journal Earth Syst.Dynam.Discussion started: 2 October 2018 c Author(s) 2018.CC BY 4.0 License.