A mechanistic particle flux model applied to the oceanic phosphorus cycle

Introduction Conclusions References Tables Figures


Introduction
The settling of organic particles to the deep sea has a profound effect on global ocean properties.It sustains complex and diverse benthic and mesopelagic food webs, sequesters vast quantities of nutrients and CO 2 away from the surface ocean and atmosphere, and creates a low-O 2 layer that restricts marine habitat.The vertical distribution of organic matter decomposition in the dark ocean together with ocean circulation determines where, when, and in what proportions regenerated nutrients and carbon will reemerge at the sea surface.If decomposition occurs deep in the water column, the regenerated nutrients may be stored for centuries, resurfacing at high latitudes where utilization by phytoplankton is relatively weak.For decomposition that occurs shallower in the Published by Copernicus Publications on behalf of the European Geosciences Union.

T. DeVries et al.: Mechanistic global ocean particle model
water column, the resupply will occur on faster timescales of seasons to decades, and may follow pathways to lower latitudes, where nutrient consumption is complete.Thus, the depths at which particles sink before being remineralized may have a large influence on marine productivity and carbon pump efficiency (Boyd et al., 2008;Buesseler and Boyd, 2009;Kwon et al., 2010).
The depth scale of decomposition depends on the ratio of the particle sinking velocity and the rate at which the particles are remineralized by bacteria (Sarmiento and Gruber, 2006).For faster sinking speeds or slower remineralization rates, nutrients will be released in deeper waters.However, these two critical rates may themselves be coupled, because particle sinking speeds are strongly influenced by particle size, and particle size is altered by biological rates of decomposition.Given the complexity of these dynamics, and the computational expense of simulating numerous particle size classes, even most sophisticated ecosystem/biogeochemical cycle models still treat particles implicitly, through highly idealized empirical relationships (Martin et al., 1987;Armstrong et al., 2002).Quantitative and predictive understanding of the underlying biological transformations and their environmental sensitivities remain quite primitive.
We present here a size-resolved model of marine particle dynamics to predict how the particle size spectrum changes as it sinks, depending on the characteristics of surface particles and their alteration by subsurface microbial decomposition.The model is based on a general mechanistic equation governing particle dynamics, which has been widely applied in meteorology for precipitating clouds (Hu and Srivastava, 1995), and in oceanography for both size-resolved bubble populations (Liang et al., 2013) and sinking particles (Burd and Jackson, 2002;Stemmann et al., 2004a;Kriest andEvans, 1999, 2000;Gehlen et al., 2006).The latter models compare predicted particle fluxes to measurements of particle mass or flux from sediment traps.We take a different approach, and evaluate the particle flux model using climatological nutrient distributions, which reflect the long-term spatial patterns and rates of remineralization.We do this by embedding the particle model in a 3-D ocean biogeochemistry model of the marine phosphorus (P) cycle to investigate the influence of particle sinking and respiration on global nutrient distributions and fluxes.The role of remineralization depth has been examined in a variety of models before (e.g., Kwon and Primeau, 2006;Kwon et al., 2010;Kriest et al., 2010Kriest et al., , 2012)).Our study has three advantages.First, it uses a mechanistic formulation of particle dynamics, so that the parameters can be interpreted and validated against laboratory and field observations.Second, we use a circulation model whose ventilation rates are constrained by radiocarbon and CFC observations, which allows errors in nutrient fields to be attributed to biases in biogeochemical processes, and not physical ones (Doney et al., 2004;Najjar et al., 2007).Finally, we perform a large number of steady state simulations, so that the sensitivity and uncertainty can be well characterized.

