On the treatment of particulate organic matter sinking in large-scale models of marine biogeochemical cycles

. Various functions have been suggested and applied to represent the sedimentation and remineralisation of particulate organic matter (POM) in numerical ocean models. Here we investigate some representations commonly used in large-scale biogeochemical models: a constant sinking speed, a sinking speed increasing with depth, a spectrum of particles with different size and different size-dependent sinking velocities, and a model that assumes a power law particle size distribution everywhere in the water column. The analysis is carried out for an idealised one-dimensional water column, under stationary boundary conditions for surface POM. It focuses on the intrinsic assumptions of the respec-tive sedimentation function and their effect on POM mass, mass ﬂux, and remineralisation proﬁles. A sinking speed appropriate for simulations exceeding a the sedimentation proﬁle is not with observed is we present an analytic integral over depth and size that can explain regional variations of remineralisation length scales in response to regional patterns in


Introduction
The sinking and remineralisation of particulate organic matter (POM) in the ocean creates vertical gradients in dissolved inorganic tracers, and affects the air-sea gas exchange of CO 2 and O 2 between the ocean and the atmosphere. A synoptic and coherent view of the ocean's distribution of biogeochemical tracers and their exchange with the atmosphere is usually achieved by simulations of basin-wide or global biogeochemical circulation models.
Production of POM is confined to the surface layer with light levels sufficient for photosynthesis. Models that calculate the flux of POM out of this surface layer account for POM sedimentation in different ways: early models parameterised an increase in POM sinking speed with depth by applying the empirically derived algorithm of Martin et al. (1987;e.g. Najjar et al., 1992;Maier-Reimer, 1993). A three-dimensional application of the model by Fasham et al. (1990) employed a constant detritus sinking speed in the upper 123 m and an instant sedimentation and remineralisation profile according to Martin et al. (1987) below (Sarmiento et al., 1993). Recently, global models have been presented that either explicitly prescribe an increase of POM sinking speed with depth (Schmittner et al., 2005), or partition POM into two different size classes with different constant sinking speeds (e.g., Gregg et al., 2003). Other approaches have suggested an effect of mineral ballast on the remineralisation length scale (Armstrong et al., 2002;Francois et al., 2002;Klaas and Archer, 2002;Gehlen et al., 2006).
The choice of constant sinking velocities may be justified by observations of individual particles (e.g. Smayda, 1970;Kriest, 2002, and citations therein). We must, however, distinguish between the properties of individual particles and the property of an aggregated POM compartment as commonly simulated in numerical models: POM (here: phytoplankton and detritus) consists of many different particles, Published by Copernicus Publications on behalf of the European Geosciences Union. which may vary in many aspects: their constituents (e.g., calcifiers or diatoms vs. flagellates), age, origin, etc. Armstrong et al. (2002) have ascribed differences in POM sinking to the variation in particle composition, and Boyd and Trull (2007) present a detailed overview over the different models of ballast-associated export and their rationale. Another important aspect, on which the present work focuses, is (phyto) plankton particle size, which, in the ocean, ranges from ≈0.2−1000 µm.
Generally, we can expect the sinking speed of an individual particle to increase approximately proportional to its diameter (Smayda, 1970). What effect does this have on the sinking speed of total POM? -Consider an ensemble of particles of different size at a given depth, z 0 , that starts its journey downwards: if individual particle sinking speed increases with its size, but remineralisation rate is constant, we can expect the average POM size and sinking speed to increase with depth, because predominantly the particles with large size and high sinking speed reach the deep ocean; the small (i.e. slow) ones will remineralise in the upper layers.
Empirical and theoretical studies indeed suggest such an increase of POM sinking speed with depth: Banse (1990Banse ( , 1994 proposed an exponential function for the description of mass flux with depth, but also suggested that the exponent (i.e., remineralisation rate over sinking velocity) of this function should be depth dependent -however, he did not comment on the exact form of the depth dependence. Lutz et al. (2002) accounted for different remineralisation length scales by fitting a sum of two exponential functions to observations of sedimentation. Martin et al. (1987) found that profiles of sedimentation collected with sediment traps could best be fitted by a power law (hereafter Martin curve), F (z)=F (0) (z/z 0 ) −0.858 , which either implies a decrease of remineralisation rate with depth (r ∝ 1/z), or an increase of (average POM) sinking speed with depth (w ∝ z; see below for derivation). Berelson (2002) analysed arrays of sediment traps and showed that the sinking speed of POM increases with depth.
Given the sensitivity of biogeochemical model results to the parameterisation of sinking speed (Heinze et al., 2003;Howard et al., 2006;Gehlen et al., 2006), in this paper we investigate two functions commonly applied (the Martin function and constant sinking speed) and their effect on the representation of POM profiles, sedimentation, and remineralisation. We do this by means of analytic solutions for the above mentioned functions, assuming stationary and seasonally varying boundary conditions for POM sinking out of the surface layer.
We do not attempt to describe in detail a particular group of particles, such as zooplankton fecal pellets or phytoplankton aggregates, but instead focus on the relatively simple, yet efficient parameterisations commonly applied in threedimensional large-scale marine biogeochemical models. In doing so, we consider sinking organic matter to be a mixture of (unspecified) particles with certain characteristics. In particular, we contrast two simple parameterisations (constant POM sinking speed, and sinking speed increasing with depth) with a model that simulates a discrete POM size spectrum, in which all size classes have a size-dependent but depth-independent sinking speed. We finally examine if, and to what extent, we can predict deep sedimentation from the size distribution of POM in the surface layer.

Model setup and results
For all of the following representations, we consider a water column of 4000 m depth below the base of the euphotic zone (located at depth z 0 ), which is not affected by horizontal processes, or by vertical mixing, with the z-axis pointing downwards. z ′ is the depth referenced to z 0 (z 0 + z ′ is the total distance from the sea surface). For the sake of simplicity, we first consider constant upper boundary conditions of POM mass: M(z 0 )=M 0 =1 mmol N m −3 (see Table 1).
For the first two models (constant sinking speed and sinking speed varying linearly with depth) we set the sinking speed of POM at z 0 to w 0 =3.52 m d −1 which is in the range of numerical models (e.g., Doney et al., 1996;Oschlies and Garçon, 1999, ; see also Table 1). (The value of 3.52 m d −1 corresponds to the average POM sinking speed of the model with 198 size classes described below.) This results in a nitrogen export out of the euphotic zone into the model domain of 3.52 mmol N m −2 d −1 which is about 2-10 times higher than global mean new production (range of observational estimates and box models: 0.27-1.53 mmol N m −2 d −1 ; Oschlies, 2001) and is supposed to represent highly productive regimes. We further assume that the remineralisation rate r is constant: r=0.0302 day −1 . The choice of this value is explained below; it is in the range of remineralisation rates applied in other biogeochemical models. The sensitivity of the models to this value, and the rationale for choosing sizeindependent remineralisation, will be explored in the discussion section.
The third model resolves a discrete POM size spectrum of 198 classes. We first define the particle characteristics (size range and the parameters b 1 , w 1 , ζ and η; see below for definition). We then define the spectral exponent of an assumed power law size-distribution of POM at the upper model boundary, ǫ 0 =ζ + η + 1.01 (see Table 1 and below for the choice of parameters). This results in an average POM sinking speed at z 0 of 3.52 m d −1 .
The exponent ǫ 0 of oceanic particle size spectra, after conversion to the distribution defined in Eq. A1, in the ocean ranges from ≈2.2 (Marañón et al., 2004) to ≈5.2 (Cavender- Bares et al., 2001), but many observations suggest that sizespectra are more or less "flat", i.e., mass is distributed evenly among logarithmically increasing size classes. In our case, assuming that the exponent that relates particle mass to its diameter (ζ ) is set to 2.28 (Mullin et al., 1966) a "flat" mass distribution would be characterised by ǫ 0 =3.28, whereas a "flat" volume distribution would yield ǫ 0 =4. For our analysis, we have chosen ǫ 0 =ζ + η + 1.01 ± 0.5=4.46 ± 0.5, where η is the exponent that relates a particle's sinking speed to its diameter, and is set to 1.17 (Smayda, 1970). The value of the standard run (ǫ 0 =4.46) was chosen because it allows the direct evaluation of the incomplete gamma function (i.e., a F >0, see Eq. A16), while still being in the range of observed spectral exponents; two experiments (ǫ 0 =4.46 ± 0.5) explore the sensitivity of the function to alterations of the exponent.
The fourth model is a size-continuous model that applies the same POM power law boundary condition as the sizediscrete model. The integration of a (continuous) size-range results in a slightly higher (3.73 m d −1 ) average POM sinking speed at z 0 .

Constant POM sinking speed (w=const)
First, consider one class of particulate organic matter of mass M that consists of particles of uniform size, having the same sinking velocity w and remineralisation rate r, which do not change with depth or time. This assumption implies that, as the particles remineralise, they do not get smaller or less dense. The time rate of change is then The mass concentration at equilibrium (∂M/∂t=0) at any depth z ′ (referenced to z 0 ) is given by likewise the POM mass flux in equilibrium is given by On a logarithmic scale, the distribution of mass and sedimentation with depth naturally turns out to be a straight line.
For the given parameters the function implies an e-folding length scale for mass and sedimentation of 142 m, and causes the POM concentration and flux to decrease by about two orders of magnitude within the upper 300 m (Fig. 1).

The Martin curve (w ∝ z)
Now we assume that POM sinking speed increases linearly with depth, similarly to Schmittner et al. (2005): given the function w(z)=a z + b, with b=0 the time rate of change for POM is In equilibrium (∂M/∂t=0) and, likewise, for the flux With z 0 + z ′ =z and r/a=0.858, function 6 represents the parameterisation of Martin et al. (1987). For our parameters w 0 =a z 0 =3.52 m d −1 and a reference depth of z 0 =100 m this results in a=0.0352[d −1 ] and, for r/a=0.858, yields a remineralisation rate of r≈0.0302 [d −1 ]. Note that for b>0 (as, e.g., in Schmittner et al., 2005, who used b≈4.54 m d −1 and a≈0.009 d −1 below 50 m), b/a will have to be added both to the enumerator and the denominator of the base of functions 5 and 6.
The decline of POM with depth in the upper 300 m is quite strong, but then quickly ceases (Fig. 1). POM at 1000 m depth is orders of magnitude higher than in the model with constant sinking speed. On the other hand, the vertically increasing sinking speed with depth causes a much slower decrease of mass flux with depth.
In this function the increase of average sinking speed with depth is not based on mechanistic rules, but deduced from observed profiles. One possible reason for this increase of sinking speed with depth can be found in an increase in average particle size, as investigated in the following paragraph.

A spectrum of 198 discrete POM size classes
We now consider a POM size spectrum (size measured as equivalent spherical diameter), from some lower boundary d 1 to an upper boundary d L , which is divided into 198 size classes of equal width, d. In our example, we consider a size range of 20−2000 µm with d=10 µm. The entire mass of POM, M is given by where M i is the mass in a class i. The time rate of change in each size class i is given by where M i =M i (t, z). The remineralisation rate is assumed to be independent of size. We assume that the sinking speed w i of particles of each size class i is determined by the size of its lower boundary, d i : is the sinking speed of the smallest particle, and η determines the dependence of the particle's sinking speed on its diameter (Smayda, 1970, see Table 1 for parameters). We assume that the coefficients w 1 and η of this function do not change with depth or time. This assumption implies that the size of individual particles does not decrease-in terms of diameter or weight-due to remineralisation. Instead all mass losses in a size class are concentrated in a few selected particles that disintegrate immediately. The analytic solution over z for each individual size class is then the same as for the one-size-class model: i.e., Analogously, total sedimentation is given by The average POM sinking speed at at any depth is then given by Biogeosciences, 5, 55-72, 2008 www.biogeosciences.net/5/55/2008/ We now assume that particles at the upper model boundary are distributed according to a power law, and that the coefficients of this distribution, ǫ 0 and A 0 do not change with time. Thus, M 0,i =const. is defined by: with C=b 1 /d ζ 1 , b 1 being the biomass of the smallest particle. ζ is the exponent that determines the relationship between a particle's diameter and its mass, and is set to 2.28 (Mullin et al., 1966). With total particle mass M 0 and the parameters given in Table 1, this results in total flux F 0 =3.52 mmol N m −2 d −1 (see also Fig. 2, upper black line, for the distribution of particle mass).
As already outlined in the introduction, large particles will travel further downwards, while particles with small size dissolve predominantly in the upper layers, as illustrated in Fig. 2. According to the model's prerequisites the model starts from a particle size distribution that is linear on a loglog scale (upper black line in Fig. 2). Because especially the small, slow particles are remineralised when they travel through the water column, the size distribution becomes unimodal with increasing depth (e.g., red line for 100 m depth below the euphotic zone, in Fig. 2). The diameter of maximum mass increases with depth.
As a consequence of the small particles becoming less abundant with depth, the average sinking speed of POM increases with depth ( Fig. 1), but not linearly (as for the Martin curve). POM mass decreases quickly in the upper 300 m. The sedimentation profile is similar to the POM profile, i.e., the increase of POM sinking speed with depth does not compensate the decrease of POM mass with depth. As neither the sedimentation nor POM profiles are straight lines on the log plot ( Fig. 1), they cannot be represented by an exponential function. Instead, the appropriate algorithm would rather be a sum of exponential functions of z, each term with its own coefficients, as in Eqs. (11) and (10). Summarising, accounting for the development of particle size distribution with depth, that arises solely from differential sinking and constant remineralisation, has a very strong effect on simulated POM concentration and its size distribution, as well as on sedimentation and sinking speed. To some extent, the resulting mass and mass flux profiles resemble the empirical Martin curve.

A continuous size spectrum of POM
The size-discrete model presented above makes an implicit assumption about the particle size distribution within the size classes. It further assumes that all particles within the size classes can be characterised by a single sinking speed. A continuous size range and analytic integration over the entire size range can provide further insight if, and how much, the discretisation of the particle length scale affects the model solution.
Again we assume that particles at the upper model boundary are distributed according to a power law, this time on an infinitely fine size grid, with d→0 for the entire size range from d 1 to d L (see Appendix A, Eq. A10). The model applies the same size-dependency of sinking speed as the discrete model. Given M 0 and the parameters in Table 1, F 0 =3.73 mmol N m −2 d −1 . This is slightly higher than the input flux of the discrete spectrum, because now the particles' sinking speed increases continuously with size.
Despite the different setup, the analytic solution ( Fig. 1, solid lines) cannot be distinguished from the size-discrete model. This indicates that the fine discretisation of the length scale presented in the previous section has only little influence on tracer distributions and fluxes.
However, the advantage of the continuous solution is not only an exact representation of deep particle mass and fluxmore importantly, we get an idea about how deep sedimentation might depend on parameters of the surface POM size distribution. The deep mass flux F (z ′ ) in the size-continuous model for ǫ>ζ + η + 1 is a F >0 where a F =(ǫ 0 −ζ −η−1)/η, X=r z ′ /w 1 , x=r z ′ /w L and γ (a F , x) and γ (a F , X) are incomplete gamma functions, which can be solved numerically (Press et al., 1992, see Appendix A for derivation). For a F <0 (ǫ<ζ + η + 1, i.e., rather "flat" distributions at the upper model boundary) we can apply the recursion formula for the incomplete gamma function. The function shows that the deep flux depends not only on the (constant) sinking and remineralisation parameters and depth, but, in addition, on the exponent of the size distribution of particles at the surface (ǫ 0 ).

Discussion
Besides the obvious finding that the parameterisation of vertically increasing sinking speed can have a strong effect on the vertical distribution of biogeochemical tracers, the results presented so far suggest, that (1) the particle size distribution in the ocean interior might be a unimodal function of diameter, even if the upper boundary conditions were characterised by a power law and (2) that the vertical distribution of tracers and fluxes might depend on the surface size distribution.
3.1 The size distribution of particles: power law or unimodal?
Surface particle distributions are often described by power laws, inferred from straight lines in log-log plots of observational data sets (e.g., Jackson et al., 1997;Gin et al., 1999;Cavender-Bares et al., 2001;Gilabert, 2001;Quinones et al., 2003;San Martin et al., 2006). Sheldon et al. (1972) observed that the distribution of plankton particles especially in the deep ocean could be represented by a power law suggesting equal biomass in logarithmically increasing size classes (a so-called "flat" distribution). On the other hand, some theoretical and empirical evidence points towards unimodal (or sums of unimodal) size distributions (Lambert et al., 1981;Jonasz and Fournier, 1996). Note that roughly linear particle number spectra (on a log-log plot) can be quite deceptive, because even small deviations from a power law in the particle size spectrum may imply unimodal mass spectra (Jackson et al., 1997). The simple model of 198 size classes suggests that in the absence of any other size-dependent process beside sinking, intermediate depths will be characterised by unimodal particle mass distributions. The deeper the water, the bigger the dominant particle size. However, the model results presented here rely on the assumption, that individual particles do not change their properties (e.g., volume or density) during remineralisation and/or sinking. A different approach was carried by Zuur and Nyffeler (1992), who assumed that particle radius changes during remineralisation. Starting from a bimodal (volume) distribution as upper boundary condition, the dominant mode shifted from the smaller size towards larger size when integrating to 2000 m; nevertheless, it did not approach a power law.
If observed particle size distributions in the ocean interior can indeed be represented by a power law, we have to think of processes that especially remove the peaks of the theoretical size spectrum, that results from sinking and remineralisation alone, either by removing the most abundant particles, or by changing the particle properties (e.g., density, diameter) on their way downwards.

Processes affecting the size distribution of POM
To consider the effect of the different processes on particle properties and the size spectrum, it might be necessary to distinguish between different particle types, namely aggregates (formed via collision) and biogenic particles such as individual phytoplankton cells, fecal pellets, and other detrital material. Stemmann et al. (2004a,b) have extensively reviewed and investigated the different processes that might lead to changes in the particle size spectra and properties, especially with respect to aggregates. To name but a few, settling, coagulation and physical fragmentation might lead to changes in the particle size spectra. While coagulation and fragmentation might play a smaller role in the mesopelagic realm, where shear rates are low, zooplankton feeding and microbial degradation might play a larger role (see Stemmann et al., 2004a,b).
Zooplankton feeding can affect the size spectra in two directions: breakup during feeding will reduce the average particle size, while the ingestion of small cells, and the production of large fecal pellets would shift the main mode towards the upper end of the size spectrum. Both processes will depend on the size structure of the zooplankton community, the animals' preference for certain food size, and the animals feeding mode. For example, flux feeders will focus on the large, fast settling particles, whereas the filter feeders ingestion will depend on the size and geometry of their feeding apparatus. Stemmann et al. (2004a,b) found that especially flux feeding (and the associated production of fecal pellets) could be as important for variations in the mass of large particles as settling. If zooplankton grazing targets for the most abundant food and reworks the particles into smaller ones, this process may indeed remove the peaks in the particle size spectra. Aggregation of particles, on the other hand, would mostly promote unimodal distributions, by reducing the number of submicron particles (Lawler et al., 1980). Microbial degradation of biogenic particles can be an important process for the particle size distribution in the deep ocean, depending on the way the particles are remineralised. If the decay rate is constant with size, and particles are assumed to shrink during degradation, it will change the size spectra towards smaller particles (see Stemmann et al., 2004a,b). Zuur and Nyffeler (1992) in their model also assumed that particle size decreases during remineralisation, and a change in mass:volume relationships of decaying phytoplankton (indicating a decrease of the mass of individual particles) has been observed by Verity et al. (2000). Biogeosciences, 5, 55-72, 2008 www.biogeosciences.net/5/55/2008/ The effect of remineralisation on size may depend, however, on the time scales considered, and on the definition of "sinking particle". Because our model formulation is linked strongly to the way we consider POM and its decay, it seems worthwhile to take a closer look at the succession in the "detritosphere", as described by Biddanda (1988): Biddanda (1988) showed that particle decay happens in different stages: first, bacteria start to grow in the vicinity of the particle (day 0-4). Afterwards, the bacteria will colonize the particle (≈day 1-6) and convert its POC first to DOC by means of exoenzymes. Only a small fraction of the organic carbon will be incorporated into bacterial biomass, the rest will be respired (Biddanda, 1988). (Assuming fixed C: N ratios for bacteria, we might assume that an equivalent amount of organic nitrogen available to bacteria will then remain/be released as DON and/or ammonia.) The production of sticky extracellular mucopolysaccharides leads to aggregation of detritus particles and/or bacteria. At the same time there is an increase in protozoa which feed on the bacteria. The combined effects of microbial and protozoan activity finally lead to the disintegration of the detrital aggregates, which have largely disappeared by day 16.
The mass contained in detrital organic matter will thus first be converted to organic and inorganic forms, with some amount of this mass being in bacteria that are attached to the POC. I.e., in terms of mass there will be a shrinkage of the particle itself; the mass decrease will be less if we consider bacteria as a part of this particle. Considering the entire detritosphere, and neglecting diffusion out of the detritosphere, there will be a transformation of mass, but not a loss. Considering the diameter of the particle (more exactly: the detrital particle plus bacteria), the associated aggregation will lead to an increase in particle size. In the end (≈day 16), the disaggregation and consumption by protozoa will lead to the decrease in both numbers and mass of the detritus. Thus, on a roughly biweekly timescale there will be shifts in the size spectrum due to "leakage" of organic substances out of the detritus particle and subsequent remineralization, colonization by bacteria (if we consider these to be part of the particle), aggregation and disaggregation. In the end, a particle with mass m will have vanished entirely; this process, under the assumption of size-independent decay rate (Ploug and Grossart, 2000), and over a long enough time scale, will leave no imprint on the size spectrum.
Finally, the presence of calcifiers or diatoms (in contrast to phytoplankton without shells) will not only affect the density of the particles (and thus increase their sinking rate), but might also protect the organic tissue of the cells from degradation, probably leading to a decrease of remineralisation rate as the particles age. Both an increase in particle density as well as a decrease in its remineralisation rate will have an effect on flux profiles, and increase the flux of organic matter to the deep ocean. For more details on this process we would rather direct the reader to Armstrong et al. (2002) or Klaas and Archer (2002).
3.3 Imposing a power law everywhere in the water column As discussed above, both theoretical as well as observational evidence suggest power law as well as unimodal particle spectra. If particle size distributions in the ocean interior are indeed unimodal, approaches such as the power law size spectral approach by Kriest and Evans (2000) and Kriest (2002) are not fully appropriate to represent the evolution of particle size spectrum with depth, unless other processes (e.g., grazing, remineralisation) remove peaks in the size spectra. Kriest and Evans (2000) and Kriest (2002) simulated marine aggregates formed by coagulation, and assumed that aggregates at any depth were distributed according to a power law (Eq. A1). The model parameterised size-dependent sinking speed up to a size of 1 cm, and constant sinking for particles larger than this size. To investigate the effect of this power law assumption everywhere in the water column, we have calculated the "sedimentation aspect" of the model by Kriest and Evans (2000) and Kriest (2002).
Similar to Kriest and Evans (2000) and Kriest (2002), the model (hereafter referred to as K02) assumes an infinite power law size distribution of POM throughout the vertical model domain, size-dependent sinking up to a certain size, and constant sinking afterwards. POM is calculated in terms of mass (M, [mmol N m −3 ]) and particle numbers (N, [cm −3 ]), which change according to where and are functions of z and t, as defined in Kriest (2002). In particular, they depend on the size distribution parameters A and ǫ (Eq. A1), which are evaluated at any time and location from M and N, as described in Kriest and Evans (1999).
We have calculated this model numerically over time and depth, with a vertical grid of z=10 m, an upstream scheme for sedimentation, and a (forward) time step of t=0.25 h, for 1080 days. Simulated POM and size distribution are constant by this time. The upper boundary conditions for POM mass (M 0 ) and size distribution exponent (ǫ 0 ) are the same as in the analytic approaches presented above (see Table 2). Due to the infinite upper boundary of the size distribution in K02, sinking speed w 0 and input F 0 at the model boundary are now higher than those of the analytic approaches presented above, which were applied to a finite size range for POM. Parameters that describe the size dependence of particle mass (C and ζ ) and sinking speed (B and η) are the same as in the model approaches presented above (see Table 2), and also correspond to scenario "dSAM" presented in Kriest (2002).  Table 2 for descriptions of the "SAM" models.
According to K02's prerequisites, POM sinking speed will only increase up to a certain depth. Imposing a power law distribution over an infinite size range, below this depth the size spectrum is nearly "flat", and most of the mass is located beyond the upper limit for size-dependent sinking. As a consequence, the average POM sinking speed is constant (Fig. 3). The initial increase of average POM sinking speed with depth in K02 is stronger than in the models presented above.
We note that a direct comparison of this model with the analytic approaches presented above is hampered by several methodological differences: first, all of the above analytic approaches assume a finite particle size spectrum, while mass in "dSAM" is distributed over an infinite size range. Second, Kriest (2002) focused on the representation of "marine snow", which has different particle scaling characteristics than the particles presented in this work.
The effect of the infinite upper boundary is examined by a model that makes the same assumptions about particle scaling, distribution and sinking as scenario "dSAM", but assumes that POM is distributed only over a finite size range (scenario "dSAM-finite", see Table 2). The finite model's increase in sinking speed with depth is more moderate (Fig. 3); further, even at 4000 m POM does not achieve the maximum possible sinking speed of 153 m d −1 , because even with negative spectral exponents (implying rising slopes on plots of log mass vs. log size), not all particles are of maximum size. Thus, omitting the "upper tail" of the size spectrum of K02 has the effect of decreasing the sinking speed of POM and its increase with depth, the consequence being a lower normalised sedimentation.
In a second experiment with K02 (scenario "pSAM-slow") we parameterise the more porous marine snow, whose density decreases strongly with aggregate size (ζ =1.62). As a result, the increase of particle sinking speed with size is not as strong as in "dSAM", as given by η=0.62 (see also Kriest, 2002, and Table 2). This has a strong effect on simulated sinking speed and sedimentation: the increase in POM sinking speed with depth is quite low, especially in the upper few hundred meters (Fig. 3). Normalised sedimentation decreases strongly in the upper water column, and is more than an order of magnitude lower at 4000 m than in scenario "dSAM".
Summarising, imposing power law size spectra everywhere in the water column in K02 (instead of the more flexible size distributions described above) leads to a strong increase of POM sinking speed and sedimentation with depth, especially in the upper few hundred meters. This is only partly explained by the infinite upper boundary of the size spectrum in K02.

Variation of parameters and comparison to observations
Our results suggest that for a size spectrum of POM and in the absence of size dependent processes other than sinking, the mean sinking speed will increase with depth, and the depth dependence of the sedimentation flux can be described by Eq. (14). The sinking speed, and, consequently, the sedimentation profile, will further depend on the surface size distribution. Model results suggest that this may vary regionally, depending on the trophodynamic state of the ecosystem  "dSAM" and "dSAM-finite" parameterise "dense" particles, while "pSAM" parameterises more porous marine snow (see also Kriest, 2002). Spectral slope (ǫ) in "dSAM" and "pSAM-slow" is calculated as in Kriest and Evans (1999;"KE1999"). "dSAM-finite" assumes a finite boundary for the POM size range, and requires a numerical solution for ǫ.
This result agrees with that of other studies: regionally variable parameters of algorithms that describe sedimentation profiles may be necessary in order to fit observed sedimentation (Lutz et al., 2002;Francois et al., 2002) or biogeochemical tracers (Usbeck, 1999). Berelson (2001) also postulated regional variability of the exponent of the Martin function from data sets of sediment traps, but Primeau (2006) later showed that a large part of this variability could also be attributed to statistical effects. Parameterisations with regionally varying remineralisation length scales improved the simulated tracer distribution presented by Howard et al. (2006). Boyd and Trull (2007) present a comprehensive overview over the possible mechanisms that may alter the regional flux pattern, and on the methods (and their limitations) applied determine the export profile.
In this subsection we investigate the three different models (constant sinking speed, Martin's sedimentation curve, and the analytic approach of the spectral model) with respect to their sensitivity to the parameters. We further compare the simulated flux ratios (sedimentation divided by sedimentation at the upper model boundary) with observations derived from Th-export, moored and floating sediment traps. The traps were deployed at least one year in the central Arabian Sea (AS-C, Lee et al., 1998), in the North Pacific (OSP, Wong et al., 1999), at the Bermuda Atlantic Time-Series station (BATS; data after Lutz et al., 2002;Conte et al., 2001) and at the Hawaii Ocean Time-Series station (HOT; data after Lutz et al., 2002). We have further added three sedimentation profiles collected during roughly biweekly in-tervals during the North Atlantic Bloom Experiment NABE (Martin et al., 1993, available from the US-JGOFS website, http://usjgofs.whoi.edu/), and flux ratios determined from carbon flux collected with neutrally buoyant sediment traps (Buesseler et al., 2007), which were deployed at two stations in the Pacific (ALOHA, K2).
Thus, the observations span a wide range of different regimes, from mainly oligotrophic (e.g., HOT, AS-C) to bloom regimes (e.g., NABE). For the model-data comparison we have always used the flux of particulate organic carbon; we divided all observed fluxes by the shallowest observed flux (usually at 100 to 150 m depth).

Constant sinking speed
As noted in the introduction, a number of biogeochemical models have employed constant, albeit model-specific, POM sinking speeds. The choice of the particular constant sinking speed is often explained by observations of individual particles, or derived from observations of sediment traps, (e.g., located at 150, 200 and 300 m at the BATS size in the North Atlantic; e.g., Doney et al., 1996). To examine the model's sensitivity to variations in the constant sinking speed of POM, we set its sinking speed such that it matches that of the model with 198 size-classes in 300 m and 1000 m, resulting in an average sinking speed of 22.3 and 45.27 m d −1 . A change in average sinking speed simply shifts the region of mismatch with respect to the spectral model's solution or observed particle fluxes (Fig. 4). Thus, deriving a constant sinking speed from sediment trap observations at a certain depth (similarly, for location) will bias a model towards this depth, but probably be of little predictive power for domains far below or above.

Martin's curve and spectral model
To test the models' sensitivity to changes in surface biology, we changed the spectral exponent of the surface boundary condition, ǫ 0 by ±0.5. This change in surface size structure corresponds to changes in w 0 ; to see how changes in ǫ 0 convert to changes in w 0 , divide Eq. (A15) in Appendix A by Eq. (A10). It also affects the exponent of the Martin curve, r/a, because a=w 0 /z 0 , and thus r/a=z 0 r/w 0 =f (ǫ 0 ). In particular, the increase (decrease) of ǫ 0 by ±0.5 converts to exponents of 1.598 (0.358).
In both models (continuous size spectrum and Martin's curve) the steepening of the spectrum (ǫ 0 =4.96) results in a stronger attenuation of the normalised sedimentation with depth. More organic matter reaches the deep ocean when the size spectrum at the ocean surface becomes flatter (ǫ=3.96). The effect is much more pronounced in the Martin model (Fig. 5). In this model, the flux ratio at 4000 m decreases by more than an order of magnitude, when the exponent is increased by 0.5. The range of flux ratio in both models encompasses the observed ratios ( Fig. 5-6); the size spectral model additionally shows a quite good fit to observed flux ratios at BATS (Fig. 6). Summarising, a model with constant sinking speed of POM may be biased towards observations and/or the biogeochemical settings at a specific location or depth, and does not reflect observed flux ratios at all depths simultaneously. The models that simulate an increasing sinking speed with depth much better reflect the observed flux ratio. Especially the size spectral model is quite close to observations in the upper 400 m at the BATS site. Because it further shows a much lower sensitivity to variations in the surface size structure than the model with linearly increasing sinking speed, its variability is more similar to that of the (few) observations presented here.

The conjunction between remineralisation and sinking parameters
So far, we have attributed any changes in the flux profile to variations in the surface size structure, as given by changes in Biogeosciences, 5, 55-72, 2008 www.biogeosciences.net/5/55/2008/ ǫ 0 or w 0 . Depending on the model, it is, however, difficult to untangle the effect of remineralisation rate and sinking speed. With respect to the sedimentation and remineralisation profiles described by the Martin function a decrease in remineralisation rate would be equivalent to an increase in surface POM sinking speed (w 0 ) via the relation a=w 0 /z 0 (see above). In this paper, we interpret the imposed change in the exponent as a change in the surface boundary conditions (ǫ 0 or w 0 ); it could however, also be interpreted as a change in remineralisation rate r under fixed ǫ 0 , in particular r=0.056 for r/a=1.598 and r/a=0.012 for b=0.358. Thus, the experiment with different exponents for the Martin curve (Fig. 5) also describes the sensitivity of the Martin function to a change in remineralisation rate.
Things are slightly different for the size spectral model, as here r is not directly related to the initial POM sinking speed, or the upper ocean size distribution (see Eq. A17, and definitions for a F , x and X given in the appendix). Roughly speaking, under a fixed ǫ 0 =4.46 (corresponding to w 0 =3.52) an increase (decrease) in r to the values mentioned for the power law (r=0.056 and r=0.012, respectively) will decrease (increase) the flux ratio especially at greater depths (Fig. 7), whereas a change in ǫ 0 will change the flux ratio especially at shallower depths.
Finally, the size spectral model will also depend on other parameters, such as the size range considered. Theoretically, we think it is possible to disentangle the effects e.g., of ǫ 0 and the POM scaling parameters (e.g., ζ or η of the mass vs. diameter or sinking speed vs. diameter relationship). The success of such an approach will, however, strongly depend on the data set used to constrain the parameters. This will be presented elsewhere, together with a more detailed modeldata comparison.

Alternative sets of observations
It is, however, difficult to decide about the appropriate flux parameterisation from direct comparison with observations, as measurements of sedimentation are sparse and may be subject to errors. Sediment traps may miss up to 80% of local in situ flux Scholten et al., 2001). Possible cause may be the loss of POM to the dissolved phase in the collecting cup (Kähler and Bauerfeind, 2001), or hydrodynamic effects associated with sediment trap design (Gust et al., 1994). Buesseler et al. (2007) presented results from neutrally buoyant sediment traps, which are supposed to overcome some of these problems. They observed different flux attenuation profiles at two different stations in the North Pacific, Fig. 6. As Fig. 4, but for the model that applies a continuous size spectrum, together with analytic integration over depth (Eq. 14). Black line: standard scenario (ǫ 0 =4.46). Red lines: experiments with ǫ 0 =4.96 and ǫ 0 =3.96.
with station ALOHA (near Hawaii) being characterised by strong vertical flux attenuation, while the subarctic gyre (station K2) was characterised by a high transfer efficiency. They attributed the differences between the two site to trophodynamic and/or ballasting effects. Our results so far suggest that differences in the surface size structure can explain the differences in flux attenuation (Fig. 6, right panels). Especially the size-spectral model is in quite good agreement with the results obtained by Buesseler et al. (2007).
A different approach to assess model performance can be found in the comparison with nutrient profiles. Inorganic nutrients are (relatively) easy to measure, and the data sets of nutrients (and oxygen) are already rather dense. Because the flux divergence with depth and the nutrient profile are tightly intertwined, a possible solution could lie in the application of global coupled (physico-biogeochemical) models with different settling characteristics, integrated over long time scales. Different assumptions about particle settling characteristics will then translate into different regional and global nutrient profiles. These can be compared to global data sets of nutrient observations and, given a reliable circulation field, help to asses the flux algorithms and their parameters (as, e.g., in Kwon and Primeau, 2006). This will be investigated in a future study.

A climatological year and annual averages for POM and sedimentation
So far, we have only investigated systems in equilibrium, i.e., systems that fulfil the condition ∂M/∂t=0. More important for global, long-term simulations with seasonally varying forcing is the non-equilibrium case, i.e., a time-varying POM or flux concentration and ∂M/∂t =0. We investigate this by means of an upper boundary condition for POM and w 0 derived from the output of a 1-D-model simulated for a site in the northern North Atlantic (Kriest and Oschlies, 2007). A distinct ensemble of particles produced in the surface layer on a certain day (i.e., POM with a certain distribution or sinking speed; here named POM t , where t denotes the time) will travel along its characteristic trajectory downwards. We can expect a certain amount of this POM to arrive at depth z by time t + t . The time it takes for POM t to arrive at this depth, t, is determined by its sinking speed; the amount of POM t that arrives at this depth is given by its remineralisation rate. Considering long enough time scales and no horizontal processes, sooner or later every POM ensemble (more precisely: a fraction of it) will arrive at a depth z -even if it travels very slowly, and has been created years before. Considering climatological years, we can thus average the POM Biogeosciences, 5, [55][56][57][58][59][60][61][62][63][64][65][66][67][68][69][70][71][72]2008 www.biogeosciences.net/5/55/2008/ sedimentation over one year, and get the average annual flux and POM, even under time varying POM concentration and sinking speed at the surface.
To illustrate this we have taken a POM mass M 0 and spectral exponent ǫ 0 at the surface from a model simulation of size dependent phytoplankton physiology (Kriest and Oschlies, 2007). We have scaled ǫ 0 of Kriest and Oschlies (2007) to match the parameters applied in this work. From these surface boundary conditions (see Fig. 8, panel A), we have calculated (1) 198 size classes with numerical integration over depth; the model was run with a time step of ≈10 min, z=10 m for 101 years, of which we present the last year, (2) 198 classes with analytic vertical integration (Eq. 10), and (3) a continuous size spectrum with analytic vertical integration (Eq. A12). For the latter two approaches we distribute surface POM immediately over depth (similar to the approach by e.g., Maier-Reimer, 1993), whereas in reality (and in the model with numerical integration over depth) it would take some time for surface POM to reach a certain depth. This is evident from comparing the panel B of Fig. 8 with panels C and D. However, the annual averages of vertical distribution of POM, sedimentation and sinking speed are almost the same for the three approaches (Fig. 8, lower panels E-G).
While the approach to immediately distribute export flux in the vertical cannot be used to resolve the short-term (daily to seasonal) changes in sedimentation, or might not be valid in areas of high eddy activity and/or in areas of high vertical mixing, it is nevertheless very efficient in models that resolve rather coarse time and space scales. Summarising, if we disregard the temporal resolution of deep POM distribution and sedimentation, we can simulate the flux at any depth without having to evaluate POM at any depth. Applying the sizecontinuous approach (Eq. 11), even the evaluation of (many) distinct size classes would not be necessary. This can be of advantage e.g., in global models that are driven by a climatological forcing and simulated over long time-scales.

Conclusions
We examined two different descriptions of POM fluxes -namely the Martin curve and constant particle sinking speed -which have been often applied especially in threedimensional large-scale models, for their intrinsic assumptions and ability to fit observed profiles of normalised sedimentation. A third parameterisation was developed from the assumption of a power law particle size spectrum at the surface ocean, size dependent sinking and constant remineralisation, and also compared to observed sedimentation.
One result of our study is that a constant sinking speed of POM combined with a constant remineralisation rate cannot reproduce observed fluxes at a variety of locations (depths) simultaneously. Our theoretical analysis of mechanistic principles and the model-data comparison both suggest, that av- Fig. 7. Flux ratio (sedimentation divided by upper model boundary condition) for the model that applies a continuous size spectrum, together with analytic integration over depth (Eq. 14). Black line: standard scenario (ǫ 0 =4.46, r=0.03d −1 ). Red lines: experiments with changes in surface ocean size distribution, ǫ 0 , assuming constant remineralisation rate, r=0.03d −1 (see also Fig. 6 for comparison to observations). Blue lines: experiments with changes in remineralisation rate r to r=0.012 and r=0.056, respectively, assuming constant ǫ 0 =4.46. erage POM sinking speed increases with depth (analogously: remineralisation rate decreases with depth). This is accounted for by both the Martin curve or the flux curve derived from a power law surface size distribution of POM. From the comparison of these two parameterisations with observations so far we cannot decide for any of the two functions. This will be subject of a further analysis, which makes use of different particle flux models coupled to a more detailed three-dimensional physical model, and compares the results to observations. However, the size resolving model provides a mechanistic explanation for the increase of POM sinking speed with depth.
We further find that the regional variation of remineralisation length scales suggested by other authors (Usbeck, 1999;Lutz et al., 2002;Howard et al., 2006;Buesseler et al., 2007, e.g.) could be explained by variations in the surface particle size spectra. While the size resolving model provides a consistent and mechanistic link between upper ocean biogeochemistry and the deep ocean, to our knowledge such a link has not been established for the Martin curve (although via the dependence of its exponent on surface POM sinking speed w 0 a model with regional or temporal variation in w 0 could result in regionally varying remineralisation length scales). Summarising, we suggest one of the two functions (Martin or size resolving model) for the coupling between surface and deep ocean, preferably with a regional variation of remineralisation length scales. The size-resolving model presented here provides such a variation in remineralisation length scales as a mechanistic interpretation of upper ocean size effects on sedimentation profiles.

An analytic evaluation over depth and size
In principle, we follow the same approach as for the sizediscrete model, but on a continuous particle length scale: assume that we can describe the number of particles per unit length d by Biogeosciences, 5, 55-72, 2008 www.biogeosciences.net/5/55/2008/ A and ǫ are supposed to vary with depth and time. Assume that C d ζ describes the individual mass, and that the coefficients of this function are constant with depth and time. We can then represent the distribution of particle mass m per unit length Assume there is a linear decay of individual particle mass with time, and that the decay rate does not depend on diameter, time or depth: Now assume that particle sinking characteristic of a particle of size d depends on diameter: w=B d η , and that the parameters (B>0, η>0) do not change with time or space. For the time rate of change for particle mass we then get ∂m ∂t +B d η ∂m ∂z +r m=0 (A4) Note that -as for the previous models -this formulation implicitly assumes that particle (number) loss rate due to remineralisation is the same as that for mass: i.e., the particles do not get less dense (or less filled with organic matter), but all the losses of mass are concentrated in a few selected particles, that disintegrate immediately.
In equilibrium (∂m/∂t=0) and with constant boundary condition m 0 =m(t, 0)=const, we get m(z ′ )=m 0 e − r z ′ B d η (A5) Assuming a size spectrum at the upper model boundary, we can represent m 0 by Eq. (A2) as: In this case The mass of the total particles ensemble (from size d 1 to size d L ) is given by the integral over the size range: To integrate this function, we substitute the exponent of e by: With The integral term in Eq. (A11) is the difference of two incomplete gamma functions γ (a, X)−γ (a, x) (e.g., Press et al., 1992), for which we can solve numerically, provided ǫ 0 >ζ +1: Our experiments carried with constant boundary conditions, and with parameters as described in the previous sections show, that the numerical solution of the two terms in brackets takes, on average, about ≈2−3 (first term) and ≈7−8 (second term) iterations, with a maximum of 15 iterations.
In analogy, and with the same substitution, we can evaluate the flux at any depth z. Multiplying both sides of Eq. (A4) with particle sinking speed w we get the analogous solution for the flux per unit size at any depth, f (z) as function of surface flux f 0 and sinking and decay coefficients for that size: The integral Eq. (A13) over the whole size range of particles is then With the flux at the upper model boundary given by x, X as defined above and we then get (provided ǫ 0 >ζ + η + 1), In case a F <0 (ǫ 0 <ζ + η + 1), we can apply the recursion formula for the incomplete gamma function, i.e.: γ (a, x)= γ (a + 1, x) + e −x x a a (A18) Figure 1 shows that the results (POM and sedimentation) of the analytic solution agree very well with the results of the numerical, size resolved model. Thus, given a steadystate power law distribution at the base of the euphotic zone and neglecting any impact of advection and mixing on sedimentation and a relation between particle sinking speed and diameter, we can evaluate the POM and its sedimentation at any depth. Finally, dividing Eq. (A17) by Eq. (A12) gives the average POM sinking speed: