Swept under the carpet : organic matter burial decreases global ocean biogeochemical model sensitivity to remineralization length scale

Although of substantial importance for marine tracer distributions and eventually global carbon, oxygen, and nitrogen fluxes, the interaction between sinking and remineralization of organic matter, benthic fluxes and burial is not always represented consistently in global biogeochemical models. We here aim to investigate the relationships between these processes with a suite of global biogeochemical models, each simulated over millennia, and compared against observed distributions of pelagic tracers and benthic and pelagic fluxes. We concentrate on the representation of sediment–water interactions in common numerical models, and investigate their potential impact on simulated global sediment–water fluxes and nutrient and oxygen distributions. We find that model configurations with benthic burial simulate global oxygen well over a wide range of possible sinking flux parameterizations, making the model more robust with regard to uncertainties about the remineralization length scale. On a global scale, burial mostly affects oxygen in the mesoto bathypelagic zone. While all model types show an almost identical fit to observed pelagic particle flux, and the same sensitivity to particle sinking speed, comparison to observational estimates of benthic fluxes reveals a more complex pattern, but definite interpretation is not straightforward because of heterogeneous data distribution and methodology. Still, evaluating model results against observed pelagic and benthic fluxes of organic matter can complement model assessments based on more traditional tracers such as nutrients or oxygen. Based on a combined metric of dissolved tracers and biogeochemical fluxes, we here identify two model descriptions of burial as suitable candidates for further experiments and eventual model refinements.


Introduction
The relationship between organic matter degradation and oxygen consumption in the water column is of importance for the distribution of biogeochemical tracers in the global ocean and for features like oxygen minimum zones (OMZs), which, although covering only a small volume of the global ocean, are crucial for controlling the marine nitrogen inventory.OMZs are notoriously difficult to reproduce in numerical models, and their representation seems particularly sensitive to changes in the parameterizations of remineralization processes (Oschlies et al., 2008;Bianchi et al., 2013).In global models such as those used to predict the evolution of OMZ under global warming (e.g., Duteil and Oschlies, 2011;Stramma et al., 2012), remineralization at low oxygen levels has been parameterized in many different ways.
Reduction or cessation of remineralization in the absence of oxygen or other oxidants seems to be a plausible choice for models that focus on the open ocean.However, in such models sinking detritus may accumulate below suboxic zones and, according to our own sensitivity experiments, reach unrealistically high concentrations.This suggests that there is (i) too intense an export production or (ii) too little lateral mixing of the remineralization signal (or organic substrates) with surrounding water masses, a problem sometimes referred to as nutrient trapping.
Spatial reorganization of zooplankton around suboxic zones (e.g., Wishner et al., 2013), which is usually not considered in current models, may mediate the adverse effects of low oxygen concentrations, and prevent the otherwise too intense an accumulation of detritus at the seafloor.
Published by Copernicus Publications on behalf of the European Geosciences Union.
Another explanation for too intense an accumulation of remineralization products in water columns associated with eastern tropical upwelling regions could be the neglect of burial of particulate organic matter in the sediment.
Despite the fact that global burial of organic matter in the sediments is rather small compared to biogeochemical fluxes in the water column, and has often been neglected more specifically: treated as absent in global biogeochemical ocean models (e.g., Matear and Hirst, 2003), it nevertheless influences biogeochemical processes locally and, on long timescales, will impact tracer distributions in the entire water column and, on even longer timescales, may affect the oxygen content of the atmosphere (Najjar et al., 2007).
Inspection of various published biogeochemical models reveals that burial and also remineralization in low-oxygen waters are often treated in very different ways in different models.For example, in the "PISCES" model (as described in the Appendix of Aumont and Bopp, 2006), which includes both aerobic as well as anaerobic remineralization (the former depending on oxygen, the latter on both oxygen and nitrate), remineralization of organic matter ceases in the absence of oxygen.In this model sinking organic matter is ultimately buried below the sea floor, thereby preventing the accumulation of detritus at the bottom.However, current models that continue remineralization even under sub-or anoxic conditions do not seem to consider any burial, but remineralize the organic matter hitting the sea floor instantaneously (Marchal et al., 1998;Matear and Hirst, 2003;Najjar et al., 2007;Yool et al., 2011) or remineralize detritus in the last box with the normal remineralization rate of detritus in the water column (Kriest et al., 2010).In order to account for more realistic nutrient regeneration at the seafloor, Yool et al. (2013) added 2-D pools of organic and biogenic material at the seafloor to model "MEDUSA-1.0"(Yool et al., 2011).Overall, most formulations of water-sediment interactions in global biogeochemical circulation models have been introduced in an ad hoc manner, and a systematic evaluation of the impact of such formulations on the model's performance is lacking.
In this paper we investigate the impact of burial against the background of a commonly used parameterization of organic matter degradation in the water column.In particular, we examine the following questions: what is the effect of different descriptions of organic matter burial on simulated oxygen and phosphate distributions in the global ocean?How do the simulated exchanges across the sediment surface, which are implicit in the different parameterizations, compare to observational estimates?

Model
For our model experiments we use the framework of the Transport Matrix Method (Khatiwala et al., 2005;Khatiwala, 2007), coupled to the biogeochemical model NPZD-DOP described in Kriest et al. (2012, hereafter referred to as KKO12).This model is a phosphorus-based, five-component model, that simulates nutrients (as phosphate), phytoplankton, zooplankton, detritus and dissolved organic phosphorus (DOP).Detritus sinks with a speed that increases linearly with depth.Experiment "ref" applies a vertical increase in sinking speed of 0.058 (m d −1 ) m −1 .With a detrital remineralization rate of 0.05 d −1 this would, in equilibrium and in the absence of advection, result in a "Martin" flux profile (Martin et al., 1987, i.e., F ∝ z −b ) with an exponent b of 0.858 (experiment "ref" of Kriest et al., 2012).
Detritus that enters the last box remains as detritus and remineralizes with its given remineralization rate of 0.05 d −1 .Similar to Marchal et al. (1998), Matear and Hirst (2003) and Najjar et al. (2007), KKO12 assumed that for oxygen concentrations below 4 µmol O 2 m −3 remineralization of organic matter continues, but does not use any oxygen, thereby mimicking the consumption of other, nonspecified oxidants such as nitrate.
A new model feature investigated here is burial of detritus at the sea floor.In order to satisfy global mass conservation of phosphorus, buried detrital phosphorus is resupplied as phosphate via river runoff, while the non-buried detritus is resuspended in the water column.In the following we describe this model modification in more detail.