A size-resolving particle sinking and decomposition model
We begin by presenting a general framework for modeling the flux of particulate organic matter (POM), starting from a population of particles of different sizes falling through the water column.The size distribution, n (unit: number per volume per size increment), of particles with a spectrum of diameters, D, undergoing gravitational settling at sizedependent velocity, w s , and shrinking due to remineralization at rate dD/dt, evolves over time at a fixed location by: where the terms on the right-hand side represent the divergence of the particle flux, particle remineralization, coagulation, and fragmentation, respectively.The particle flux is achieved through both fluid velocity, u = [u f v f w f ], and the sinking rate, w s , of the particles.Our focus will be on particle fluxes in the context of the long-term, large-scale general circulation of the ocean interior.Accordingly, we make three simplifications.First, for particles with sinking speeds of order 10 m d −1 , the transport divergence can be reasonably approximated by the vertical particle velocities (i.e., w s w f and horizontal length scales much greater than vertical length scales).Second, we assume that the particle size distribution is in steady state throughout the water column.Finally, we focus on regions below the turbulent boundary layer where z (= z − z s ) > 0, with z s a nominal mixing depth and z defined positive downwards.Here the fragmentation and coagulation terms are relatively small and can be neglected (Boehm and Grant, 2001), although some studies suggest that coagulation can be an important process governing the vertical particle flux below the mixed layer (e.g., Stemmann et al., 2004b).Under these assumptions Eq. ( 1) simplifies to which states that the divergence of the flux of particles of a given size is balanced by the conversion of particles from larger size classes to smaller ones.Solutions to this particle equation depend on the sinking rate, the rate of mass loss (i.e., dD/dt), and boundary conditions on the particle size distribution, n, in the surface mixed layer (i.e., z ≤ 0).The sinking rate of particles depends strongly on size.Individual plankton sinking rates have been measured for a variety of species as a function of cell size, usually measured in equivalent spherical diameter (Smayda, 1970(Smayda, , 1971;;Stemmann et al., 2004a, Fig. 2).The data are commonly fit by a power law, which we adopt here: Sinking rates generally do not increase as quickly with size as do terminal velocities predicted by the Stokes law (i.e., η < 2).Several other factors also influence sinking speeds, including the density and shape of particles (McDonnell and Buesseler, 2010).Dead/senescent cells sink much more quickly than living ones, indicating an important role for motility and buoyancy regulation.Resolving all of these factors affecting particle settling speed is beyond the scope of the present study.We assume that variations in these factors are effectively averaged over the scales of interest.
The rate at which particles lose mass involves a complex set of processes by which particles are grazed by filter feeders, and colonized by free-living bacteria that hydrolyze organic matter, releasing dissolved organic matter (DOM) into the surrounding water.We simplify the biological dynamics by assuming that, in each size class, the rate of mass loss is proportional to the particle mass, an assumption that is supported by measurements on phytoplankton aggregates (e.g., Iversen and Ploug, 2013, Fig. 3).We further assume that the mass of particles increases with size according to where ζ may be less than 3 to account for the increase in fractional water content of larger particle aggregates (Alldredge and Gotschalk, 1988).Combining Eqs.
(3)-( 5) gives An analytic solution to Eq. ( 2) can be obtained assuming that the size distribution of particles in the mixed layer (above z = 0) can be described by a power law (Sheldon et al., 1972;Jackson et al., 1997), with n o a constant that determines the total mass, and the value of determining the size distribution of particles in the well-mixed surface layer.Estimates of surface particle size distribution from satellite observations use the same powerlaw formulation (Kostadinov et al., 2009), and can therefore be used to incorporate spatially variable in a global biogeochemical ocean model (see Sect. 3).Under this assumption, a solution to Eq. ( 2) can be found using the method of characteristics,  (Kostadinov et al., 2009), η = 1.17 (Smayda, 1970), ζ = 2.28 (Mullin et al., 1966), c r = 0.03 d −1 (Kriest and Oschlies, 2008).
According to Eq. ( 8), the size distribution of particles increases with decreasing particle size, and decreases with increasing depth below the mixed layer (Fig. 1a).The total mass of particles can be calculated at any depth by integrating the particle mass, m, and the particle size distribution, n.
Near the surface, the bulk of the total particle mass is contained in small particles, but this peak shifts towards larger particles deeper in the water column (Fig. 1b).At intermediate depths, the particle mass reaches a maximum at intermediate sizes (Fig. 1b).
The net conversion of mass from particulate to dissolved forms can be calculated from the mass and flux of particles integrated over the full particle size distribution.The total mass of sinking POM is thus given by and the total sinking flux (mass times velocity) of POM is given by where D S (z ) and D L (z ) are the smallest and largest particle sizes, respectively, which may vary with depth.An upper limit on the size of particles at a particular depth will be set by the size of the largest particles in the euphotic zone.The change in the upper size limit with depth can be evaluated by combining Eqs. ( 3) and ( 6) to obtain where c 1 = c r /c w .In principle, D S (z ) should vary with depth according to a formula similar to Eq. ( 11).However, we found very little sensitivity of either the total mass or the total particle flux to a change in D S with depth, and therefore for simplicity in all our calculations we set D S (z ) = D S (z = 0).Here and throughout, we assume particle sizes at the surface range from D S = 20 µm to D L = 2000 µm in the euphotic zone, and use a discretized particle size of dD = 2 µm.The total particle mass decreases strongly with depth in the first several hundred meters below the euphotic zone, due to conversion from large to small particles and the loss of particles at the upper end of the size spectrum (Fig. 1c).The total particle mass decreases approximately log-linearly with depth below about 1000 m below the euphotic zone as the mass spectrum becomes flatter (Fig. 1b).The total particle flux is heavily weighted toward the sinking flux of the largest (heaviest) and fastest-sinking particles, and therefore is not as strongly attenuated with depth as the total mass flux (Fig. 1d).
The average (mass-weighted) particle sinking velocity is defined as, For the particular combination of parameters in Fig. 1, w s increases with depth up to about 2500 m below the euphotic zone, due to the shift in the mass spectrum toward larger particle sizes at depth (Fig. 1b and d).The average particle sinking velocity then begins to decrease with depth below 2500 m as the largest particles begin to disappear completely, resulting in a decrease in the upper particle size limit (see Eq. 11 and Fig. 1a) and an overall shift toward smaller and slower sinking particles (Fig. 1d).An increase in particle sinking velocity with depth is consistent with some observations (Berelson, 2002;McDonnell and Buesseler, 2010) and is implicit in the widely used power-law used to describe the attenuation of particle flux with depth (Martin et al., 1987;Kriest and Oschlies, 2008).Here we see that the particle sinking speed can decrease with depth in the abyssal ocean if large particles begin to degrade fully.The protection of sinking POM by ballast mineral assemblages (e.g., Armstrong et al., 2002) could counter this deep trend by ensuring a supply of large particles to the deep ocean and thus a continued increase in the average particle sinking speed with depth.This possibility is addressed in Sect.4.3.The decrease in POM flux with depth depends on several parameters, which can be seen by rewriting Eq. (10) in the form where c F = n o c m c w .The actual value of c F is arbitrary, and in practice we always set c F = 1/F (z = 0) so that Eq. ( 13) is normalized by the flux at the base of the euphotic zone.Thus, there are four parameters that control the particle flux profile: ζ , η, and c 1 .We varied these parameters within plausible ranges to examine the sensitivity of the particle flux profile (Fig. 2).
The POM fluxes, F (z ), show similar sensitivities to , the exponent of the surface particle size distribution, and ζ , the exponent in the relationship between mass and particle size (Fig. 2a and b).Both an increase in and a decrease in ζ have the effect of shifting the particle mass spectrum toward smaller masses, which results in more POM being respired near the surface and less POM reaching the deep ocean.F (z ) also shows similar sensitivities to variations in η, the exponent in the relationship between particle size and sinking velocity, and c r , the degradation rate of POM (Fig. 2c and d).The relative values of η and c r control the depth that individual particles sink to before decaying, and thus have similar . Schematic illustrating the transformations between the various phosphorus pools in the model: particulate organic phosphorus (POP), dissolved organic phosphorus (DOP) and inorganic phosphate (PO 4 ).The J POP term is parameterized according to PRiSM as described in Sects. 2 and 3. effects on F (z ).Increasing η results in faster-sinking particles that penetrate deeper into the ocean, while reducing the degradation rate c r has a similar effect on the particle flux profile (Fig. 2c and d).Within this parameter space, the slope of the surface particle size distribution has the largest effect on POM flux to the deep ocean.We will therefore pay particular attention to the influence of variations in this parameter on the large-scale nutrient fluxes in the global biogeochemical model that follows.

Phosphorus cycle simulations in a global ocean model
The global and long-term distribution of nutrients such as PO 4 provides a strong constraint on the patterns and rates of remineralization implied by particle flux models.To exploit the information in these observations, and to derive appropriate parameters for the particle sinking model, we incorporate it into a global ocean circulation/biogeochemistry model.This can be done by simply using the normalized particle flux profiles Eq. ( 13) at each grid point, which we refer to as the Particle Remineralization and Sinking Model (PRiSM).In PRiSM, all the essential dynamics of a sizeresolved particle spectrum are included without the computational expense of explicitly simulating that spectrum.

Model formulation
We model the internal cycling of phosphorus (P) in the ocean as it is transformed between the particulate organic phosphorus (POP), dissolved organic phosphorus (DOP) and inorganic phosphate (PO 4 ) pools (Fig. 3).Only DOP and PO 4 are explicitly carried as tracers in the model -the effects of POP formation and degradation are treated implicitly as described below.The governing equations for PO 4 and DOP cycling are where A is a matrix transport operator that represents the combined effects of advection and eddy diffusion, and is derived from a data-constrained ocean circulation model (De-Vries and Primeau, 2011;DeVries et al., 2012;DeVries, 2014).The uptake of PO 4 to form organic matter (J up ) is parameterized by restoring to observed phosphate (PO 4,obs ) in the euphotic zone (taken to be the same as the mixing depth, z s ) wherever modeled PO 4 exceeds observed PO 4 using a restoring timescale of τ b = 30 days (Najjar et al., 2007), The model circulation is steady state, and does not resolve the seasonal cycle, and so we interpolate the 2009 World Ocean Atlas annual mean objectively mapped PO 4 concentrations (Garcia et al., 2010) to the model grid to obtain PO 4,obs .A fraction σ of the production is routed directly to DOP in the euphotic zone, and the remainder is routed to POP (Fig. 3).DOP is respired to PO 4 in a first-order reaction with decay rate κ, The cycling of POP is treated implicitly in the model.The rate of POP export at the base of the euphotic zone is and below the euphotic zone POP is assumed to degrade instantaneously to DOP with a vertical distribution dictated by the particle flux profile, The remaining terms in Eqs. ( 14) and ( 15) represent the P budget of the ocean as a whole.Term J sed represents the source of PO 4 in the bottom box due to the flux of POP that hits the sea floor and is regenerated.In general, we allow a fraction (f B ) of that benthic flux to be buried permanently, so that where F R is the "rain rate" of POP to the sea floor (in mmol P m −2 d −1 ) and z is the thickness of the bottom model grid  cell.Similarly, the rate of total P loss due to burial can be computed from At steady state, the rate of P burial in the sediments is matched by allochthonous inputs of P to the ocean, so J input is required to satisfy R i must be specified in order to obtain a solution to Eqs. ( 14)-( 15).The allochthonous P inputs could include dissolved and particulate P in river runoff, aeolian deposition of P in atmospheric dust, and terrestrial P from ice-rafted debris (Wallman, 2010).For simplicity, and because the spatial distribution and magnitudes of allochthonous P inputs are poorly constrained, we assume a uniform rate of P input over the entire model ocean surface.The terms in Eqs. ( 14)-( 15) involving the remineralization (J sed ) and burial (J burial ) of organic matter in the sediments, as well as the allochthonous inputs of PO 4 (J input ), are treated differently in various different model configurations, as described in the upcoming Sects.4.1-4.3.Equations (14-22) together with Eq. ( 13) constitute a complete model of the P cycle built upon a size-resolved model of particle sinking and remineralization.To calculate F (z ) we must specify the parameters of PRiSM (Eq.13), which is described in the next section.

Model validation and parameter estimation
The parameters of the particle flux model are determined from a combination of literature values and an inverse modeling procedure.Because the particle size distribution ( ) has a strong influence on POM fluxes, we include its spatial variability as inferred from the satellite-based estimates of Kostadinov et al. (2009).The value of ranges from 3.3, in locations where large particles are relatively abundant, such as the eastern equatorial Pacific, North Pacific and North Atlantic oceans, to 5.3, in subtropical gyre regions where Table 1.Parameters of PRiSM (Eq.13) and the biogeochemical models used to simulate the oceanic P cycle: CTL (control simulation without P burial or allochthonous P inputs), BUR (with P burial and allochthonous P inputs), and BUR+BAL (as BUR, but including effects of ballast-protected sinking particles).

Model configuration CTL BUR BUR+BAL Note
Particle parameters Biogeochemical parameters smaller particles are more abundant (Fig. 4a).As expected from the PRiSM particle flux profiles (Fig. 2), this results in large spatial variability in the fraction of POM reaching the deep ocean.The normalized particle flux at 1000 m below the euphotic zone varies approximately tenfold, ranging from about 0.5 in regions of large particles ( 3.5) to less than 0.05 in regions of very small particles ( 5) (Fig. 4b).Most of the particle flux reaching the deep ocean is due to the sinking of large particles.Small particles less than 200 µm in diameter contribute less than 10 % of the total particle flux at 1000 m depth (Fig. 4c).Away from the sub-tropical gyre regions, small particles generally contribute less than 5 %, and as little as 1 %, to the total particle flux at 1000 m (Fig. 4c).The spatial patterns shown in Fig. 4b and c are robust to variations in the values of the parameters controlling the particle flux profile.
The values of the other parameters controlling the particle flux profile, η, ζ , and c r , may also vary spatially due to variability in ecosystem structure and bacterial abundance.However, lacking specific information about their spatial variability, and for simplicity, we adopt spatially uniform values of these parameters here.Ideally, we would like to determine values for all of these parameters by adjusting them to obtain an optimal fit of the modeled PO 4 distribution to the observed PO 4 distribution.However, η (the exponent in the relation-ship between particle size and sinking velocity) and c r (the degradation rate of POM) have nearly identical influences on the shape of the particle flux profile (Fig. 2c and d).For this reason, we fix the value of η at 1.17, as determined from observations (Smayda, 1970), and determine ζ and c r through an optimization procedure.
The "optimal" model is determined by minimizing the volume-weighted misfit between modeled and observed PO 4 concentrations, where V is the ocean volume and dV the discretized volumes of the individual model grid boxes.When evaluating the cost function, we exclude grid points in the Japan Sea, where the model circulation is poor due to the lack of radiocarbon or CFC observational constraints on the circulation.Each model is initialized with a set of fixed parameters for PRiSM as well as the P cycling model (see Table 1).We then iteratively vary the "control" parameters c r and ζ of PRiSM with different euphotic zone depths (73 m or 115 m, corresponding to the base of the second and third model layers, respectively).For each given set of fixed parameters, O (10 2 ) model simulations are needed to find the optimal set of control parameters.This large number of model simulations is made possible by applying a Newton-Krylov method to find the equilibrium solution to the governing Eqs. ( 14)-( 15).
To focus more clearly on the effects of sinking particles on the vertical PO 4 distribution, we also investigated distributions of "regenerated" PO 4 , which is phosphate that is derived from remineralized organic matter, rather than the "preformed" PO 4 that is transported conservatively into the deep ocean from regions of incomplete surface utilization.We estimate preformed phosphate (pPO 4 ) by solving for the equilibrium distribution of PO 4 subject to the condition that all PO 4 in the euphotic zone is preformed, and there are no interior sources or sinks.Regenerated phosphate (rPO 4 ) is then computed from Preformed PO 4 is calculated from observed and modeled PO 4 distributions in the same way, using using either observed surface PO 4 or the PO 4 simulated by the model, respectively.The concentration of rPO 4 implied by the observations depends on the depth of the euphotic zone used in the calculation, which is either 73 m or 115 m.This generates a range of "observed" rPO 4 that we use as an uncertainty estimate.
As a further check on the appropriateness of the model solution, we compare model-derived particle flux profiles to observations of particle fluxes from equatorial Pacific sediment traps (Berelson, 2001).Sediment trap data are not included as a quantitative constraint on the model solution due to the large degree of scatter in the particle trap data (cf.Gehlen et al., 2006, Fig. 3) and the lack of ancillary data (e.g., surface particle size distributions) needed for a direct model/data comparison.It is also difficult to weigh the relative strengths of the PO 4 and sediment trap data appropriately as constraints on model parameters.However, we find that the equatorial Pacific sediment traps provide a valuable qualitative check on the model solution that helps to identify significant biases in the modeled deep-ocean particle flux.This is discussed in more detail in Sect. 4.

Results
Here we discuss the results from a hierarchy of model configurations designed to evaluate the ability of PRiSM to reproduce the time-averaged distribution of PO 4 .We focus in particular on depth profiles of PO 4 and regenerated PO 4 , as these are very sensitive to the particle flux profile.

Control simulation
In the control simulation (CTL), we ignore the effects of organic matter burial.In this case, any POP that reaches the sediments is instantaneously remineralized there (i.e., f b = 0).In these simulations, the J burial and J input terms are retained, but are so small that they do not affect the distribution of PO 4 or DOP, and simply serve to set the modeled total PO 4 inventory to the observed value.This is accomplished by setting the J burial term to remove PO 4 everywhere in the ocean at a rate of r g = 10 −6 yr −1 , while the J input term everywhere restores modeled PO 4 to the observed mean PO 4 at the same rate (cf.Primeau et al., 2013;Holzer et al., 2014).
Upon optimizing the PRiSM parameters c r and ζ , we find that the model achieves an excellent fit to the observed globally averaged vertical PO 4 distribution (Fig. 5a).The volume-weighted root mean square error (RMSE) is 0.14 mmol PO 4 m −3 , which yields a normalized RMSE (RMSE divided by the average PO 4 concentration) of 0.065.For comparison, the volume-weighted RMSE for PO 4 in the suite of coarse-resolution ocean biogeochemistry models considered by Duteil et al. (2012) ranged from 0.20 to 0.40, while the best-fit model considered by Kriest and Oschlies (2013) had a normalized RMSE of 0.10 for PO 4 .The model displays a very good fit to the "observed" globally averaged vertical profile of rPO 4 (Fig. 5a).The model performs worst in the lower mesopelagic zone (∼ 500-1500 m depth), where modeled rPO 4 concentrations are slightly lower than observed.Since the model circulation is constrained to match radiocarbon and CFC-11 observations (cf.DeVries and Primeau, 2011;DeVries et al., 2012), the lower-than-observed rPO 4 in this region probably indicates too little organic matter remineralization there.The model also predicts slightly higher than observed abyssal rPO 4 concentrations (below ∼ 3000 m depth), suggesting too much deep ocean remineralization.The deficiencies in the modeled vertical distribution of rPO 4 can be seen more clearly on the basin scale (Fig. 5b and c).The slightly too high abyssal rPO 4 concentrations in the model are primarily in the deep Atlantic Ocean (Fig. 5b).Too low mesopelagic rPO 4 concentrations can be traced to the Indian Ocean in the depth range 200-2000 m (Fig. 5c).
Given the overall excellent fit of the model to observed PO 4 and rPO 4 , and particularly their vertical distributions, one might expect that the model also produces a good fit to independent observations of POM settling from suspended sediment traps.However, this is not the case.In the CTL model, the optimal values of the parameters (Table 1) produces a particle sinking profile that rapidly deflects to very low values in the deep ocean (Fig. 5c).By contrast, observations from sediment traps in the equatorial Pacific Ocean (Berelson, 2001) suggest that the particle flux remains fairly constant below about 2000 m below the euphotic zone, at between 1 and 10 % of the flux at the base of the euphotic zone (Fig. 5c).Sediment traps from other locations such as the Arabian Sea, the North Atlantic, and the Southern Ocean show similar normalized particle flux values in the deep ocean (Berelson, 2001).
Thus, in the CTL model there is a conflict between the vertical attenuation of the particle flux implied by the observed PO 4 distribution, and that measured by sediment traps.The conflict is particularly severe in the deep ocean (Fig. 5c).It arises because for the model to match deep ocean PO 4 values, it must assign a fast rate of remineralization, to prevent a large particle flux into the deep ocean.This tendency is not unique to this model.In nearly every model used in Phase 2 of the Ocean Carbon Cycle Model Intercomparison Project (OCMIP-2), PO 4 concentrations in the deep ocean were significantly overestimated (Najjar et al., 2007, Fig. 8).All of these models used a power-law depth dependence of the sinking POP flux, the so-called "Martin curve" (Martin et al., 1987), and assumed a closed P budget.Since the Martin curve was derived from particle flux profiles from sediment traps, models using the Martin curve naturally overestimate remineralization in the deep ocean if burial of POM is not allowed (Kriest and Oschlies, 2013).Given the long residence times of abyssal waters, small errors in deep-ocean POM remineralization rates may accumulate into large biases in PO 4 .
One obvious solution to these biases is to allow part of the benthic particle flux to be buried permanently, rather than remineralized to PO 4 at the sea floor.This solution was examined by Kriest and Oschlies (2013), who found that explicitly modeling organic P burial was necessary in order to achieve a good fit to benthic PO 4 concentrations in a model that used the Martin curve parameterization for sinking POM.However, biases in circulation could also contribute to the deep-ocean PO 4 bias.In the following sections we explore whether adding organic matter burial can resolve the conflict between the particle flux attenuation implied by PO 4 observations and that derived from sediment traps, in a data-constrained circulation model.

Including the effects of organic matter burial
We now consider a model (BUR) in which we include the burial of POP in the sediments.We assume that the fraction of POP that is buried in sediments, f B , can be related to the "rain rate" at which POP is delivered to the sea floor, F R , following the relationship Equation ( 25) is similar to the relationship used by Burdige (2007) and Kriest and Oschlies (2013), except that here we apply the tanh function to the right-hand side to ensure that the burial efficiency f B does not exceed 1.We jointly optimized the parameters c r and ζ , along with the new parameters α, β, and R i (the rate of allochthonous P inputs needed to match POP burial) using the same procedure as for the CTL model.and ζ , are very similar for the BUR and CTL models (Table 1).The optimal values of α and β are about 0.8 and 1.45, respectively (Table 1).The optimal rate of allochthonous P inputs is about 70 Gmol P yr −1 , and about 20 % of organic matter reaching the sediments is buried there, although the uncertainty on these quantities is about 100 % (Table 2).Overall there is very little difference between the PO 4 distribution in the BUR and CTL models (compare Figs. 6 and  5).There is a slight improvement over most ocean basins in the fit of the model to the observed rPO 4 (Fig. 6b and c).The particle flux profiles in the BUR and CTL models are also very similar (Fig. 6d and Fig. 5d).The misfit between the particle flux predicted by the model and that observed from sediment traps in the deep ocean is still very evident (Fig. 6d).Thus, we conclude that the addition of organic matter burial by itself is not sufficient to resolve the conflict between the particle flux attenuation implied by the PO 4 observations, and that implied by the sediment trap observations.
The reason that burial alone cannot reconcile the nutrient distributions with sediment trap data is that the burial rate is proportional to the benthic flux of POM, and thus decreases rapidly with depth.Burial of P in deep sediments permits a larger particle flux to reach the deep ocean without creating a surplus of PO 4 .However, because of the rapid particle flux attenuation, this would require even more P removal from shallower depths, creating a low PO 4 bias there.To fit the PO 4 globally, the model therefore must maintain a low rate of PO 4 burial overall.This trade-off between PO 4 biases in the deep and mid-depth water column would be less stringent if the flux of POM did not decrease so rapidly with depth.This suggests that one solution to the apparent discrepancy between the particle flux data and the PO 4 distribution is for a fraction of the sinking flux of particulate P to be relatively recalcitrant.This would allow its flux to decrease less rapidly downward, so that burial from the deep sea could be achieved without slowing PO 4 regeneration too much in the thermocline.The need for a component of POM that resists degradation has been proposed as an explanation for the constancy of the deep particle fluxes, with protection of organic matter from bacterial degradation by inclusion in ballast mineral assemblages as a specific mechanism for it (e.g., Armstrong et al., 2002;Francois et al., 2002;Klaas and Archer, 2002).We test this hypothesis in the model as described in the next section.

Including the effects of ballast-protected organic matter
Here we separate the flux of sinking POM into two pools with different time scales of degradation to investigate the effects of a slowly degrading P pool on the total particle fluxes and PO 4 distributions.As a basis for this separation, we adopt the hypothesis that mineral ballast acts to protect some organic matter from bacterial degradation (Armstrong et al., 2002;Francois et al., 2002;Klaas and Archer, 2002).
Other mechanisms for creating a more slowly degraded component of particulate P are also possible, however.In particular, the recent discovery of significant concentrations of polyphosphates in organic matter in both the water column and sediments (Diaz et al., 2008) will be discussed below as an alternative interpretation of the model results.
The exact mechanism of ballast-mineral protection is not completely understood, but laboratory experiments suggest that ballast minerals may be scavenged onto particle aggregates during the sinking process (Passow and De La Rocha, 2006).This mechanism was explored in the model of Gehlen et al. (2006), who found that the combined effects of particle aggregation and mineral ballasting resulted in large particle fluxes to the deep sea.Here we do not explicitly simulate scavenging of ballast minerals onto sinking organic matter particles, but make the simplifying assumption that the fraction of POM that is ballast protected is proportional to the flux of ballast minerals out of the euphotic zone (Armstrong et al., 2002).In this case we can express the total flux of POP as the sum of an unprotected component and a ballastprotected component, where f P is the fraction of the POP produced in the euphotic zone that is routed to the ballast-protected POP pool, F U (z ) is the particle flux profile for unprotected POP, and F P (z ) is the particle flux profile for ballast-protected POP.We assume that f P is proportional to the "ballast ratio", R B , which is the mass ratio of the sinking flux of ballast minerals to the sinking flux of organic carbon at the base of the euphotic zone, Our P cycle model does not simulate ballast mineral or organic carbon fluxes, and so we use values of R B from the Geophysical Fluid Dynamics Laboratory Earth System Model (GFDL-ESM) to calculate f P (see Appendix A and Fig. A1).Following Armstrong et al. (2002), the value of ρ is assumed to be 0.05.With these values for ρ and R B , the value of f P varies between 0.04 and 0.25, with a mean value of 0.11.We model the sinking of unprotected and protected POP separately.For protected POP, we use the same parameters as for unprotected POP, except that rather than solving for c r and ζ as part of the inversion, we fix these parameters at values that give reasonable ballast mineral flux profiles.According to Armstrong et al. (2002), F P ≈ 0.4 in the deep ocean.Assuming ζ = 2.28 (Mullin et al., 1966), and for an average value of 4.2, a value of c r = (365 d) −1 gives a value of F P = 0.4 at about 5000 m below the base of the euphotic zone.These then are the parameter values we specify for sinking POP that is ballast protected (Table 1).Because the protected POP is protected from bacterial degradation, we expect it to be buried with a much greater efficiency.Therefore, we use different values of α and β in Eq. ( 25) for protected and unprotected POP.
The resulting model that includes both organic matter burial and ballast mineral effects (BUR+BAL) has seven unknown parameters: c r , ζ , α, β, R i , and α P and β P (for protected POP).We jointly optimized these parameters using the same procedure as for the CTL and BUR models.The optimal value of c r is about (21 d) −1 , lower than the ∼ (30 d) −1 in the CTL and BUR models (Table 1).This degradation rate for unprotected POP is nearly 20 times faster than that assumed for ballast-protected POP.On average, 95 % of ballast-protected POP reaching the sediments is buried there, while only 7 % of unprotected POP reaching the sediments is buried (Table 2).The absolute rates are discussed in the next section.
The BUR+BAL model shows improvement compared to the CTL and BUR models in the overall fit to observed PO 4 .The globally averaged normalized RMSE is 0.060 for the BUR+BAL model, and 0.065 for the CTL and BUR models (Fig. 7).There is also significant improvement in the fit to regenerated PO 4 on the global and basin scales, with fits improving by up to 15 % compared to the BUR model.Most significantly, the particle fluxes in the BUR+BAL model are consistent with the flux estimates from sediment traps in the deep sea (Fig. 7d).The sinking flux of unprotected POP deflects to even lower values in the deep ocean than in the CTL and BUR models.However, about 10 % of the sinking POP is protected by ballast minerals, allowing the total POP sinking flux to reach the observed values of about 0.01-0.1 of the flux at the base of the euphotic zone in the deep ocean (Fig. 7d).We therefore conclude that in order to simulate both POM flux to the deep ocean and the remineralization of PO 4 in the deep ocean correctly, both organic matter burial as well as ballast mineral protection must be modeled.
The magnitude of organic P production is relatively constant among the model configurations (Table 2).Total organic P production ranges from 13.8 ± 2.3 Tmol P yr −1 in the BUR model to 14.4 ± 2.8 Tmol P yr −1 in the BUR+BAL model (Table 2).For comparison, Dunne et al. (2007) used satellite chlorophyll observations and empirical models to estimate that 9.6 ± 3.6 Pg C yr −1 is exported out of the euphotic zone as particles.Assuming a C : P ratio of 106 : 1 in fresh organic matter (Anderson, 1995) and that 80 % of organic matter production is exported as sinking particles (Hansell et al., 1997), yields a rate of organic P production of 9.4±3.5 Tmol P yr −1 .This is on average smaller than the rate of P production derived here, but the estimates agree within their relatively large uncertainties.
In contrast to total POM production, the model configurations yield substantial differences in the latitudinal patterns of P fluxes within the ocean, and in the total input/output budget of P.These are discussed in the next two sections.

Implications for P cycling
Here we consider the influence of two characteristics of the surface particle distribution -the surface particle size distribution and the ballast ratio -on the internal cycling of P. We find that both of these effects lead to a reduction in the latitudinal variation of export and subsequent deep remineralization.
The largest difference in organic P production among the model configurations is caused by the inclusion of ballastmineral protection in the model BUR+BAL.Relative to the models without ballast-mineral protection, production is reduced in the Southern Ocean and increased in the tropical and sub-tropical oceans (Fig. 8a).This is due to the effect of ballast minerals on the remineralization profile of sinking organic matter.In the Southern Ocean, the ballast ratio is relatively high due to high production rates of biogenic silica associated with diatom-dominated communities.Because the ballast ratio is high, particles sink deeper on average before remineralizing, and therefore the supply of remineralized nutrients to surface waters, which can fuel new production, is reduced.The opposite effect occurs in the tropical and sub-tropical ocean, where the ballast ratio is low (Appendix Fig. A1).
The surface particle size distribution may also have a significant effect on organic matter production.To evaluate this, we re-ran each of the models in the CTL and BUR + BAL configurations using a spatially uniform surface particle size distribution, with = 4.2, in place of the spatially variable used in the standard configuration (see Fig. 4).The results show a significant effect of the particle size distribution on organic matter production rates.With a spatially variable surface particle size distribution, organic P production rates are reduced in regions of high productivity, and en- hanced in regions of low productivity in both the CTL and BUR+BAL models (Fig. 8b).This effect occurs because regions of high productivity tend to be dominated by larger particle assemblages, while regions of low productivity are characterized by smaller particles (Kostadinov et al., 2009).Larger particles on average sink deeper before remineralizing than smaller particles, and therefore the supply of regenerated nutrients to the euphotic zone is reduced when large particles are produced, reducing new production.These results suggest that models using a spatially uniform particle sinking speed, or spatially uniform particle remineralization profile, will overestimate production in high-productivity areas such as coastal upwelling regions, and underestimate productivity in low-productivity regions such as the oligotrophic sub-tropical gyres.
The remineralization of PO 4 in the sediments is controlled by the rain rate of organic P to sediments and by the POP burial efficiency.However, the rate of benthic remineralization in the deep ocean does not strongly depend on whether organic P burial or ballast effects are explicitly modeled.This is because in each model the parameters controlling the particle sinking profile and the burial efficiency are adjusted to achieve the best possible fit to observed PO 4 , and the deep ocean PO 4 concentrations provide a strong constraint on the rate of benthic remineralization.Benthic remineralization below 2000 m in the CTL and BUR models are nearly identical, at 82 ± 9 and 79 ± 8 Tmol P yr −1 , respectively (Table 2).Deep-ocean benthic remineralization in the BUR+BAL model is slightly lower, at 59 ± 35 Tmol P yr −1 .The main difference in deep-ocean benthic remineralization among the model configurations tested here results from adding mineral ballast effects (model BUR+BAL), which changes the spatial structure of benthic remineralization (Fig. 8c).Because ballast-protected POP is buried more efficiently in deep-ocean sediments, benthic remineralization in the BUR+BAL model is reduced in areas with high ballast ratios, compared to the BUR model (Fig. 8c).We also calculated benthic remineralization in the CTL and BUR+BAL models with a uniform surface particle size distribution with = 4.2.Compared to the case in which varies spatially, benthic remineralization rates are decreased nearly everywhere in both models (Fig. 8d).This is because the POP flux to the deep sea tends to be dominated by large particles (cf.Fig. 4).Imposing a uniform surface particle size distribution tends to reduce particle sizes in high-productivity regions, and ultimately less POP is delivered to the sea floor, resulting in lower benthic remineralization rates.This suggests that models that do not consider spatially variable particle sizes and particle sinking rates will underestimate benthic remineralization rates in the deep ocean, particularly under high-productivity regions.

Implications for the P budget: a more dynamic marine P cycle?
The model fit to observed PO 4 , and to observed particle flux profiles from sediment traps, is best when both organic mat-ter burial and mineral ballast effects are included (model BUR+BAL).In that case, the optimal rate of P burial in the sediments is 698±137 Gmol P yr −1 .If these P burial rates are correct, they would imply a much more active marine P cycle than previously thought.The oceanic residence time of P derived from the BUR+BAL model is about 3400-5200 yr.This suggests that the marine P cycle may be as dynamic as the marine N cycle, since marine fixed N has a mean residence time of about 3500-5000 yr (Eugster and Gruber, 2013;De-Vries et al., 2013).
The large burial flux of organic P in the BUR+BAL model is driven almost exclusively by the burial of ballast-protected POP.The total burial of ballast-protected POP is 684 ± 150 Gmol P yr −1 , while the total burial of unprotected POP is nearly negligible at only 14±20 Gmol P yr −1 (Table 2).This difference in burial rates can be traced to the much higher burial efficiency of ballast-protected POP.The burial efficiency of ballast-protected POP is about 95 %, while the burial efficiency of unprotected POP is only about 5 % (Table 2).In the BUR+BAL model, we find that burial efficiencies generally increase with rain rate for unprotected POP (Fig. 9a, dashed red curve), but that the burial efficiency of ballast-protected POP is relatively constant with rain rate (Fig. 9a, dashed blue curve).This leads to very different depth dependencies of burial for the unprotected and ballastprotected POP fractions (Fig. 9b).Unprotected POP is preferentially buried in shallow sediments, where POP fluxes are relatively high, while the burial rate of ballast-protected POP decreases only slightly with depth due to the decrease in particle flux with depth.
A difference in burial efficiency of ballast-protected vs. unprotected POP has not to our knowledge been measured, and therefore cannot be confirmed.However, the finding that a substantial fraction of recalcitrant P is buried in the deep ocean is consistent with findings that the C : P ratio of organic matter burial is lowest in low-sedimentation rate pelagic environments (Ingall and Van Cappellen, 1990).Moreover, we can compare model fluxes to observations of sedimentary rain rate and burial of organic carbon, which were compiled by Kriest and Oschlies (2013) (their Table 4), from about a dozen locations ranging from shelf sediments to abyssal plain sediments.We converted the C rain rate to P rain rate using a C : P ratio of 110 : 1 (Wallman, 2010), and converted C burial rates to P burial rates using ratios given by Wallman (2010) (his Fig. 3) for different marine environments: the C : P of organic matter burial is 32 : 1 in shelf sediments (taken here to be less than 200 m in depth), 23 : 1 in slope regions (taken here to be greater than 200 m and less than 2000 m in depth), and 15 : 1 in abyssal sediments (taken here to be greater than 2000 m in depth).The results show that the burial efficiency of P can vary widely for equal rain rates (Fig. 9a).This variability occurs approximately within the limits of the burial efficiencies derived for unprotected and ballast-protected POP in the model, and Previous estimates of organic matter burial in abyssal sediments vary widely.On the one hand, the empirical formulations of Dunne et al. (2007) yield a carbon burial flux of 0.012 ± 0.02 Pg C yr −1 which, using a C : P ratio of 15 : 1 for organic matter burial in abyssal sediments (Wallman, 2010), yield a P burial of 67 ± 111 Tmol P yr −1 .On the other hand, the empirical formulations of Muller-Karger et al. ( 2004) yield a burial flux of 0.09 Pg C yr −1 below 2000 m, which yields a P burial of 500 Tmol P yr −1 .Finally, observations of the P content of marine sediments suggest that only about 80 Gmol P yr −1 accumulates in deep-sea sediments (Baturin, 2007;Wallman, 2010).

Conclusions and caveats
We present a model of the ocean P cycle based on a sizeresolved and spatially variable model of particle fluxes (PRiSM) embedded in a data-constrained ocean circulation model.From a hierarchy of model configurations, we find that the size distribution of particles exiting the surface ocean, and the ballasting of exported organic matter are important controls on P fluxes within the ocean and its longterm burial in the deep ocean and sediments.The strength of these results rests on the use of a mechanistic formulation of particle dynamics, and an ocean circulation model that is able to match tracers of ocean ventilation rates.Still, each of these components contains simplifications that could influence the results.
First, the ocean circulation model lacks a seasonal cycle.Our diagnostic approach to export fluxes based on nutrient restoring should provide a good estimate of the export fluxes from the upper ocean, and the integrated rate indeed matches other empirical estimates.However, any covariation between particle size distributions, ballast content, and export flux are not represented.It is unclear what the effect of such seasonal and higher frequency covariations would be.
Second, the particle model used to drive the global P cycle simulations makes several simplifications about particle dynamics and the associated biological rates.The use of a single sinking speed for each particle size, the neglect of coagulation and fragmentation below the turbulent boundary layer, and of environmental effects on the intrinsic (per mass) rates of particle decomposition, are all simplifications that need further investigation.Given the relative homogeneity and quiescence of the water column below 2000 m, it seems unlikely that any of these simplifications could reconcile the apparent conflict between the sediment trap and nutrient data at those depths, so as to obviate the need for a dynamic P budget.
The most important caveats then, concern the factors that give rise to the high P burial rates in the deep sea, implied by the BUR+BAL model (> 500 Gmol P yr −1 ).Here we have used a simple formulation for ballast protection based on the ratio of the sinking flux of ballast minerals to organic carbon in the euphotic zone.An alternative origin of less degradable P in organic particles is non-reactive detrital P (Paytan and McLaughlin, 2007), such as the polyphosphates observed in organic matter in both water and surface sediment material (Diaz et al., 2008).Polyphosphates are produced primarily by diatoms, so that their contribution to total organic P export may have a similar spatial pattern to that of ballast.Moreover, the proportion of polyphosphates (7-8 % of organic P) is similar to the fraction of ballast-protected carbon estimated from observations and model simulations, and used in our calculations (about 10 %, see Appendix A).Given these similarities, and the large uncertainties associated with both mechanisms, we view either of them as providing a plausible interpretation for the BUR+BAL model results.
While the ballast-protected organic matter formulation appears to match well with observations of deep-sea particle fluxes in the equatorial Pacific, there are several sources of uncertainty that we have not accounted for.First, we have assumed a uniform proportionality constant (ρ = 0.05) between the fraction of ballast-protected POP and the ballast ratio in the euphotic zone.Armstrong et al. (2002) found a mean value of ρ = 0.05 for the equatorial Pacific, but also report values of ρ ranging from 0.027 to 0.065 in the Southern Ocean.Since the flux of ballast-protected POP to the deep ocean scales linearly with ρ, a factor of two uncertainty in ρ should lead to a factor of two uncertainty in the burial rate of POP in the deep ocean.
Second, we have assumed a degradation rate, c r , of (365 d) −1 for ballast-protected POP, which for a typical value of (4.2) matches the fraction of ballast-protected POM reaching the deep ocean (0.4) estimated by Armstrong et al. (2002).However, given the spatial variability in , the fraction of ballast-protected POP reaching the deep sea in the model ranges from about 0.1 to 0.7.Thus, uncertainty in c r for ballast-protected POP probably contributes to an additional factor of two uncertainty in the deep-sea POP flux and burial rate.
Third, the sediment trap data and the estimates of ρ are based on C fluxes to the deep ocean.However, measurements of particle C : P from the European continental margin indicate that the ratio of C : P in particles appears to increase with depth, when one considers solubilization of particles within the sediment traps (Antia, 2005).If this relationship holds globally, then we would expect the flux of POP to the deep ocean to be reduced by a factor of two relative to the estimates here, reducing P burial accordingly.
Another possibility for the large discrepancy between our model-based rate and the sediment-based rate of organic P burial is that the sedimentary records do not adequately sample regions with large fluxes of ballast-protected POP to the deep ocean.In the BUR+BAL model, approximately 50 % of ballast-protected POP burial in the deep ocean occurs in the Southern Ocean (south of 30 • S), a region that is very poorly sampled (cf.Palastanga et al., 2011, Fig. 8b).If polyphosphates play a key role in P burial, we would also expect large burial fluxes in the Southern Ocean, where diatom production is high.For this reason, the estimates based on sedimentary data may significantly underpredict burial of organic P in the deep ocean.To ultimately reconcile the model-predicted deep-ocean P burial rates and those derived from sediment data will require much more high-quality deep-ocean sediment trap data.With sufficient spatial coverage, it should be possible to constrain the parameters of the ballast-protected sinking POP fraction better, such as c r and ρ, rather than specifying them based on limited data, as we have done here.This would allow a more accurate determination of POP delivery to the deep ocean, and the deep ocean PO 4 data would then be better able to constrain the POP burial in sediments.Until these questions can be resolved empirically, we regard the high rates of P burial implied by the model hierarchy as intriguing but somewhat tentative.

Figure 1 .
Figure 1.(a)The size distribution, n, of particles as a function of depth and particle size class as predicted by Eq. (8) for the surface size range D S = 20 µm to D L = 2000 µm.Black x's mark the upper size limit of particles at each depth from Eq. (11).(b) The mass of particles within each size class with depth.(c) The normalized integrated particle mass and particle flux with depth using the size distribution from (a).(d) The mass-weighted particle sinking velocity with depth.All calculations used the following parameter values: c w = 2.2 × 10 5 m 1−η d −1(Kriest and Oschlies, 2008), = 4.2(Kostadinov et al., 2009), η = 1.17(Smayda, 1970), ζ = 2.28(Mullin et al., 1966), c r = 0.03 d −1(Kriest and Oschlies, 2008).

Figure 2 .
Figure 2. The sensitivity of the normalized particle flux predicted by the Particle Remineralization and Sinking Model (PRiSM, Eq. 13) to the four parameters controlling the particle flux: (a) , the exponent in the relationship between particle size and particle size distribution, (b) ζ , the exponent in the relationship between particle mass and particle size, (c) η, the exponent in the relationship between particle sinking velocity and particle size, and (d) c r , the degradation rate of sinking particles.Solid curves use the parameters from Fig. 1.

Figure 4 .
Figure 4. (a) Spatial variability in the exponent for the surface particle size distribution used in PRiSM.Lower values of indicate larger particles, and higher values of indicate smaller particles.(b) Spatial variability in the particle flux at 1000 m depth resulting from variability in the surface particle size distribution.(c) The fraction of the flux at 1000 m due to small particles (less than 200 µm in diameter).For (b) and (c) we used the same parameters of PRiSM as in Fig. 1.

Figure 5 .
Figure 5. (a) Globally averaged depth profile of total PO 4 (red curve) and regenerated PO 4 (rPO 4 , blue curve), for the CTL model and the observations (black x + error bars).(b) Modeled (curves plus shading indicating uncertainty) and observed (black x + error bars) rPO 4 averaged over the Pacific (red) and Atlantic (blue) oceans.(c) Same as (b) for the Indian (red) and Southern (blue) oceans.(d) Modeled particle flux profile (red curve plus shading indicating uncertainty) for = 4.2, and observed particle fluxes from sediments traps in the equatorial Pacific (symbols).Printed on (ac) is the normalized root mean squared error (RMSE divided by the average PO 4 or rPO 4 concentration) for the CTL model for the region and data type displayed.

Figure 6 .
Figure 6.As Fig. 5, but for the BUR model.

Figure 7 .
Figure7.As Fig.5, but for BUR + BAL model.In (d), the average sinking fluxes of unprotected POP (dashed red curve) and ballastprotected POP (dashed blue curve) are shown separately.The uncertainty on the total POP sinking flux (magenta curve + shading) reflects uncertainty on the parameters c r and ζ of PRiSM, and spatial variability in the ratio of ballast mineral to organic carbon flux at the base of the euphotic zone (see Fig.A1).

Figure 8 .
Figure 8.(a) New production by latitude for the CTL, BUR, and BUR+BAL models (averaged over the four different model configurations).(b) As (a), but comparing new production rates for the CTL and BUR+BAL models using the standard spatially varying values, and a spatially uniform value of = 4.2.(c) As (a), but showing benthic remineralization below 2000 m.(d) As (b), but for benthic remineralization below 2000 m.

Figure 9 .
Figure9.(a) POP burial rate vs. the rain rate (rate of delivery of POP to sediments) for the BUR and BUR+BAL models (averaged over the four different model configurations).For the BUR+BAL model, the relationship between burial and rain rate differs substantially for ballast-protected (blue dashed curve) and unprotected POP (red dashed curve).Observations (marked with a *) are taken from Table4ofKriest and Oschlies (2013) using C : P ratios from Fig.3ofWallman (2010).(b) POP burial rate as a function of depth in the BUR and BUR+BAL models (averaged over the four different model configurations).

Table 2 .
Rate of P production, benthic remineralization, and burial from the three different P cycle models considered here.