On the variability of the exponent in the power law depth dependence of POC flux estimated from sediment traps

Using the power law relationship J ð z Þ ¼ J 0 ð z = z 0 Þ (cid:1) b relating the particulate organic carbon (POC) ﬂux in the water column at depth z to the export ﬂux J 0 at z 0 ¼ 100m, I reanalyzed the sediment-trap data set compiled by Berelson [2001. The ﬂux of particulate organic carbon into the ocean interior: a comparison of four U.S. JGOFS regional studies. Oceanography 14, 59–67]. The goal is to better understand the variability of the exponent b . I show that the usual approach of estimating parameters for each station separately and then pooling the estimates from all stations confounds within-station estimation errors with true regional variability and so tends to inﬂate the apparent variance of b . Furthermore, I show that the correlation between b and J 0 observed by Berelson [2001. The ﬂux of particulate organic carbon into the ocean interior: a comparison of four U.S. JGOFS regional studies. Oceanography 14, 59–67] is spurious and attributable to the fact that the estimation errors for b and J 0 are correlated. A two level mixed-effects regression model is introduced to properly take into account the contribution to the variability of the POC ﬂuxes of within-and between-station effects. The analysis shows that the variability of b across stations is only a small contributor to the POC ﬂux variance in the sediment trap data set. In fact the additional POC ﬂux variance captured by a model with regionally varying b in comparison to a simple model with a ﬁxed b is not statistically signiﬁcant. The best estimate for the ﬁxed b applicable to all stations is 0 : 70 (cid:2) 0 : 08 (95% C.I. 90 d.f.).


Introduction
The depth dependence of the downward flux of particulate organic carbon (POC) in the ocean, JðzÞ, is often parameterized by a power law where z 0 is some reference depth, typically near 100 m, at which the export flux J 0 Jðz 0 Þ can be explicitly modeled or measured.The functional form in Eq. ( 1) was initially proposed by Martin et al. (1987) as a way of summarizing sediment-trap data, but has since been adopted as a useful parameterization of carbon export by sinking particles in many global biogeochemistry and ecosystem models (e.g.Sarmiento et al., 1993;Yamanaka and Tajika, 1996;Najjar and Orr, 1998).The parameter b controls the depth of particle remineralization.For large values of b there is a rapid decrease of the particle flux with depth, and most of the remineralization happens in the top of the water column just below the euphotic zone.
For smaller values of b the particle flux decreases more slowly, and a larger fraction of the flux reaches the sea floor.In models, the efficiency of the biological pump and the resulting partition of carbon between the atmosphere and ocean is sensitive to b (Kwon and Primeau, 2006).Constraining b from data is therefore important.
In a summary of JGOFS sediment-trap data using Eq.(1), Berelson (2001) found variability in the estimates of b from different stations and suggested that this variability might point to regional differences in the biogeochemical controls of the attenuation of the POC fluxes with depth.Furthermore, he noted that a significant fraction of the variability of b b (hats are used to indicate estimates of the ''true'' value) can be explained by variations in b J 0 .Taken at face value, this result implies that regions of higher export have a higher fraction of the flux remineralizing in the upper water column.A similar observation was made by Lutz et al. (2002), who suggested that the flux of POC to depth was more efficient during periods of low primary production after noting a correlation between the export flux J 0 and the ratio, JðzÞ=J 0 (what they call the s ratio ).Berelson suggests several possible mechanisms that might explain the correlation: (1) varying rate constants for POC degradation associated with differences in molecular compositions; (2) varying kinetics due to different oxygen, temperature and mineral distributions in the water column; (3) differences in microbial and heterotrophic communities; and (4) variations in settling velocities, all of which could plausibly correlate with the spatially varying biological production in the top 100 m.
Here I show that estimating pairs of b and J 0 values separately for each station and then pooling the resulting estimate pairs as was done by Berelson (2001) and previously by Martin et al. (1987) tends to inflate the apparent station-tostation variability of b because it confounds errors in the estimates of b with true regional differences.I show also that such an analysis produces estimates of b and J 0 with correlated errors, making it difficult to interpret the apparent covariance between J 0 and b noted by Berelson (2001).One of the goals of this article is therefore to introduce an improved estimation method that can properly distinguish within-station errors from across-station variability in b and J 0 .I introduce a two-level mixed-effects regression model (Raudenbush and Bryk, 2002) to properly take into account the combined contribution of within-and between-station effects to the variability of the POC flux measurements.Withinstation effects include measurement errors and other effects that contribute to deviations from the power-law flux profile.Between-station effects include differences in J 0 and b from station to station.A statistical analysis based on the new model shows that with some caveats concerning the North Atlantic Bloom Experiment (NABE) station, the JGOFS sediment-trap data set compiled by Berelson (2001) is consistent with the hypothesis that there is no detectable covariance between b and J 0 .Moreover, the analysis shows that a model in which b varies from station to station does not account for a significantly larger fraction of the POC flux variance than a simpler model with a fixed b.