Parameterization of benthic burial
Let 0 ≤ f B ≤ 1 be the fraction of organic matter reaching the sea floor which is buried permanently in the sediment.For f B = 0 (no burial), detritus remains in the last box, where it slowly remineralizes, further depleting the oxygen (e.g., in Schmittner et al., 2008;Kriest et al., 2010).In the opposite case, f B = 1 (e.g., in the PISCES model; Aumont and Bopp, 2006) implies complete burial of organic matter reaching the sea floor.Between these two extremes f B may vary, among other things, with flux to the sea floor, water depth, or bottom water conditions.For example, Burdige (2007) suggested relating the so-called burial efficiency, f B , to the amount of organic matter falling onto the sea floor (F R ; here in mmol P m −2 d −1 ) via according to which the amount of burial in the sediment, F B , can be described via (2) In order to obtain an observational estimate of the parameters α (in (mmol P m −2 d −1 ) 1−β ) and β (dimensionless) of Eq. ( 2), we first compiled a data set that contains observations for burial, F B , and flux onto the sediment F R (sometimes also called "rain rate"; note that this usually refers to the sum of sedimentary remineralization and burial).Details about the data set can be found in Appendix A. The regression of the log-transformed data of burial (F B ) vs. rain rate (F R ) gives α = 1.6828 (mmol P m −2 d −1 ) 1−β and β = 1.799, with r 2 = 0.80. Figure 1 compares this relationship to other estimates of burial vs. rain rate.Flögel et al. (2011) divided the ocean into two domains, namely the coastal and the open ocean, for which they derived two separate empirical relationships between F B and F R .The two descriptions by Flögel et al. (2011) have very similar exponents β, but very different offsets α (see Fig. 1).Restricting our regression to regions with F B < 250 mmol C m −2 yr −1 , i.e., placing more emphasis on the open ocean, results in a much weaker dependence of burial on rain rate (α = 0.0176 (mmol P m −2 d −1 ) 1−β ), β = 1.022, r 2 = 0.61) than the regression for the full domain, and is comparable to the open ocean estimate by Flögel et al. (2011) with α = 0.024 (mmol P m −2 d −1 ) 1−β (converted with a C : P ratio of 106) and β = 1.05.
Our compilation over the full range of rain rates is, at the upper and lower end of range, similar to the functionalempirical relationship given by Dunne et al. (2007, see also our Fig.1), which relates the burial to a sigmoidal function of rain rate: . At low rain rates Dunne's data set agrees quite well with the openocean estimate by Flögel et al. (2011) and with the regression derived from our of our "open ocean" data compilation.At higher rain rates, the estimate by Dunne agrees with the shelf estimate by Flögel et al. (2011).Few observations are available in the mid-range of the rain rates between ≈ 0.01-0.1 mmol P m −2 d −1 , where the function by Dunne deviates from the other parameterizations.
For our model experiments we use the regression over the full data compilation, restricted to a maximum value of f B = 1.We further test a regression over the lower end of flux rates (i.e., for F B < 250 mmol C m −2 yr −1 ; see also Table 1).We also present some results of experiments carried out with Dunne's parameterization of burial.

Parameterization of the benthic nepheloid layer
For incomplete burial (f B < 1) the model also has to describe what happens to the organic matter that is not buried.Here we can distinguish two cases: (i) under oxic conditions organic matter can be instantaneously remineralized, implying a very "active" sediment, e.g., via bottom fauna, which instantaneously consumes and remineralizes any input of organic matter.(ii) Organic matter can be resuspended back into the overlying water where it is subject to mixing and advection with the bottom water.Observational evidence suggests that the deposition of organic matter onto the sea floor is not necessarily the ultimate sink of these particles, but that the (freshly) deposited organic particles may be easily remobilized, and become again subject to horizontal or even vertical (by mixing and upwelling) transport (e.g., Bacon and van der Loeff, 1989).2007) is denoted as green line.Thin black lines denote constant burial efficiencies of 100 (straight), 10 and 1 % (both dotted).Observations (symbols) and lines for algorithms by Dunne et al. (2007) and Flögel et al. (2011) were given in carbon units, and have been converted to phosphorus using a constant C : P ratio of 106.
We test the effect of instantaneous remineralization case (i) in a simulation where we assume instantaneous remineralization of all organic matter that hits the seafloor (experiment "INST", see text below and Table 1).As shown below, we find little difference with respect to a simulation where detritus remains suspended in the deepest grid box of the water column.For all burial scenarios, we therefore follow option (ii) and assume resuspension of the non-buried detritus, thereby mimicking a benthic nepheloid layer.

Phosphorus budget closure via river runoff
The above parameterization of permanent burial of organic matter in marine sediments removes P from the oceanic pelagic system.To account for mass preservation on long timescales, we re-supply it again via river runoff, using a subset of river runoff data of Table 2 of Perry et al. (1996), i.e., the volumetric flow rates of the largest rivers in the world.We did not use direct P-or N-discharge rates, because these may exhibit a strong anthropogenic influence, e.g., due to fertilizers, and therefore may not correspond to the steady-state system envisaged here.
Instead, we calculated the volumetric runoff of each river as a fraction of global runoff, and distribute it vertically over www.biogeosciences.net/10/8401/2013/Biogeosciences, 10, 8401-8422, 2013 the entire water column of the model boxes that receive river runoff, i.e., in the vicinity of a river mouth (thereby mimicking unspecified processes, such as frontal dynamics).To combine P loss due to burial and P gain due to river runoff, during runtime the loss of P due to burial in the sediment is integrated over total model area and -for computational efficiency -over one year.The resulting global annual loss is then resupplied continuously over the following year.We note that although the immediate resupply of buried phosphorus via river runoff may appear as a very fast "rock cycle", under equilibrium conditions addressed here, phosphorus sinks and sources have to be of identical size and response timescales become irrelevant.More details about the river runoff fluxes used can be found in Appendix B. Kriest et al. (2012) used circulation fields derived from a relatively coarsely resolved (2.8 • × 2.8 • , 15 vertical levels) circulation model (here referred to as MIT2.8) to examine the effect of different biogeochemical parameterizations on the spatial distribution of phosphate and oxygen.Because ventilation pathways and coastal processes may be important for the formation and persistence of OMZs, we here focus on biogeochemical models simulated with transport derived from the Estimating the Circulation and Climate of the Ocean (ECCO) project, which provides circulation fields that yield a best fit to hydrographic and remote sensing observations over the 10-year period 1992 through 2001 on a spatial resolution of 1 • × 1 • horizontal resolution with 23 vertical levels (Stammer et al., 2004).

Model experiments
We first examine the effect of the finer resolution and different circulation inherent in the ECCO matrix via comparison to the results obtained with the "coarser" MIT2.8 matrix.For both configurations, ECCO and MIT2.8, we have carried out a set of seven numerical experiments, where we varied the flux exponent from 0.429 up to 1.287 in model NPZD-DOP of KK012 (hereafter named CTL; see also Table 1).In these initial model simulations, burial is not included and resuspension of any detritus hitting the bottom is assumed, i.e., we apply option (ii) of Sect.2.2.
We then evaluate different biogeochemical models in the ECCO configuration.Starting from the model CTL, in setup BUR we introduce burial via Eq.(2), with our best fit to our sediment data compilation, α = 1.6828 (mmol P m −2 d −1 ) 1−β and β = 1.799.We further test the effect of a weak relationship between rain rate and burial using the open-ocean composite described above (α = 0.0176 (mmol P m −2 d −1 ) 1−β , β = 1.022).This setup is named WBUR.We finally carried out experiments with the algorithm suggested by Dunne et al. (2007), here denoted by DUNNE.To examine the effect of instantaneous benthic remineralization, as described above we also present experiments where all organic matter hitting the sea floor is immediately remineralized ("INST").
As in KKO12, for each of the setups we have carried out a set of four experiments where we varied the flux exponent b from the reference ("ref") experiment (b = 0.858) upwards to 1.0725 (experiment "s2" of KKO12; hereafter named experiment "slow") and 1.287 ("s1" of KKO12; hereafter named "very slow") and downwards to 0.644 ("s3" of KKO12; hereafter named "fast") and 0.429 ("s4" of KKO12; hereafter named "very fast"; see also Table 1).Note that "fast" corresponds to particle sinking speed, or a deep penetration of flux, whereas "slow" indicates slowly sinking particles, or a shallow remineralization of organic matter.As im KKO12, each of the model experiments was spun up for 3000 yr, using eight time steps per day for ocean tracer transport, and another eight biogeochemical time steps within each ocean time step.Annual mean concentrations and fluxes of year 3001 were then used for model evaluation and analysis (see also Kriest et al., 2010, for transient of similar model experiments).

Model assessment
Results from model experiments are compared to the gridded (1 • × 1 • ) analyzed set of observations of phosphate and oxygen compiled by Garcia et al. (2006a) and Garcia et al. (2006b, in the following referred to as "WOA05").We mapped the data set to the respective model resolution.By doing so, we release the coarser model from the penalties imposed by its lack of resolution, and investigate just the errors due to its numerical diffusion (e.g., due to upstream schemes in particle flux; see Kriest and Oschlies, 2011), circulation and biogeochemistry.To assess the models' performance with respect to dissolved tracers, we apply a global misfit function based on the volume-weighted root mean square error, as described in KKO12.Data sets used for comparison to simulated pelagic and benthic fluxes will be described in the respective section.

Impact of spatial resolution and physics
The volumetric distribution of global phosphate is quite similar in the two control model configurations (Fig. 2, upper panels).The overall mismatch of ECCO with respect to phosphate is lower than that for MIT2.8, but the sensitivity to changes in the remineralization length scale is very similar between the two different models driven by different circulation fields (Fig. 2b).The same increase in sinking speed (e.g., "ref" to "fast") has a larger effect in MIT2.8 than in ECCO.In the former, this increase results in a much more dispersive distribution of phosphate, with an overestimate of water volume with rather low and rather high concentrations.This effect is not so pronounced in ECCO (Fig. 2a, c).
Greater differences appear when comparing the oxygen distributions (Fig. 2, lower panels).First, the remineralization length scale shows a strong impact on the simulated oxygen distribution, with a strong volumetric overestimate of low oxygen regions when sinking is fast.Second, in MIT2.8 we find a bimodal volumetric distribution of oxygen, whereas the oxygen distribution in ECCO is more or less unimodal and, in this respect, in better agreement with WOA05.Particularly the fast sinking scenarios of MIT2.8 show oxygen concentrations below 50 mmol O 2 m −3 in large parts of the deep ocean, causing a peak of ocean volume with rather low concentrations.Reasons for differences among the two model configurations may be differences in the circulation fields, but also the differences in the vertical resolution, causing a higher numerical diffusion in MIT2.8 (see also Kriest and Oschlies, 2011).
For the biogeochemical model variations considered, ECCO fits observed oxygen much better than MIT2.8, while its results show the same sensitivity to changes in the remineralization length scale.In the following sections we will focus on results from simulations of the ECCO configuration.

Impact of burial and benthic remineralization on pelagic tracers
Similar to results obtained by Kwon and Primeau (2006) and Kriest et al. (2012), simulated nutrient concentrations in older waters such as in the mesopelagic and bathypelagic North Pacific are very sensitive to changes in the remineralization length scale in our models CTL and INST.Further, Kwon and Primeau (2006) found that nutrients in the eastern equatorial Pacific (hereafter named EEP) are particularly sensitive to changes in organic matter production and decay terms.As we will show below, in our model experiments these two regions also appear particularly sensitive to assumptions about burial and remineralization.Besides investigating the nutrient and oxygen concentrations above the sea floor, we will therefore mostly focus on phosphate and oxygen in these two regions in the following.
We show the effects of benthic burial and remineralization on the nutrient and oxygen concentration near the sea floor, as depicted in Fig. 3 for models INST, CTL and BUR simulated with moderate (scenario "ref") particle sinking speed.Both model INST and CTL are very similar, as they predict far too high bottom-water phosphate concentrations in the Pacific Ocean, and, correspondingly, far too low bottomwater oxygen concentrations.This mismatch gets worse for faster sinking speed (e.g, for scenario s4; not shown), resulting in deep phosphate concentration above 3 mmol m −3 over wide areas of the North Pacific, and correspondingly low (< 30 mmol m −3 ) oxygen concentrations.Introducing burial improves the match to observed pelagic tracer fields, both with respect to phosphate and oxygen.The differences among the different scenarios decrease for slower sinking speeds, where the models show a similar fit for phosphate, and even slightly too high oxygen (above 180 mmol m −3 in scenario s1) in many regions of the deep Pacific Ocean (not shown).
The introduction of burial also affects the simulated distribution of phosphate and oxygen in upper parts of the water column, as evident from a section along the "conveyor belt" (Fig. 4).The slow sinking scenario of model CTL underestimates mesopelagic (0-2000 m) phosphate in the North Pacific, and overestimates its concentration in the deep waters of this region (Fig. 4, panel A1).Increasing transport of organic matter to the deep ocean, as in scenario "fast" (panel A2 of Fig. 4) further reduces shallow and mesopelagic, and increases deep (> 2000 m) phosphate in this region.At the same time, phosphate in the North Atlantic is reduced.Therefore, in accordance with Kriest et al. (2012) remineralization length scale in model CTL has the effect of shifting nutrients downwards, and towards older waters along the conveyor belt.Introducing burial in model BUR leads to "slow" and "fast" model solutions that are more similar to each other than the corresponding solutions of model CTL.Introducing burial has little effect on shallow and mesopelagic nutrients in both scenarios "slow" (shallow remineralization) and "fast" (deep remineralization), but reduces deep nutrients in the northern North Pacific of the "fast" scenario making them more similar to those of the "slow" scenario (panels A3 and A4 of Fig. 4).Burial thus lowers the sensitivity of deep model nutrients to changes in remineralization length scale.Likewise, the sensitivity of simulated deep phosphate in the eastern equatorial Pacific (EEP) to changes in sinking speed is reduced, as evident in Fig. 5, panels A1 to A4.As a consequence, the models' misfit to observed phosphate (calculated as global average of the root mean square error RMSE; see also Kriest et al., 2010) becomes less sensitive to the remineralization length scale when introducing burial (Fig. 6).Instantaneous remineralization of organic matter at the sea floor (as introduced in model INST) shows almost the same distribution and sensitivity of phosphate and oxygen as in model CTL (see Figs. 6 and 7).
Interestingly, the impact of burial on the oxygen distribution is even more pronounced than its impact on phosphate.While the fast-sinking scenario of model CTL shows a severe underestimate of simulated oxygen in the deep North Pacific and in the EEP, this underestimate almost disappears with the introduction of burial.Burial also considerably improves the deep oxygen in the EEP and in the northern North Pacific of the fast sinking scenario.This is also evident from the volume distributions of oxygen for the different models: in contrast to that of phosphate, which is very similar among the different models and scenarios, oxygen distributions change considerably among the burial and non-burial models (Fig. 7).While scenario "fast" of model CTL predicts far too large a volume of water with low (≈ 100 mmol m −3 ) oxygen, the detrimental effect of fast sinking disappears almost entirely when burial is considered (BUR).If burial is only weakly related to sediment input (WBUR), simulated oxygen again becomes more sensitive (compared to scenario BUR) to changes in the remineralization length scale.The lower sensitivity of simulated oxygen to changes in the remineralization length scale in model BUR is also evident from plots of the model misfit to oxygen vs. flux exponent (Fig. 6).
To summarize, adding burial yields a strong improvement of the modeled pelagic biogeochemical tracers for fast sinking speed (comparable to a long remineralization length scale, or deep penetration of flux).The effect is most pronounced in old waters such as in the northern North Pacific, or in the EEP.As a consequence, simulated dissolved tracers become less sensitive to changes in remineralization length scale, thereby reducing the potential of biogeochemical water-column data to constrain this parameter.Given this weak sensitivity of both nutrients and oxygen simulated by model BUR, we ask to what extent possible impacts of model deficiencies may be "hidden" in the sediment, and if comparison to benthic observations can compensate for the reduced sensitivity of model-data misfits in terms of pelagic oxygen and phosphorus distributions to remineralization length scale.In the next subsection we therefore examine the potential of additional diagnostics of particle fluxes and benthicpelagic coupling to identify model deficiencies and to better constrain the respective model parameters and parameterizations.(2008).The models reflect the general pattern of high particle flux in the high northern latitudes and the upwelling regions, and low particle flux in the oligotrophic subtropical gyres.As expected, decreasing the flux exponent increases particle flux at 2000 m, especially in the eutrophic regions.Simulated global particle flux as well as its pattern is quite insensitive to changes in the strength of burial (see Table 2).The slow sinking scenarios yield the lowest RMS misfit to observed flux of 112 and 113 mmol C m −2 yr −1 for runs CTL and BUR, respectively (see Table 2).However, the moderate to fast sinking scenarios of these models overestimate particle flux in the Southern Ocean (see Fig. 8).