Expected variance and covariance between the estimates of the power-law exponent and POC export flux
Although Berelson (2001) fitted the power law directly to the flux measurements, it is useful first to transform the measured fluxes to a logarithmic scale so that the resulting model becomes linear.Following Martin et al. (1987) and Armstrong (2002), we take the logarithm of (1) to obtain log JðzÞ ¼ log J 0 À b log z=z 0 . (2) Applying Eq. (2) to the discrete set of sediment-trap data from one particular site yields a linear estimation problem of the form where where z i are the sediment-trap depths, and e i are the errors assumed to be independent and identically distributed with variance s 2 .The least squares solution of Eq. ( 3) is b and the expected variance-covariance matrix of the parameter estimates is s 2 ðX T X Þ À1 (e.g.Faraway, 2005).Explicitly, the variance-covariance matrix for the parameter estimates is where When parameter estimates obtained from different locations are combined, (e.g.Berelson, 2001;Martin et al., 1987) it is important to take into account the within J 0j obtained from different station is the sum of the across-station variance-covariance plus the within-station variance-covariance of the estimates.Only the across-station variance-covariance represents true geographical differences in b and J 0 .The within-station variance-covariance, given by the average across all stations of the matrix expression in (6), represents a bias due to within-station estimation errors.
From the off-diagonal terms in (6) we see that the sign of the bias in the covariance depends only on the depths of sediment traps in relation to the reference depth z 0 .Because the fluxes are referenced to the export flux at 100 m near the top of the water column, the depths z i are generally much greater than z 0 , and as a result, the expected covariance is positive for all sites.Importantly, the bias does not decrease as the number of stations in the ensemble increases.It decreases only as the number of flux measurements within each station increases.Because the number of flux measurements at any given site in the Berelson (2001) data set is small (n ranges from 3 to 12), the bias is potentially significant.
In Fig. 1, I show plots of the parameter estimates for the data set compiled by Berelson (2001) separated by region, but estimated separately for J 0 .From Eq. ( 6), the slope of the tilt is given by which is fairly uniform for all stations even though the trap deployment depths vary somewhat between stations.Because the uncertainties in the parameter estimates are positively correlated at each site, we expect the combined parameter estimates from all four regions (17 stations) to be positively correlated even if the true parameters are independent.Berelson (2001) computed a regression line for the scatter plot of his estimates of b against J 0 and found what appeared to be a statistically significant positive slope.However, Eq. ( 6) shows that the errors in the estimates are not independent and the usual significance tests are not applicable.The slope of a regression line through a scatter plot is given by (e.g.Janke and Tinsley, 2005) where d

Cov and d
Var are used to denote the sample covariance and variance of the 17 pairs of parameter estimates.If the true b and J 0 are independent, the numerator in (9) will be given by the average over the 17 stations of the off-diagonal term in (6).As already noted, this term is positive for all stations, so the average will also be positive.That a positive slope in the scatter plot of b vs. J 0 is clearly visible in Berelson's (2001) analysis is expected from the correlation of the errors.This feature of the covariance of b and log J 0 demonstrates the need for a more careful analysis to ascertain whether b is spatially varying and whether it co-varies with J 0 .
It is easy to understand the origin of the biased covariance between b b and b J 0 .A random error that leads to an over prediction in the export flux at the surface leads to a positive error in the b-estimate because the POC flux profile has to decrease more rapidly in order to be close to the deeper flux measurements.Similarly, a random error that leads to an under-prediction in export production at the surface leads to a negative error in the b-estimate because in this case the flux profile has to decrease more slowly if it is to be close to the deep flux measurements.Mathematically this effect is captured by the non-orthogonality of the columns in a regression design matrix X in Eq. ( 4).By choosing a deeper reference depth, ðz 0 % 670 mÞ it is possible to make the columns of X nearly orthogonal, but the model is then conceptually less useful because it relates the sediment flux profile to a depth that is far from where the particles originate.More fundamentally, the approach of first estimating parameters separately for each station and then combining the resulting estimates to study their variability confounds within-station errors with the ''true'' regional variability of the parameters.Unless the withinstation error variance is much smaller than the sought after variance-covariance of the ''true'' parameters, the usual estimation approach is inadequate for studying the spatial variability of POC flux profiles.
In the next section I introduce a two-level mixedeffects model (Raudenbush and Bryk, 2002) to properly take into account the contribution of within-station and between-station effects to the variability of the POC flux measurements.Using the two-level model, I reanalyze the sediment-trap data set used by Berelson (2001) to: (i) test if there is a statistically significant correlation between b and J 0 , and (ii) test how much of the often noted station-tostation variability in the estimates of b is due to true geographical differences in b.

Two-level mixed effects model
To properly treat the uncertainty associated with modeling POC fluxes measured from sediment traps using Eq.(1), we introduce a two-level mixed-effects regression model.A comprehensive introduction to this type of model can be found in Raudenbush and Bryk (2002).The two-level mixed-effects model I use takes into account uncertainty in the predicted fluxes due to the across-station variability of b and J 0 in addition to uncertainty due to the usual within-station effects such as measurement errors.In the usual parameter estimation problem based on ordinary least squares (OLS) the coefficients b and J 0 are assumed fixed, whereas here we treat them as random variables.Mixed-effects models have become well established in several research fields and routines for estimating parameters are standard in many statistical software packages.Mixed-effects models are also known as random-effects models in the biometric literature, as random-coefficients regression models in the econometrics literature and as covariance components models in the statistical literature (Raudenbush and Bryk, 2002).
To keep the notation manageable I denote the ith POC flux measurement at the jth station by J ij and define The first level of the model (level-1) describes the within-station flux variability where g 00 and g 10 are, respectively, the acrossstation mean log J 0 and b values and where u 0j and u 1j are random effects that vary from station to Note that estimating the parameters by ordinary least squares is not appropriate, because the error terms, u 0j þ u 1j Z ij þ r ij , do not have equal variance either within or across stations and moreover are not independent within stations.A proper estimation of g 10 (i.e. the mean b value) in Eq. ( 14) takes into account the different precisions with which b is determined at each station.It thus produces a better estimate than simply averaging the b values obtained from different stations (e.g.Berelson, 2001) or as in the case of Martin et al. (1987) by pooling all the measurements and estimating only a fixed J 0 and b.I estimated the parameters in (14) using a maximum likelihood approach as implemented in the subroutine lme (Bates and Pinheiro, 1998) in the statistical software package R.The parameter estimates are listed in the row labeled full in Table 1.The numbers in parenthesis are the 95% confidence intervals for the estimated parameters.The estimate for the fixed b value is b g 10 ¼ 0:72 AE 0:1 (95% C.I.), and the estimated standard deviation of b across stations is sd½u 1j ¼ 0:14, with lower and upper 95% confidence limits of 0.06 and 0.31.In comparison, The full and null model parameters were estimated from the full data set.The full Ã and null Ã model parameters were estimated from the data set in which the two shallow POC fluxes at the NABE station (Buesseler et al., 1992) were excluded.In the null and null Ã models, b is fixed and is not allowed to vary across stations.J 0 is expressed in mmol C m À2 day À1 .The quantities in the parentheses are the lower and upper limits of a 95% confidence interval.The full data set contains 17 stations and a total of 108 flux measurements.