Benthic remineralization
Because our models do not represent the sediment explicitly, defining the model's counterpart to observations of benthic remineralization is not straightforward.We may, however, approximate benthic remineralization using the following assumptions: Over long timescales every (simulated) detritus particle in the model's deepest grid box will have encountered the sediment at least once.Although the model's upstream numerics for the description of sinking assumes that this resuspended detritus is equally distributed within the grid box, we may also view it as sediment "fluff", i.e., as organic matter that is only loosely associated with the sediment, but can be easily remobilized by currents, animals, etc.We therefore define detritus remineralization vertically integrated over the deepest model box as an (upper) model estimate of benthic remineralization per sea-floor area.
Figure 9 shows simulated benthic remineralization of model BUR, together with observations as compiled by Seiter et al. (2005, we   area which otherwise would be only sparsely represented in the observations.In general, the model scenarios reflect the very high benthic remineralization rates observed along the coasts, and predict very low remineralization in the subtropical gyres.For the Northern Hemisphere the simulated low subtropical values are supported only by the few observations by Seiter et al. (2005).However, in the South Pacific, Fischer et al. (2009) observed some sites with rather high benthic fluxes, which are not simulated by the model.Another severe model-observation mismatch can be found in the Arabian Sea.Here the observations suggest extremely high benthic fluxes, whereas the model predicts rather low benthic remineralization.As the sites in the Arabian Sea are quite far from continental slopes, lateral advection does not serve to explain the apparent mismatch.We note that the very high fluxes shown by five data points collected in the Arabian Sea cannot be simulated by any of our model experiments, and therefore strongly deteriorate the model's fit to observations.However, the observations in this region have been derived from a local model fit to observed sediment properties (Luff et al., 2000), which may not correspond to assumptions inherent in the global model used in this study.In addition, mesoscale processes are not resolved in our coarse, 1 • × 1 • model, but their representation has been shown to be very important in this region (Kawamiya, 2001;Resplandy et al., 2011).Global z level models often fail to represent the complex hydrodynamic and biogeochemical structures of the Arabian Sea (Dietze, personal communication, 2013).We thus have skipped these observational estimates from model comparison.The fit of the models with pronounced burial (BUR, DUNNE) to observations of benthic remineralization improves with faster sinking speed (Table 2 and Fig. 9).This is in contrast to the fit to observations of the particle flux at 2000 m described above, which tends to be best for slow sinking speeds.This can be ascribed to higher benthic remineralization with increasing flux to the ocean floor, as well as to a steeper gradient between regions of simulated low and high benthic fluxes.If burial is weaker, as in model WBUR, the misfit to observations of benthic respiration is larger (Table 2), with a minimum misfit at medium sinking speed.Without burial (CTL), the faster sinking leads to an overestimate in benthic remineralization in many regions, resulting in a worse fit (Table 2).

Burial in the sediment
We compare simulated burial flux to the observations described in Appendix A, which provides more details about the individual data sets.As in both the model and the observations, burial is defined by the fact that organic matter is removed from the system for a very long time, the modeldata comparison appears more straightforward than for benthic remineralization.
When simulating medium-to fast sinking speed, model BUR often overestimates observed burial (Fig. 10), except for the coastal regions off Washington, California and Mexico.The overestimate is most severe in the EEP.In addition, slow sinking together with the insufficient representation of shelf regions leads to a quite low simulated burial of organic matter off California and Washington, and therefore to an increase of the RMS misfit.Therefore, although simulated burial mitigates the otherwise too high an oxygen consumption by remineralization in the EEP under medium-to fast sinking speed, it results in unrealistically high burial of organic matter below the sediment surface especially in this region.Compared to models with weak (WBUR) or -not surprisingly -no (CTL) burial, models BUR and DUNNE are in better RMS agreement with observational estimates of global burial (Table 2).

Global fluxes
Despite the sometimes considerable model-data mismatches in various regions, some of the models perform quite well when compared to global fluxes of organic matter (Table 2).Benthic and riverine fluxes in all models increase with decreasing flux exponent (= increasing particle sinking speed).The strength of burial has only little influence on particle flux in 2000 m, but global burial and therefore riverine input is higher in models with strong burial (BUR and DUNNE),  while benthic remineralization is higher in the model with weak burial (WBUR).Despite the relatively fine (1 • × 1 • ) spatial resolution, the models do not resolve the shelf very well.This is reflected in a smaller difference between total and deep (> 2000 m) simulated benthic remineralization and burial than found in the observations (Table 2).While the models simulate 40-66 % of total benthic remineralization, and 32-80 % of total burial to occur below 2000 m, observations suggest lower percentages for the "open ocean" (i.e., shelf and slope region omitted), of 26 and 29 % for remineralization and burial, respectively.Because of this deficient representation of coastal and shelf areas, in the following we focus on the deep (< 2000 m) fluxes.

Biogeosciences
Global particle flux in the water column at 2000 m depth is relatively insensitive to the implementation of burial, thus all models show about the same response to changes in particle sinking speed (Table 2).Particle flux at 2000 m is reproduced best with slow sinking speed (scenarios "very slow" and "slow"), a finding that agrees with results obtained with a coarser resolution of this model (see Kriest et al., 2012).
The implementation of burial has a strong impact on the model's representation of global burial below the sea floor.Model WBUR and, of course, the no-burial models CTL and INST, underestimate global burial.At the same time, these models for moderate-to fast sinking speeds overestimate benthic remineralization.Taken together, this points towards too weak a burial, or too strong a "resuspension" and subsequent remineralization in these models.In contrast, the base scenarios of BUR and DUNNE match observed global burial flux much better, but underestimate global benthic remineralization (Table 2).The sinking speed plays a strong role for these models' agreement with observed fluxes of burial.

Sensitivity of dissolved tracers to particle flux and pelagic-benthic coupling
A finding that is, at first sight, surprising for a model with fixed stoichiometric relations is that simulated oxygen seems to be much more affected by changes in model structure and sinking speed than phosphate.The reason for the lower sensitivity of phosphate may be found in the fixed P-inventory imposed onto the model: although an increase in particle sinking speed will result in higher phosphate concentration in deeper and older waters, at the same time this increase is, in our spun up steady-state solution, limited by the concomitant phosphate decrease in younger waters (see also Kriest et al., 2012).Thus, the misrepresentation of phosphate in the model is limited by the constraint of phosphorus mass conservation.We expect results to be very similar in models with dynamic phosphorus inventory, because of the residence time of phosphorus being long compared to the ocean overturning timescale (Paytan and McLaughlin, 2007).In contrast, the oxygen inventory can change due to the combined effects of biogeochemistry, circulation and mixing, and air-sea Biogeosciences, 10, 8401-8422, 2013 www.biogeosciences.net/10/8401/2013/gas exchange with the atmosphere, so that the model cannot only fail in the internal distribution of oxygen, but also in its average concentration.This is consistent with the much shorter residence time of oxygen than phosphorus, with the former being of the order of that of the overturning circulation, as indicated by mean transit times from the surface to the deep North Pacific of 1360 ± 350 yr (derived from observations; Khatiwala et al., 2012).Model simulations also suggest a mean age between ≈ 1000-1600 yr for this region (Khatiwala, 2007).The introduction of burial reduces the sensitivity especially of simulated oxygen distribution to changes in the particle flux exponent.The models with fast particle sinking now perform much better with respect to oxygen, the reason being most likely, that now "excess" organic matter, whose remineralization would strongly decrease oxygen, is buried in the sediment.In other words, it is "swept under the carpet", and only reappears as phosphate in regions far away from the areas of burial.Note that many of the worlds major rivers discharge into the Atlantic Ocean (Fig. 13).While we cannot exactly define the pathways of phosphorus atoms between the Pacific and Atlantic Ocean, we can calculate the loss and gains through boundaries for each basin: about 1 / 3 (28.7 Gmol) of phosphorus buried in the Pacific (north of 40 • S) is not resupplied into this ocean.This corresponds to an oxygen demand of 4.9 Pmol (if all of the organic matter buried was respired instead).In contrast, the Atlantic Ocean receives more that twice as much phosphorus via river runoff (177.4Gmol) than it looses via burial (77 Gmol).Alltogether, the Atlantic receives about 67 % of global runoff, but contributes only to 29 % of global loss via burial.These figures suggest that a fraction of phosphorus buried in the Pacific and elsewhere is "teleported" to the younger waters of the Atlantic by the way we parameterize these processes (and mass conservation) in our model.

Comparison to vertical fluxes
This reduced sensitivity of simulated dissolved tracers in the burial models poses problems for any method that aims at constraining the model parameters via calibration against global pelagic tracer distributions.A potential solution to this problem could be to use additional observations, such as pelagic and benthic fluxes of organic matter, that are more closely related to particle flux and remineralization.However, some peculiarities like mismatches in the spatial and temporal scales, and methodological constraints can make such a direct comparison between models and observations difficult.
It is beyond the scope of this paper to discuss the methodological problems associated with sediment traps (but see, e.g., Honjo et al., 2008, and citations therein).However, even with a perfect particle interceptor trap we would still have to deal with the spatial hydrographic variability around many trap locations, which may hinder a direct model-data com-parison in regions of strong horizontal gradients (e.g., Siegel et al., 2008).However, the general overestimate of particle flux by the fast sinking models as evident in Fig. 8 indicates -together with the other metrics -that these are less likely to represent pelagic biogeochemical processes well (see also Kriest et al., 2012, for the relation between fit to estimated global flux and dissolved tracer).
Benthic C-org remineralization or dissolved oxygen utilization ("DOU") is measured either (a) from pore water profiles and diffusion models, or (b) via the oxygen decline over time in benthic chambers.While method (a) can lead to decompression and handling artifacts (much of which has been overcome by using in situ probes and measurements), (b) can introduce some variability through benthic fauna activities (e.g., Glud et al., 1994).Further, while (a) represents a steady-state system (due to the intrinsic assumptions), the benthic chamber used in (b) represents a snapshot of sediment respiration in an isolated system that excludes both transfer of organic matter to the sediment, as well as exchange with the surrounding water and dissolved oxygen.The latter may be of importance when observing systems at low to vanishing oxygen concentrations, where oxygen decline may affect oxygen sensitive processes and rates.However, assuming that during the incubation time of typically a few days, the rate of respiration is neither carbon nor oxygen limited, this observation may well represent a steady state system.
Burial is usually estimated directly from 14 C age, mass accumulation rate and % C org of the sediment cores, meaning that these flux estimates assume steady state over multiple thousands of years.Perhaps counter-intuitively, simulated and observed burial estimates compare quite well, strengthening the case for both the validity of the steady-state assumption over these timescales and fidelity in the models.Unfortunately, as we have seen above (Fig. 10) these observations are very sparse, and biased towards coastal areas and/or highly productive regions.In addition, there may be complications and artifacts due to sediment focusing and erosion (e.g., Kienast et al., 2007), or potential age offsets among different sediment weight fractions (Heinze et al., 2009, and citations therein).
Overall, the spatial sparsity of the flux observations, their potential bias towards certain regions, as well as the sometimes difficult assignment of proper model counterparts to the observations make it difficult to use them as strong constraints on model performance.However, the models' overestimate of simulated burial especially under medium-to high sinking rates in the eastern equatorial Pacific (EEP) points towards a potential lack or misrepresentation of processes by the models in that region.
An explicit representation of sediment processes, inventories and fluxes, as model level 3 or 4 suggested by Soetaert et al. (2000), or as implemented in global models by Heinze et al. (1999)  fluxes, and provide more insight on the impact of benthic processes on the (pelagic) ocean, and vice versa.Likewise, a better representation of bottom topography would likely improve the simulation of benthic processes, especially on the coastal shelf.While this is beyond the scope of this study, we consider this approach very useful.We thus regard the flux observations as useful additional model constraint, that may complement any model assessment based on more "traditional" observations such as nutrients or oxygen.

Model performance in the eastern equatorial Pacific and in coastal areas
A potential explanation for the apparent misfits in the EEP could be the neglect of oxygen sensitivity of remineralization in the model.This neglect currently allows for complete consumption of oxygen at unaltered remineralization rates, whenever the supply of organic matter is high enough.
Thereby the current model assumes a form of "implicit" denitrification, as remineralization continues even in the absence of oxygen.An alternative explanation for the mismatches in the EEP could be a misrepresentation of the Equatorial Intermediate Current System in this region, resulting in too little oxygen (Dietze and Loeptien, 2013).According to that study, and citations therein, the problem is "endemic even to eddyresolving circulation models".Improved representations of physical processes in regional or (nested) global models might help to resolve this problem.
Considering only areas deeper that 2000 m, the models with strong burial seem to be better suited to represent global vertical fluxes, although results for particle flux at 2000 m, benthic remineralization and burial point towards different "best" candidates.It remains to be investigated, whether a better topographic representation of coastal and shelf areas helps to achieve a closer match to local observations of benthic processes, as well as to global estimates of burial and benthic remineralization.

Combining model metrics
So far it seems that the addition of burial has been of little benefit to the model in terms of its ability to reproduce watercolumn distributions of biogeochemical tracers.A main effect of burial in the model is a weakening of the constraints provided by water-column tracer distributions on the remineralization length scale, thereby making the model more robust with respect to errors in the parameterization of particle flux.At moderate to slow sinking speeds all models represent dissolved tracers about equally well.
In order to identify a "best" model type, we have evaluated model performance using a combined metric of dissolved tracers and vertical fluxes.The metric consists of terms for normalized root mean square errors (RMSE) for phosphate (r P ), oxygen (r O ), particle flux (r F ), benthic remineralization (r R ) and burial (r B ).We further complemented the analysis of model skill by calculating the relative deviation of simulated to observed global oxygen inventory (d O ; expressed as global average oxygen) as well as simulated to observed global particle flux (d F ), benthic remineralization (d R ), burial (d B ) and river runoff (d r ).As noted above, because of mass conservation, the global phosphate inventory does not help to constrain the model, and is therefore not used for model assessment.To add quantities of different units, we have normalized the RMSE and deviations of oxygen by the respective observed average global concentration.RMSE and deviations of vertical fluxes have been normalized by the respective global fluxes (see Table 2).Total model skill S is thus evaluated as Table 3 and Fig. 11 depict S of the different model experiments, as well as its components.As already noted above, the misfit to dissolved tracers shows only small differences among the models, particularly for phosphate, and for slow to moderate sinking speeds.Differences between simulated and observed oxygen inventory (d O ) indicate a considerably better performance of models with slow sinking speed, even for the burial models.
For the two models CTL and WBUR with no or only weak burial, metrics according to particle flux (RMSE as well as deviation of global flux) coincide with those from dissolved tracers, similar to the results obtained by Kriest et al. (2012).Models perform almost identical with respect to particle sinking speed.Thus, comparison to particle flux supports the decision for slower particle sinking speed, but does not help to decide among different model types.
Models BUR and DUNNE with rather strong burial in most cases perform better with respect to the RMSE for benthic fluxes than do models CTL and WBUR with no or weak burial.The pattern is more complex for the deviation to global fluxes.Here, models BUR and DUNNE when simulated with fast or very fast sinking speed exhibit a strong overestimate of global burial, and therefore a stronger mismatch than model CTL or WBUR, but perform much better than CTL when simulated with the medium sinking speed of scenario "ref".At the same time, fast sinking speed in BUR and DUNNE leads to a better agreement with observed global benthic remineralization.
The misfits for burial r B and global burial d B are quite large and therefore contribute strongly to the overall metric, S. As a result, given the good fit to burial of the reference scenarios of models BUR and DUNNE, together with the quite similar fit of all models to dissolved tracers, scenario "ref" of BUR and DUNNE performs best with respect to S.Even the second best scenario "slow" of these two model types performs much better than any experiment of CTL or WBUR.
Summarizing, despite overly high burial rates in certain regions, we regard models BUR or DUNNE as more suitable than CTL or INST without any pelagic-benthic coupling, or Biogeosciences, 10, 8401-8422, 2013 www.biogeosciences.net/10/8401/2013/model WBUR with only weak burial at the seafloor.The little difference between simulations BUR and DUNNE may be attributed to the fact that both burial functions are quite similar in regions of high particle rain, which strongly impact global bulk fluxes and misfit.
6.5 Deep detritus: an "escape route" for excess organic matter?
So far, we have investigated dissolved tracers in the water column, as well as vertical and benthic fluxes.What remains is to examine the last potential "escape route" for excess organic matter that sinks to the sea floor, namely detritus in the deep ocean.
Our suite of model experiments includes the most extreme case for benthic remineralization, namely model "INST", in which all organic matter arriving at the sea floor is remineralized instantaneously.This model setup is the same as model "level 2" presented by Soetaert et al. (2000).Its dissolved pelagic properties are quite close to that of model CTL (Fig. 3), indicating that in the long term these two exhibit similar bulk fluxes in the bottom layer and at the benthicpelagic interface.
Differences arise, however, when implementing benthic burial, where models with burial show much lower bottom phosphate, and higher oxygen concentrations near the sea floor (Fig. 3).In models BUR and DUNNE we allowed particulate organic matter to be locked away in the sediment, thereby mimicking burial and/or the activity of benthic organisms.The lower phosphate and higher oxygen concentrations in the bottom box the two "best" models DUNNE and BUR come not only at the cost of partly too high a burial, but are also accompanied by elevated deep detritus concentrations.However, we argue that these are not unrealistically high.While model INST due to its inherent assumptions shows decreasing detritus concentrations with depth (Fig. 12), model CTL exhibits elevated (in contrast to the layer above) detritus at the sea floor.In some areas detritus increases towards the seafloor even more than 10-fold (Fig. 12).Burial reduces this increase towards values that are four to seven times that of the overlying layer, thereby mimicking a benthic nepheloid layer (BNL).
BNLs are indicated by elevated turbidity and have been observed in many different regions (Vangriesheim et al., 2001;Inthorn et al., 2006;Lukashin and Shcherbinin, 2007;Capello et al., 2009;Fischer et al., 2009).Typically, the BNL has a thickness up to tens to hundreds of meters (e.g., Turnewitsch and Springer, 2001;Inthorn et al., 2006).The amount and spatial extension of these remobilizations depends not only on physical processes, but also on the megafauna composition, the type of the material deposited (e.g., Vangriesheim et al., 2001;Lampitt et al., 2001), or the abundance and activity of groundfish (Yahel et al., 2008).
The large relative increase of detritus in the bottom water compared to the overlying grid box exhibited by model CTL seems to be rather high when compared to observations, but model BUR with its less pronounced particle increase seems to be more in line with observations of a several-fold increase of particle number or mass taken e.g by Boetius et al. (2000), McCave et al. (2001), or Lukashin and Shcherbinin (2007).Note that despite the several fold increase of detritus concentration relative to the overlying water, the absolute concentrations in deep bottom waters are mostly at nanomolar levels.The relative increase of detritus in the bottom layer, compared to the layer above, declines with slower sinking speed (no figure), and increases with faster sinking speed, particularly for model CTL.Thus, differences between CTL and BUR decrease for slow sinking speeds.
An alternative to the reflection of detritus at the sea floor could be its immediate remineralization, as exemplified by model INST, and suggested by Soetaert et al. (2000, their model level 2).As shown above, using the current, rather coarsely resolved model, that continues remineralization even under suboxic conditions, its impact on pelagic dissolved tracers is quite small.It remains to be investigated, how the difference between types INST and CTL would impact model outcome when applied with oxygensensitive remineralization (where suboxic conditions would hinder aerobic remineralization).Further, cross-or alongshelf transport of organic matter might play a larger role in connecting the shelf and the open ocean in models with a finer spatial resolution.
As shown above, some scenarios of models BUR and DUNNE perform best with respect to a combined metric of dissolved tracers and fluxes.By avoiding extreme tracer concentrations near the seafloor, they might therefore serve as a good starting point for a later parameterization of oxidant sensitive processes such as oxygen-dependent remineralization or denitrification.

Conclusions
In contrast to the constraint of mass conservation imposed on phosphorus, exchange of oxygen with the atmosphere results in variable oxygen inventory.Therefore, simulated oxygen is more affected by changes in model structure and sinking speed than phosphate.The introduction of burial reduces the sensitivity of simulated oxygen and, to a lesser extent, phosphate distributions to changes in the particle flux.
In some regions remineralization of organic matter, in particular in conjunction with high sinking speed, would cause a severe oxygen deficit.In models that allow burial of detritus in the sediment, high concentrations of organic matter in deep layers, associated oxygen deficit and therefore high sensitivity towards changes in remineralization length scale are all "swept under the carpet".et al. (1996) for rivers south of 60 • N and mapped onto the topography of the MIT2.8 degree model as described in the text.Symbols denote original locations from the data set by Perry et al. (1996).Crosses: omitted from runoff calculation (see text).Open squares: included in runoff calculation.

Biogeosciences
53 Fig. 13.Distribution of river runoff (as fraction of total runoff, colored boxes), as calculated from Perry et al. (1996) for rivers south of 60 • N and mapped onto the topography of the MIT2.8 degree model as described in the text.Symbols denote original locations from the data set by Perry et al. (1996).Crosses: omitted from runoff calculation (see text).Open squares: included in runoff calculation.
The resulting weakened sensitivity of simulated dissolved tracers to the remineralization length scale can be partly overcome by using vertical fluxes as model constraints, but mismatches in the spatial and temporal scales, and methodological constraints may hinder the direct comparison between models and observations.However, observational flux estimates may serve as useful indicators of potential model deficiencies.For example, the strong overestimate of simulated burial under medium-to high sinking rates in the eastern equatorial Pacific (EEP) points towards a potential lack or misrepresentation of processes by the models in that region.Potential candidates for model improvements particularly in the EEP could be an improved representation of physical pro-cesses in finer scale regional or nested global models, and/or the (explicit) parameterization of oxidant-sensitive processes under oxic and suboxic conditions.
Given the robustness of the burial models to observed dissolved tracers, and their relatively good match to many of the observed fluxes, we regard the burial models using parameterizations BUR and DUNNE as most suitable candidate for further studies investigating processes such as oxygendependent remineralization or denitrification.
A thorough and comprehensive examination of the impact and potential feedback processes of the sediment on ocean inventories and fluxes can only be carried out by coupling an explicitly simulated sediment to a global biogeochemical model (as, e.g., in Heinze et al., 1999;Maier-Reimer et al., 2005).This is beyond the scope of this study.However, given the importance of exchange processes across the lower model boundary demonstrated in this study, and the benefit gained from comparison to benthic observations, we consider this to be a useful addition to the model.Honjo et al. (1992).
In addition to these eight data sets, we have further added the data set by Berelson et al. (1997), which gives burial and remineralization rates for the equatorial Pacific along 140 • W. The data set by Kienast et al. (2007) only gives burial rates for the Panama Basin.As this region is especially sensitive to ventilation and biogeochemical parameters (especially DOP decay parameters; see Kwon and Primeau, 2006), we have added this data set as a valuable constraint on simulated burial.We supplemented it by carbon rain rate to the sea floor, derived from sediment traps deployed in the vicinity of this site (Honjo et al., 1992).Taken together, we therefore have 10 different data sets: 1. Bender and Heggie (1984) -Pacific and Atlantic: we used their Table 12 and station locations from their Table 1.For the site in the eastern equatorial Atlantic we only considered the five stations nearest the equator.We only considered C org -Oxidation from O 2 .
2. Berelson et al. (1987) -Californian Borderland (S.Pedro and S. Nicolas): we used their Table 8 for remineralization and burial rate of C org .Station locations were determined from their Fig. 1, and station depths from their Table 1.Station locations were calculated from averages of all stations listed in their Table 1.
4. Jahnke (1990) -Californian Borderland (S.Monica and S. Catalina): we used Table 3 from Jahnke (1990), but only values for S. Catalina and S. Monica.Station locations were read from their Fig. 1.For depth we used the value of the nearest isoline.
5. Berelson et al. (1997) -Equatorial Pacific: station locations and fluxes were taken from their  12.For remineralization we used the "NO − 3 Remin" value.For a range of fluxes, we took the average (center) of this value.For burial, we used the first ( 230 Th-normalized) rate.
8. Devol and Hartnett (2001) -off Washington and Mexico: Table 1 from Devol and Hartnett (2001) contains 22 entries.We used averages over i = 1 : 13 for the region off Washington State, and for i = 14 : 22 for the region off Mexico, for location, depth and fluxes of rain rate, burial and remineralization.9. Stahl et al. (2004) -Porcupine Abyssal Plain (PAP): we read flux for station PAP from their Fig. 6.Station location was taken from their Table 1 (average over all values).10.Kienast et al. (2007) and Honjo et al. (1992) -Panama Basin: burial data (Acc.rate) were read from Fig. 8, carbon flux from mass accumulation rate.Station locations were taken from their Table 1.
For core P7, we have added the trap-derived sedimentation rate of 10.87 mg C m −2 d −1 at 2869 m observed by Honjo et al. (1992, their Table 5; "Honjo -A").Note that this value is probably an upper estimate of the carbon rain rate to the sea floor, in the vicinity of this site.
Table A1 shows the fluxes converted to mmol C m −2 yr −1 .For the plot shown in Fig. 1 and corresponding regressions, we used remineralization + burial as total flux to the sea floor.In case where there was no remineralization available (data by Lohse et al., 1998 andKienast et al., 2007) we used the rain rate to the sediment as upper estimate for remineralization rate.Note that in these two data sets, burial is just a small fraction of rain rate -therefore, subtracting it from the rain rate would have made only a small difference to the regression.

Data set for river runoff
We use the data set of volumetric flow rates given in Perry et al. (1996).Because little seems to be known about their nutrient contents of Arctic rivers (but see the recent paper by Holmes et al., 2012), and because during model simulation the nutrient supply via these rivers may get trapped in the Arctic, we excluded the 14 rivers that discharge north of 60 • N, namely Yenisei, Lena, Ob, MacKenzie, Yukon, Pechora, Severnaya Dvina, Khatanga-Popigay, Kolyma, Pyasina, Indigirka, Taz, Kuskokwim, and Copper.For the remaining 61 rivers, we calculated the minimum distance of their mouth's location to the corresponding model location in the MIT2.8 configuration.We only included rivers which have a corresponding "wet" point in the model within a distance of two times the model grid resolution (2.8 • ), times √ 2. Because the MIT2.8 model does not resolve many marginal seas, this results in the omission of the following rivers from the data set: Danube, Dniepr (Black Sea); Nile, Po, Rhone (Mediterranean Sea); Nelson, La Grande (Hudson Bay); Neva (Baltic Sea); Shatt al-Arab (Persian Gulf); Huang Ho (Yellow Sea).
In order to account for the potentially large fan of river runoff, we first calculated the number of all surface model boxes that are influenced by a river in the modified data set.In case a river's discharge exceeds half of the corresponding model surface box volume per year, we extended the region over which this discharge is distributed equally over several surface boxes.In the MIT2.8 model geometry this only affects the Amazon outflow, whose discharge is being distributed over 11 horizontal grid points around its mouth.The resulting runoff (discharge volume of a river, expressed as fraction of total discharge of all rivers), and all river locations south of 60 • N from the original data set are shown in Fig. 13.
For the parameterization of runoff fields for ECCO, we use the MIT2.8 runoff volumes, remapped onto the ECCO geometry.By doing so, we assure that the riverine input of buried P happens at approximately the same location in both model configurations.Burial and runoff for both model types show about the same magnitude and distribution, with some difference along the coasts.Although not shown here, the general response of the MIT2.8 model to the introduction of burial and runoff is the same.

Fig. 1 .
Fig. 1.Burial vs. rain rate onto sediment for various approaches.Thick red symbols: data compilation for regions with burial < 250 mmol C m −2 yr −1 .Thin red symbols: data compilation for regions with burial ≥ 250 mmol C m −2 yr −1 .Straight red line: regression for data compilation over all regions (used in model scenario BUR).Dashed red line: regression for data compilation over regions with burial < 250 mmol C m −2 yr −1 (used in model scenario WBUR).Regressions from Flögel et al. (2011) are shown as straight (open ocean) and dashed (coastal) blue line.The algorithm by Dunne et al. (2007) is denoted as green line.Thin black lines denote constant burial efficiencies of 100 (straight), 10 and 1 % (both dotted).Observations (symbols) and lines for algorithms byDunne et al. (2007) andFlögel et al. (2011) were given in carbon units, and have been converted to phosphorus using a constant C : P ratio of 106.

Fig. 2 .Fig. 2 .Fig. 3 .
Fig. 2. Volume distributions of global phosphate (upper, A and C) and oxygen (lower, D and F) for model CTL in the MIT2.8 (left) and ECCO (right) configuration.Grey bars denote the corresponding observations (WOA05).Lines denote the different experiments.Thin lines: fast sinking profile.Thick lines: slow sinking profile.Medium line: reference scenario (see Table 1 for model configurations).Mid panel: Misfit of phosphate (upper) and oxygen (lower).Thin lines: MIT2.8 configuration.Thick lines: ECCO configuration.Horizontal lines: Global spatial variance (square root) of observations.42

5
Particle flux in the water column and at the sediment-water interface 5.1 Particle flux in the water column In Fig. 8 we compare the particle flux in 2000 m simulated by the different sinking scenarios of model BUR against observations compiled by Honjo et al.

Fig. 6 .Fig. 6 .
Fig. 6.Misfit function, calculated as RMSE, divided by observed global standard deviation, for phosphate (upper panels) and oxygen (lower panels).Left: Misfit for northern North Pacific (north of 40 • N).Right: misfit for eastern equatorial Pacific (between ±10 • and east of 140 • W).Middle: global misfit.Black lines: CTL; Thin black lines with open circles: INST; Green small crosses and lines: WBUR; Blue large crosses and lines: BUR.Turquoise pluses and lines: DUNNE.

Fig. 7 .
Fig. 7. Volume distributions of global phosphate (upper panels) and oxygen (lower panels) for different model experiments.Grey bars denote the corresponding observations (WOA05).Lines denote the different experiments.Thin lines: fast sinking profile.Thick lines: slow sinking profile.Medium lines: reference scenario.

Fig. 7 .
Fig. 7. Volume distributions of global phosphate (upper panels) and oxygen (lower panels) for different model experiments.Grey bars denote the corresponding observations (WOA05).Lines denote the different experiments.Thin lines: fast sinking profile.Thick lines: slow sinking profile.Medium lines: reference scenario.

Fig. 8 .
Fig. 8. Upper panels: particle flux in 2000 m, for experiment BUR with burial, and observations.Sinking speed in models increases from the left (slow) to the right (fast).Small colored circles denote observations from the data set by Honjo et al. (2008), on the same color scale as the model results.Lower panels: simulated vs. observed particle flux.Lines indicate 1 : 2, 1 : 1 and 2 : 1 relationship.The panels also denote the RMSE of the modeled vs. observed flux, as well as the number of data points sampled.

Fig. 9 .Fig. 10 .
Fig. 9. Upper panels: benthic remineralization (see text) for experiment BUR with burial, and observations.Flux exponent in models decreases from the left (slow) to the right (fast).Colored symbols denote observations from the data set by Seiter et al. (2005, circles) and Fischer et al. (2009, diamonds), on the same color scale as model results.Lower panels: simulated vs. observed benthic remineralization.Different colors and symbols denote different depths and/or data types: green small circles: 1000-2000 m.Turquoise large circles: 2000-3000 m.Blue triangles: 3000-4000 m.Black crosses: > 4000 m.Pink stars: Arabian Sea (> 3000 m).Data byFischer et al. (2009, South  Pacific)  are denoted by large blue triangles (3000-4000 m) and large black crosses (> 4000 m).The panels also denote the RMSE of the modeled vs. observed benthic remineralization (data set bySeiter et al., 2005), and the number of data points sampled (upper numbers: without Arabian Sea data, lower numbers: total data set).
or Maier-Reimer et al. (2005) can facilitate and improve the comparison of simulated to observed benthic www

Fig. 11 .
Fig. 11.Metrics for different model types plotted vs. different sinking speed.Panel B: sum of misfit (RMS, divided by observed global average concentration) for phosphate and oxygen; Panels C-E: normalized misfit to data sets for particle flux in 2000 m (C), benthic remineralization (D), benthic burial (E).Bottom panels F-J: normalized deviation between simulated and observed global inventory of oxygen (F), between global fluxes of organic particles (H), benthic remineralization (I), benthic burial (J) and global river runoff of phosphate (G).Panel A shows the sum over all panels B to J. Normalization of biogeochemical fluxes has been carried out by global integrated fluxes listed in table 2. Colours as in figure6.

Fig. 11 .
Fig. 11.Metrics for different model types plotted vs. different sinking speed.Panel (B) sum of misfit (RMS, divided by observed global average concentration) for phosphate and oxygen; panels (C-E) normalized misfit to data sets for particle flux in 2000 m (C), benthic remineralization (D), benthic burial (E).Bottom panels (F-J) normalized deviation between simulated and observed global inventory of oxygen (F), between global fluxes of organic particles (H), benthic remineralization (I), benthic burial (J) and global river runoff of phosphate (G).Panel A shows the sum over all panels (B) to (J).Normalization of biogeochemical fluxes has been carried out by global integrated fluxes listed in Table 2. Colors as in Fig. 6.

Fig. 12 .
Fig. 12. Detritus concentration in the bottom layer of reference scenarios models INST (A), CTL (B), and BUR (C), expressed as multiple of detritus concentration the layer above.Note the different color scales.

Fig. 13 .
Fig. 13.Distribution of river runoff (as fraction of total runoff, coloured boxes), as calculated fromPerry et al. (1996) for rivers south of 60 • N and mapped onto the topography of the MIT2.8 degree model as described in the text.Symbols denote original locations from the data set byPerry et al. (1996).Crosses: omitted from runoff calculation (see text).Open squares: included in runoff calculation.

Table 3 .
Different misfit metrics, normalized by global average concentration of dissolved tracers (2.17 mmol P m −3 for phosphate, 174.31 mmol O 2 m −3 for oxygen) or by global annual flux (see Table2), for scenarios "slow" to "fast" (s2 to s3, corresponding to increasing remineralization depth) of different models.See text and Eq.(3) for more information.

Table A1 .
Latitude ( • N), longitude ( • W), depth (m), rain rate, remineralization and burial (all in mmol C m −2 yr −1 ), data set Id and author/location for the data set used for regression and model assessment.
* Rain rate determined from trap-derived sedimentation rate by Table 1, but only for stations where all fluxes (rain rate, remineralization and burial rate) are available.If there is more than one value, we calculated the average.6. Lohse et al. (1998) -OMEX: we used Table 6 of Lohse et al. (1998), averaged over regions "shelf" (A), "upper slope" (I, B, II), "lower slope" (C, F, III) and PAP (E).Station locations, listed in Table 1, were averaged accordingly.7. Sayles et al. (2001) -Ross Sea: averages of station locations were taken from Table 1.Fluxes were extracted from Table