ARTICLE IN PRESS
the standard deviation of the b values estimated in Section 2 by ordinary least squares is 0.22, approximately 60% higher.This is consistent with the fact that the variance of the estimates obtained by ordinary least squares is equal to the sum of the estimation error variance plus the across-station variance of the true b.The variance of the OLS estimates of b is 0.047 and the across-station average of the estimation error variance for b from the OLS formula is 0.031.If we try to correct the OLS variance by subtracting the average of the estimated within-station error variance we obtain an across-station standard deviation of b sd½b ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0:047 À 0:031 p ¼ 0:13 which is closer to our estimate based on the mixed-effects model, but appears to be biased low.The plot of b b j b g 10 þ b u 1j vs. c J 0j expfb g 00 þ b u 0j g estimated from the two-level mixed-effects model is shown in Fig. 2.
Since the standard deviation of b must be positive, the confidence interval for the acrossstation standard deviation of b excludes zero.To test the statistical significance of including the across-station deviations in b, we cannot simply check to see if its confidence interval excludes zero.Instead we construct a null model in which acrossstation deviations of b are set identically to zero (i.e.u 1j 0) and compare it to a full model by a likelihood ratio test (e.g.Faraway, 2005).This test yielded a p-value of 0.15, indicating that there is a 0.15 probability that the additional POC-flux variance captured by the full model is due to chance alone.Furthermore, a closer look at the data shows that the already weak evidence for a variable b is strongly influenced by a single station.The point for the North Atlantic Bloom Experiment labeled NABE in Fig. 2 stands out as a clear outlier.At NABE, two shallow flux estimates, at 150 and 300 m, from Buesseler et al. (1992) based on the thorium method are a factor of 3.5 and 6.5 greater than those obtained by Martin et al. (1993) 1).A likelihood ratio test, comparing the null Ã and full Ã models yielded a pvalue of 0.7807.Without the Buesseler et al. (1992) shallow NABE fluxes based on the thorium method there is no evidence for rejecting the null model with a fixed b.
A plot of the b and J 0 values estimated from the data set excluding the two shallow thorium-based fluxes at NABE is shown in Fig. 3. Without the high shallow fluxes, the NABE station no longer stands out as an outlier.I hypothesize that including the Buesseler et al. (1992) thorium-derived flux data at NABE without having thorium-based estimates of the flux deeper in the water column biases the exponent towards a higher value because the power   On the grounds of simplicity, the best model appears to be the null model in which b is fixed.The best estimate for the fixed b is 0:70 AE 0:08 (95% C.I. 90 d.f.).The value is robust to the inclusion or exclusion of the level-2 random effects for b or to the exclusion of either the shallow NABE fluxes based on the thorium method or the NABE station data altogether; all cases produced fixed b estimates well within the quoted uncertainty (Table 1).

ARTICLE IN PRESS
From the full model I also estimate t 01 t 10 ¼ 0:67 with a 95% confidence interval ð0:03ot 01 o0:92Þ.Recall that t 01 is the true across station covariance between log J 0 and b.Since the confidence interval is positively bounded, there appears to be some evidence that stations with high 100-m export fluxes tend to have higher b exponents as observed by Berelson (2001).However, the correlation is entirely due to the influence of the two shallow fluxes based on the thorium method at the NABE station.This is evident from Figs. 2 and 3, where the apparent correlation disappears when the shallow thorium-based fluxes are excluded from the NABE station.Without the two shallow fluxes based on the thorium method at the NABE site (null Ã in Table 1), the estimate for the covariance is t 01 ¼ 0:13, with a ðÀ0:86; 0:92Þ 95% confidence interval, confirming that the correlation is not robust.The tendency for faster remineralization with depth in regions of high export flux is not systematic across the entire data set and is sensitive to the possible inconsistency of using different flux measurement methods at different depths.
I have highlighted the possible inconsistency of mixing POC fluxes based on sediment traps and the thorium method at the NABE site.One should not infer from this that the NABE site was exceptional in Berelson's (2001) data set for having thoriumbased shallow fluxes; in fact all stations in the data set included thorium-based shallow POC fluxes.I have no evidence for discounting the thorium-based shallow fluxes at the NABE site in favor of the sediment-trap measurements.The point I wish to make is that the apparent correlation between d log J 0 and b b is sensitive to inconsistencies in the shallow vs. deep flux measurements, and that the correlation noted by Berelson (2001) is not a systematic feature of the parameters once the errors in these parameters are properly represented.

Discussion
I have shown that the usual approach of fitting power law exponents to the vertical profile of POC flux measured from sediment traps overestimates the variability of the exponent because it confounds estimation errors with regional differences in the true b.I introduced a two-level mixed-effects regression model to properly model the uncertainty associated with predicting POC fluxes with a powerlaw profile.The mixed-effects model takes into account uncertainty associated with the acrossstation variability of b and J 0 as well as the usual within-station uncertainty associated with measurement errors.Contrary to previous studies I found that a model with a fixed exponent b applicable at all sites adequately fits the data.The additional variance in the POC flux captured by a model in which b varies across stations is not statistically significant.Our best estimate for the exponent is b ¼ 0:70 AE 0:08 (95% C.I. 90 d.f.).
There is of course little reason to expect that a universal b value would be applicable at all stations.In fact the existence of a universal b would be surprising given the complex transformations that affect particles as they rain down from the surface.Our result should not be interpreted as meaning that there is a universal b value.What it tells us is that the variability in the POC flux data set collected from sediment traps is dominated by effects other than the variability of b across stations.This is an important point to bear in mind if, as suggested by Berelson (2001), differences in estimated b values are to be used as guidance for seeking regional differences in the processes that might systematically affect the attenuation of POC fluxes with depth.To detect significant differences in the power law exponent a larger number of measurements per station is needed than was typical in the Berelson (2001) data set.Another possibility for decreasing the exponent's estimation variance would be to pool stations into subgroups and test for statistically significant differences between the groups.Finally, working directly with the time-series data as opposed to the annually averaged fluxes might reveal differences in b as a function of time.For the purpose of the analysis I presented here, any such temporal variability was lumped into the withinstation measurement error.
Our estimated b-value is lower than the average value (b ¼ 0:82 AE 0:16 1 s.d.) quoted in Berelson (2001).The reason for this is that our estimation procedure based on the two-level mixed-effects model explicitly takes into account the differences in precision with which b is determined at each station.Bishop (1989) reviews several parameterizations of POC flux, and our estimate of b falls in the range of those that are cast in terms of the power law profile ( b b ¼ 0:628-1.0).Without proper uncertainty estimates for the b values estimated in previous studies it is difficult to judge if the differences in b values are significant.For example, our estimate is lower than the Martin et al. (1987) open-ocean composite (OOC), b ¼ 0:858, but without a reliable uncertainty estimate for the OOC it is difficult to assess the significance of the difference.Based on the published results of Martin et al. (1987), I estimated the uncertainty of b based on two methods: one from the quoted r 2 value for the OOC and one from the standard deviation of the b values fitted for each station separately.Appendix A show how I backed out a formal standard error for the residual of the OOC using the quoted r 2 value in Martin et al. (1987).The result is s OOC ¼ 0:059.Which yields a 95% confidence interval of ð0:7395obo0:9765Þ.This overlaps with our 95% confidence interval ð0:62obo0:78Þ.However, the formal standard error based on the OLS formula is likely to be an underestimate because it is based on the incorrect assumption that the OOC residuals are independent.As our two-level mixed-effects model makes explicit, the residuals for the OOC fit depend on a variable J 0 value that varies across stations but not within stations.Martin et al. (1987) analysis for the OOC assumes a fixed J 0 value for the pooled data.This suggests that the difference is probably not significant.At the other extreme, if I use the standard deviation of the b values for the individual stations that went into the OOC, I obtain a standard error of 0.12, which would suggest that the difference between the OOC and our estimate is insignificant.However, this uncertainty is an overestimate since it does not properly take into account the larger number of degrees of freedom in the pooled data set.While inconclusive the above considerations make it clear that proper uncertainty estimates are needed in order to detect real regional differences in the power-law exponent.Multi-level mixed effects regression models are ideally suited for this task.
The 95% confidence interval for our estimate, ð0:62obo0:78Þ, is relatively large suggesting that carbon cycle models have some latitude in tuning the b value in order to improve the model performance.Global models might be able to provide additional constraints on the power law exponent by optimizing the exponent so as to best reproduce global dissolved inorganic carbon, phosphate or alkalinity fields (e.g.Usbeck et al., 2003).Ultimately a better understanding of the processes that affect POC fluxes in the water column should produce a mechanistically based model to replace the empirical power law parameterization (e.g. the mineral-ballast model of Armstrong et al., 2002).I believe that mixed-effects regression models can also be usefully applied to test the statistical significance of parameters in the mineral-ballast model since it can easily be extended to take into account withinand across-station differences in mineral fluxes.deviation of the residual can be backed out from the quoted r 2 value of 0.81, and the approximate depths, z i of the 48 sediment traps used in the OOC, from the following relationships I obtained approximate trap depths within AE25 m from Fig. 5 in Martin et al. (1987).The resulting standard error for the OOC b value is b s b % 0:059.I tested that the estimate is not very sensitive to errors in our estimation of the OOC trap depths by adding noise with variance of 635 m 2 to the z i values and repeating the calculation.
-station sampling variability due to random errors in the sediment-trap flux measurements before attributing the variability of b b or the covariance of b b with b J 0 to any systematic differences in the processes happening in the water column at different geographical locations.The variance-covariance one computes by combining estimates b b j and b

Fig. 1 .
Fig. 1.Plot of the parameter estimates b b and b J 0 together with their 95% confidence intervals separated by region.The units of b J 0 are in mmol C m À2 day À1 .Each point was estimated separately for each station by the ordinary least squares formula.The POC flux data set compiled by Berelson (2001) was used for the estimation.The data set includes sediment-trap fluxes from the Southern Ocean (SO), Equatorial Pacific (EqPac), North Atlantic Bloom Experiment (NABE), and Arabian Sea (AS).
where r ij is the within-station error term and where b 0j denotes the logarithm of the 100 m reference POC flux for the jth station, and where b 1j denotes the b exponent for the jth station.At the second level (level-2), the model describes the across-station variability of the b coefficients (i.e. the station-to-station variability of log J 0 and b): from floating traps at the same depths.To test the influence of these two flux measurements I repeated the estimation after excluding them.Without the two shallow thorium based fluxes at the NABE station I estimated ( b b ¼ 0:68, b sd½b ¼ 0:07) and ( b b ¼ 0:68) for the full and null models (denoted full Ã and null Ã in Table

Fig. 2 .
Fig. 2. Plot of the exponent b b and export flux b J 0 estimated from the Berelson (2001) data set and the two-level mixed-effects model.The units of b J 0 are mmol C m À2 day À1 .The outlier labeled NABE corresponds to the North Atlantic Bloom Experiment.Note that b J 0 is plotted on a logarithmic scale.

Fig. 3 .
Fig.3.Same as Fig.2except that the two shallow NABE fluxes fromBuesseler et al. (1992) measured by the thorium method were excluded fromBerelson's (2001) data set for the estimation.Note that the axis scales are different than Fig.2.
t 00 and t 11 represent, respectively, the variances of log J 0 and b across stations, i.e. the ''true'' stationto-station variance of the parameters as opposed to the variance of the estimates.Similarly t 01 ¼ t 10 represents the ''true'' covariance between log J 0 and b.

Table 1
Parameters estimated from the Berelson (2001) POC flux data set for the two-level mixed-effects regression model