- Email: [email protected]
Quantification of co-occurring reaction rates in deep subseafloor sediments
Available online at www.sciencedirect.com
Geochimica et Cosmochimica Acta 72 (2008) 3479–3488 www.elsevier.com/locate/gca
Quantification of co-occurring reaction rates in deep subseafloor sediments Guizhi Wang a,*, Arthur J. Spivack a, Scott Rutherford b, Uri Manor c,1, Steven D’Hondt a a
Graduate School of Oceanography, University of Rhode Island, South Ferry Road, Narragansett, RI 02882, USA b Department of Natural Sciences, Roger Williams University, One Ferry Road, Bristol, RI 02809, USA c Parks College of Engineering, Saint Louis University, Saint Louis, MO 63130, USA Received 10 August 2007; accepted in revised form 24 April 2008
Abstract Net rates of biogeochemical reactions in subseafloor sediments can be quantified from concentration profiles of dissolved reactants or products and physical properties of the sediment. To study net rates of microbial activities in deep sediments, we developed a robust approach that is well suited to use over a broad range of sediment depths. Our approach is based on a finite-difference solution to a continuity equation that considers molecular diffusion, sediment burial, fluid advection, and reaction under the assumption of steady state. Numerical procedures are adopted to identify the maximum number of depth intervals with statistically different reaction rates. The approach explicitly considers downcore variation in physical properties and sample spacing. Uncertainties in the rate estimates are quantified using a Monte Carlo technique. We tested our approach using synthetic concentration profiles generated from analytical solutions to the continuity equation. We then applied the approach to concentration profiles of dissolved sulfate, sulfide, methane, and manganese in the 420-m thick sediment column of eastern equatorial Pacific Ocean Drilling Program Site 1226. Our results indicate that (i) sulfate reduction and iron reduction occur at most sediment depths, (ii) net methane production occurs in discrete depth intervals and (iii) manganese reduction occurs near the seafloor and deep in the sediments. These results provide quantitative evidence that multiple respiration pathways co-exist in the same depth intervals of these deep subseafloor sediments. Ó 2008 Elsevier Ltd. All rights reserved.
1. INTRODUCTION This study has two principal objectives. The first objective is to provide a rigorous technique for quantifying net rates of biogeochemical reactions throughout deep subseafloor sediment columns. The second objective is to quantify the extent to which multiple energy-yielding microbial reactions co-occur in deeply buried marine sediments. Previous studies have used concentration profiles of dissolved chemicals to infer, but not quantify, the co-occurrence of multi-
*
Corresponding author. E-mail address: [email protected] (G. Wang). 1 Present address: Department of Biology, Johns Hopkins University, 3400 N. Charles, Baltimore, MD 21218, USA. 0016-7037/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.gca.2008.04.024
ple metabolic reactions in deep subseafloor sediments, including sulfate reduction and methane production (Mitterer et al., 2001; Bralower et al., 2002; D’Hondt et al., 2002) and sulfate reduction, metal reduction and methane production (D’Hondt et al., 2004). Profiles of dissolved chemical concentrations in sedimentary pore fluids supply fundamental information about sedimentary processes. They are used to identify and quantify net rates of microbial reactions (e.g., Berner, 1964; Froelich et al., 1979; Devol and Ahmed, 1981; Schulz, 2000; D’Hondt et al., 2002, 2004; Burdige, 2006; Wortmann, 2006), net rates of diagenetic reactions, such as CaCO3 precipitation/dissolution caused by metabolic CO2 production (e.g., Emerson and Bender, 1981; Hales et al., 1994), and net rates of fluid advection through sediments and oceanic crust (Baker et al., 1991). This type of data is
3480
G. Wang et al. / Geochimica et Cosmochimica Acta 72 (2008) 3479–3488
particularly important for quantifying metabolic rates in deeply buried subseafloor ecosystems, where microbial activities are otherwise difficult to quantify (e.g., D’Hondt et al., 2004; Parkes et al., 2005). Two principal approaches have been used to estimate net reaction rates from dissolved chemical concentration data. The first approach uses these data to constrain rate parameters in an analytical solution of the continuity equation where reaction rates vary as a simple function of depth (e.g., Berner, 1964; Westrich and Berner, 1984). If reaction rates are not a simple function of depth, the usefulness of this approach is obviated. In deep subseafloor sediments, dissolved chemical profiles provide qualitative evidence that at least some reaction rates are not simple functions of depth. For example, profiles of dissolved iron and dissolved manganese concentrations indicate that net production of these dissolved chemicals occurs in discrete zones deep beneath the seafloor at some Ocean Drilling Program Sites (D’Hondt et al., 2004). The second approach uses finite-difference methods to directly calculate net reaction rates at discrete depths (e.g., Berg et al., 1998). Several studies have used a finitedifference approach to quantify reaction rates in deep subseafloor sediments (Bottcher et al., 2004; Schippers et al., 2005; Arndt et al., 2006; Wortmann, 2006). In this study, we use a finite-difference method. Our approach differs from most previous finite-difference studies of subseafloor reaction rates by explicitly allowing reaction zones of variable thickness and variable sample spacing and by accommodating vertical variations in porosity, tortuosity, pore water burial rate, and diffusion coefficient. It differs from all previous studies by quantifying uncertainties of the calculated rates, using a Monte Carlo technique. We provide a Matlab program for this approach, which users can download and, if necessary, adjust for their purposes. Briefly, our approach finds the maximum number of constant reaction rate depth intervals that is statistically significant. A Monte Carlo method (Bevington and Robinson, 2003) is applied to estimate the uncertainty of the calculated rates due to error associated with sampling and chemical measurements. In our study, the computer program is written in Matlab, but the principles are general. It consists of two sub-programs. Sub-program NRR1, provided in Electronic annex EA-1, estimates a reaction rate profile from measured concentrations and physical properties. Sub-program NRR2, provided in Electronic annex EA-2, quantifies uncertainties of the rates. The manual for the program is provided in Electronic annex EA-3. In this paper, we describe the approach, test its reliability with synthetic data, and apply it to pore fluid profiles from eastern equatorial Pacific Ocean Drilling Program (ODP) Site 1226. 2. PROCEDURES 2.1. Continuity equation For steady state conditions in one dimension, the mass balance of a pore-water chemical species is expressed as
( ) o DðxÞ o½/ðxÞCðxÞ 2 þ ½/ðxÞbðxÞ þ /ðxÞvðxÞCðxÞ ox ox h ðxÞ þ RðxÞ ¼ 0
ð1Þ
where C(x) is concentration in pore water, x is depth below the seafloor, /(x) is porosity, which is the ratio of the volume of pore water to the volume of the whole sediment (e.g., Berner, 1980), D(x) is the molecular diffusion coefficient, h2(x) is tortuosity, b(x) is pore water burial velocity, v(x) is externally driven flow velocity, R(x) is the rate of diagenetic reaction per unit volume of sediment, and (x) signifies that the parameter is a function of sediment depth. Positive values of R(x) refer to net production, while negative values indicate net consumption. To simplify the notation, in later equations, (x) will be omitted in most cases. Based on scaling analysis, the steady state approximation is reasonable for sediments where D/(Ls0) > 1 and R(x) does not vary on time scales shorter than L2/(2D), where L is the length scale of the sediment column and s0 is the sedimentation rate near the seafloor. 2.2. Diffusion coefficients Diffusion coefficients are calculated for in situ temperatures, based on the Stokes–Einstein relation (Bockris and Reddy, 1970) DðxÞ ¼
T ðxÞ gðT r Þ DðT r Þ g½T ðxÞ T r
ð2Þ
where D(Tr) is a solute’s diffusion coefficient at a reference temperature Tr (kelvin), T(x) is the temperature (kelvin) at depth x, g[T(x)] is the viscosity of the pore water at T(x), and g(Tr) is the viscosity of the pore water at Tr. A compilation of diffusion coefficients of dissolved chemicals at specific temperatures is in Schulz (2000). The viscosity of pure water at in situ temperature is calculated according to the empirical relationship in Lide (2004). We assume that the viscosity of seawater is the same as that of pure water. 2.3. Tortuosity and formation factor Tortuosity, h2, is calculated by the relation, h2 ¼ /f
ð3Þ
where f is the formation factor, which is the ratio of the electrical resistivity in sediment to the resistivity in pore water extracted from the sediment (e.g., McDuff and Ellis, 1979). Following the substitution of Eqs. (3) into (1), we have o D oð/CÞ ð/b þ /vÞC þ R ¼ 0 ð4Þ ox /f ox The parameters in this equation can be in any units as long as they are consistent. In our study, the following units are used: for D, m2 yr1; for C, mol m3; for x, meters below seafloor (mbsf); for b, m yr1; for v, m yr1; and for R, mol m3 yr1. When formation factor is not measured, it can be approximated by an empirical relationship with porosity. A commonly used relationship is f = a/n, where a and n are
Quantification of co-occurring reactions
constants. Several studies have derived the two constants for different sediment compositions from measurements of f and / (e.g., Berner, 1980; Boudreau, 1997; Breitzke, 2000). 2.4. Pore water burial rate and externally driven flow If there is negligible compaction, which is often a good approximation for shallowly buried sediments, the pore water burial rate is the same as the sediment burial rate. If, however, compaction is non-negligible, interstitial waters are displaced and pore water and sediment are buried at different rates (Berner, 1980). If steady state compaction and constant sedimentation rate are assumed, then the product, /(x)b(x), is constant (Berner, 1980; Burdige, 2006) and bðxÞ ¼
where 2 Df ð/j bj þ /0 v0 ÞDx jþ1 j AAj ¼ 2 D D ðDxÞ þ 2 f ðDxÞ2 þ ð/jþ1 bjþ1 /j bj ÞðDxÞ3 2 f jþ1 j D D 2 f 2 f þ ð/jþ1 bjþ1 þ /0 v0 ÞDx j jþ1 BBj ¼ 2 Df ðDxÞ2 þ 2 Df ðDxÞ2 þ ð/jþ1 bjþ1 /j bj ÞðDxÞ3 jþ1 j 2 Df 2 Df þ ð/j1 bj1 þ /0 v0 ÞDx j j1 þ 2 D D ðDxÞ2 þ ð/j bj /j1 bj1 ÞðDxÞ3 2 f ðDxÞ þ 2 f 2
D f
j
j1
and
/L ð1 /0 Þ s0 /ðxÞð1 /L Þ
ð5Þ
where /L is the porosity at a depth where compaction becomes negligible, and /0 is the porosity near the seafloor. Multiplying both sides of Eq. (5) by /(x), we have /ðxÞbðxÞ ¼ /L
3481
ð1 /0 Þ s0 ð1 /L Þ
ð6Þ
The right hand side of Eq. (6) is a constant, which indicates that the product of porosity and pore water burial rate is depth-independent although both can vary with depth. Additionally, assuming that pore water behaves as an incompressible fluid, the product, /(x)v(x), is also constant whether or not there is compaction. Thus vðxÞ/ðxÞ ¼ v0 /0
ð7Þ
where v0 is the flow rate near the seafloor. 2.5. Numerical determination of reaction rates Reaction rates are estimated numerically by the following steps: (a) data filtering; (b) numerical approximation of the reaction rate profile; (c) optimization of the number of reaction zones; and (d) estimation of uncertainty. (a) Data filtering: To minimize the effect of outliers, C, f, and / are first smoothed using a Gaussian 5-point filter. (b) Numerical approximation of the reaction rate profile: To solve for R(x), the sediment column is divided into n layers, each with the same length, the boundaries of each corresponding to discrete depth points (x0, x1, . . ., xn). The thickness of each layer is Dx. C, D, f, /, and b are interpolated to these depths. At depth xj (j = 1, . . ., n 1), Eq. (4) is approximated using a centered finite-difference method as h i h i D oð/CÞ D oð/CÞ ð/b þ /vÞC ð/b þ /vÞC /f /f ox ox jþ1=2
j1=2
Dx
j
ð8Þ Following finite-difference approximations of the two terms of the numerator [see Wang (2006) for details] and substituting /(x)v(x) with /0v0 based on Eq. (7), Eq. (8) is converted into a set of equations ð9Þ
j1
Therefore, for the sediment column there are (n 2) equations of C vs. R, which can be expressed in matrix notation: AC ¼ S where
BB1 AA1 AA2 CC 2 BB2 A ¼ CC n2 BBn2 CC n1 R1 CC 1 C 0 R2 S ¼ Rn2 Rn1 AAn1 C n
ð10Þ ; AAn2 BBn1
and
To solve for R, a least-square minimization method is used. R(x) is assumed to be a piecewise constant function. The sediment column is divided into N zones with constant reaction rates. The length of each zone is determined by a specified number of chemical measurements. In other words, the reaction rate profile over the depth range 0 to L is described by N constant rates. From this rate profile, we generate [R1, . . ., Rn1] corresponding to depths [x1, . . ., xn1]. The smoothed measured concentrations (Cm) at the top and bottom of the sediment column are taken as the boundary conditions so that C0 and Cn are known. The concentration estimated by [R1, . . ., Rn1] is ^ ¼ A1 S C
¼ Rj
AAj C jþ1 þ BBj C j þ CC j C j1 ¼ Rj
2 Df þ ð/j bj þ /0 v0 ÞDx j1 j CC j ¼ 2 D D ðDxÞ2 þ ð/j bj /j1 bj1 ÞðDxÞ3 2 f ðDxÞ þ 2 f 2 Df
where C1 C2 ^ ¼ C C n2 C n1
ð11Þ
3482
G. Wang et al. / Geochimica et Cosmochimica Acta 72 (2008) 3479–3488
We used the Thomas algorithm (Von Rosenberg, 1969) ^ in our Matlab program as A is a tri-diagonal to solve for C ^ is then interpolated to the depths of the measurematrix. C ments and the sum of squared errors (SSE) between Cm and ^ is calculated. SSE is then minimized to determine C [R1, . . ., Rn1] (and hence a reaction rate profile with values of N rates) that provides a best-fit concentration profile to the measured concentration profile. We used Fminsearch, a multi-dimensional minimization function in Matlab that uses the simplex search method of Lagarias et al. (1998), for this purpose, though other minimization techniques could be used. The minimization procedure we wrote in Matlab is provided in Electronic annex EA-4. For a total number of concentration data points, M, each at a different depth, the maximum number of reaction zones is equal to (M 2) because at least three points are needed to approximate a reaction rate. Using (M 2), however, most likely will result in relatively noisy estimated rates, of which many are not statistically significant, as is the case when reaction rates are directly calculated using a finite-difference method at discrete depths. It also results in long runtimes of the program. From our tests with synthetic concentration profiles, M/3 is a reasonable compromise, for the maximum number of zones, between runtime and potential number of reaction zones that are statistically significant. That is, at least three concentration data points, at different depths, must be available within each reaction zone to determine a constant rate. In this portion of the model, N is initially set at 1 and then increased to the maximum number of reaction zones, e.g., M/3. The minimization is repeated for each N, and subsequently, a set of rate profiles and corresponding concentration profiles are generated. (c) Optimization of the number of reaction zones: Two steps are used to optimize the reaction rate profile. First the largest number of statistically significant reaction zones is determined from the set of rate profiles generated above, using partial F-tests. Second, adjacent zones with similar rates are combined into zones with constant rates if combining does not decrease the statistical significance. For the first step, the null hypothesis tested here is, given a number of reaction zones, K (i.e., specific number of variables), the addition of more zones (i.e., more variables) does not significantly improve the fit of modeled concentrations to the measured concentrations. The F value for the hypothesis is calculated as (Kleinbaum and Kupper, 1978): F ¼
SSEK SSEL LK SSEL ML1
ð12Þ
where L is a number of reaction zones (L > K). A probability is calculated for each F value and compared with the significance level the user chooses (0.05 is commonly used) to accept or reject the hypothesis. In the second step, adjacent zones with the most similar rates are combined and then the minimization is repeated. The merging of adjacent zones and re-determination of reaction rates are repeated until there is only one zone for the entire data set. The new set of fit profiles is compared again using partial F-tests to determine a final rate profile.
(d) Uncertainty estimates: A Monte Carlo method is used because the uncertainty limits in the estimates of R due to analytical and sampling uncertainty cannot be expressed in a simple formula. The Monte Carlo method we use in NRR2 makes use of randomly determined numbers to simulate processes that involve an element of chance, as there is in the errors associated with analysis and sampling. In our application, for each sample, random chemical measurements (Cr) are assumed to be described by a Gaussian distribution, with a mean of Cm and a relative deviation due to uncertainties associated with chemical measurements and sampling (e). A random measurement can then be expressed as: C r ¼ C m ð1 þ e rÞ
ð13Þ
where r is a random number generated from the standard Gaussian distribution N(0, 1). In our approach, a set of randomized concentration profiles is generated, based on the smoothed original (measured) concentration profile, Cm. A rate profile is estimated for each randomized concentration profile, using the same procedures and the same values of other parameters as in NRR1. A standard deviation is calculated at each depth of the estimated rates from this set of rate profiles, taking the rates calculated from the original concentration profile as the mean. Using a Monte Carlo method to determine the uncertainty of estimated rates, the greater the number of random concentration profiles generated, the closer the calculated uncertainty is to the true uncertainty. However, the more profiles generated, the longer it takes to run NRR2. From our tests, generation of 50 random profiles provides a good compromise between length of runtime and uncertainties. 3. ASSESSMENT 3.1. Reliability test The approach was tested following three steps: (i) creating a concentration profile with an analytical solution and a specified rate profile (hereafter referred to as analytical rates), (ii) determination of rates from the analytical concentration profile using the program, and (iii) comparison of the model-derived rates with the analytical rates. We tested the approach with constant analytical rates, variable analytical rates, and varied sample spacing. In all of these tests (Table 1), the approach reproduced the analytical rates very well (Fig. 1). As an illustrative example, we present Test 1, which is a test with a relatively complicated reaction rate profile where there is an increasing rate of depletion of a dissolved species with depth and then an abrupt decrease in rate to zero at 50 m. This test works well in demonstrating the strength of our approach in identifying both smoothly varying and abrupt changes of rates. We parameterized the consumption rate, under the assumption that it increases exponentially with depth, RðxÞ ¼ Rð0ÞeðaxÞ
ð14Þ
until the concentration reaches 0 at 50 mbsf below which the rate is zero, where R(0) is the consumption rate at
Quantification of co-occurring reactions
3483
Table 1 Tests of the numerical method Parameters R (mol m3 yr1) 2
1
D (m yr ) f / v (m yr1) b (m yr1) BCa: C (mol m3), x (m)
Test 2
Test 3
Test 4
1.5 104
1.5 104
0
0.018 2.2 0.8 1.0 105 2.5 105 C(x = 1) = 30 C(x = 50) = 0
0.018 2.2 0.8 0 0 C(x = 1) = 30 C(x = 100) = 25
0.018 2.2 0.8 1.0 105 2.5 105 C(x = 1) = 30 C(x = 100) = 25
0.018 2.2e0.002x 0.8 0 0 C(x = 1) = 30 C(x = 200) = 25
BC is the abbreviation for boundary conditions.
a
0
Depth (mbsf)
a
Test 1 1:5 104 e0:01x ; 0 6 x 6 50 x > 50 0;
20
b
analytical concentrations modeled rates
40
analytical rates
60
1 σ uncertainty envelope
80 100 0
Depth (mbsf)
c
10 20 30 -3 -2 -1 0 Concentration Reaction rate (mol m -3 ) (10-4 mol m-3 yr -1)
-2 -1 10 20 30 0 Reaction rate Concentration (10-4 mol m-3 yr -1) (mol m-3 )
0
d
0 20
0
50
40 100 60 150
80
200
100 0
10 20 30 0 -2 -1 Concentration Reaction rate (mol m-3 ) (10 -5 mol m-3 yr -1)
0
10 20 30 -2 0 2 Concentration Reaction rate (mol m-3 ) (10-5 mol m-3 yr -1)
Fig. 1. Method assessment results. Rates and uncertainties estimated by the numerical method are compared to the specified analytical rate profiles, which are described in detail in Table 1. (a) Test 1: variable reaction rates with depth, (b) Test 2: constant reaction rate with depth and negligible pore water burial rate and externally driven flow rate (i.e., b, v = 0), (c) Test 3: constant reaction rate with depth, and pore water burial rate and externally driven flow not equal to 0, and (4) Test 4: no reaction (i.e., R = 0), b and v equal to 0, and f changing with depth.
x = 0, and a is a constant. Then Eq. (4) was solved analytically assuming D, f, and / are constant with depth /f ðvþbÞ x D
C ¼P þQe
D
a2
Rð0Þ f eax / f ðv þ bÞ a
ð15Þ
where P and Q are constants. From this analytical solution, we generated a concentration profile for a 100-m thick sediment column. The sediment depth evaluated started at 1 mbsf. The concentration at the top of the sediment column was set to 30 mol m3, D = 1.80 102 m2 yr1, f = 2.2, / = 0.8, v = 1.0 105 m yr1, 5 1 b = 2.5 10 m yr , R(0) = 1.50 104 mol m3 yr1, and a = 0.01 m1. P and Q were solved from the boundary
conditions in Table 1. A variable depth increment was used to generate a concentration profile that has sample spacing similar to real measurements. Along with this analytical concentration profile, constant D, f, /, v, and b were input into NRR1. The chosen minimum number of concentration data points in each reaction zone was 3, since the analytical concentration profile is smooth (Fig. 1a), and a significance level of 0.05 was used. NRR1 was run with a maximum of 12 zones and identified 7 statistically significant reaction zones. In the subsequent run of NRR2, 50 random concentration profiles were generated, in which a relative standard deviation in concentration of 1% was assumed.
3484
G. Wang et al. / Geochimica et Cosmochimica Acta 72 (2008) 3479–3488
The results, compared with the analytical rates, are shown in Fig. 1a. The estimated rates clearly identify and reproduce the exponential increase in magnitude in the analytical rate. Moreover, the estimated rates explicitly reflect the abrupt change of the analytical rates near 50 mbsf. Additionally, the 1r uncertainty envelope closely brackets the analytical rates. In each reaction zone, the estimated constant rate is approximately equal to the mean of the analytical rates. The results of this test indicate that our approach performs well in identifying variable and abruptly changing rates from a concentration profile and in providing an accurate estimate of rates with objective uncertainty limits. 3.2. Application of the approach In most deeply buried sediments, metabolic activities are too low to be quantified with radiotracer methods. Additionally, low spatial resolution in sampling has previously limited the applicability of pore fluid modeling in deep sediments. ODP Site 1226 provides high-resolution sampling of a deep sediment column, which makes it ideal for the application of our method. At this site, co-occurrence of multiple metabolic activities is suggested from qualitative analysis of chemical profiles (D’Hondt et al., 2004). The application of the method described in this paper to data from Site 1226 illustrates the power of the method for quantifying metabolic activity. These results provide quantitative evidence of the co-existence of multiple metabolic activities through a deep sediment column and can be used to test basic aspects of metabolism in this environment. In particular, this work quantifies rates of sulfate, iron, and manganese reduction along with methane generation as a function of depth. Site 1226 is located in the eastern equatorial Pacific at a water depth of 3297 m (Shipboard Scientific Party, 2003). The sediment column at this site is 420-m thick. Seawater flows through the underlying basalt driving sulfate back up toward seawater concentrations at the bottom of the hole (D’Hondt et al., 2003). High-resolution measurements are available for dissolved sulfate, total dissolved sulfide, dissolved methane, dissolved manganese, porosity, and formation factor throughout the sediment column (Fig. 2) (Shipboard Scientific Party, 2003). In situ down-hole temperature (t, °C) is given by: t = 1.74 + 0.0572x (Shipboard Scientific Party, 2003). The sedimentation rate at Site 846 (Shipboard Scientific Party, 1992), which is within 100 m of Site 1226, was taken as the sedimentation rate at Site 1226. The effect of compaction on the pore water burial rate was considered in our calculation. There is no identifiable advection at this site (Shipboard Scientific Party, 2003). A significance level of 0.05 was used in the calculation of the significance of reaction rates. The combined analytical and sampling uncertainty for each chemical species was estimated from the deviations of the measured profiles from their smoothed profiles. These are 1%, 10%, 7%, and 1% for sulfate, sulfide, methane, and manganese, respectively. For each chemical species, 50 random concentration profiles were generated by NRR2 to estimate uncertainty by the Monte Carlo method.
The number of statistically significant reaction intervals varied between chemical species. For dissolved sulfate (Fig. 2c) and methane (Fig. 2d), there were 5 intervals. For dissolved sulfide (Fig. 2e) and manganese (Fig. 2f), there were 6 and 14 intervals, respectively. The confidence limits were 6%, 22%, 15%, and 75% for sulfate, sulfide, methane, and manganese, respectively. The concentrations predicted from the calculated rates fit the observed data very well. Overall, sulfate is the dominant terminal electron acceptor in the sediment column although in some intervals, its reduction rate is below the limit of detection. The depthintegrated net sulfate reduction rate from 3 to 418 mbsf is 1.8 103 ± 1 104 mol m2 yr1. The depth-integrated net rate of dissolved sulfide production from 0.59 to 410 mbsf is 7.6 105 ± 1.7 105 mol m2 yr1. The depth-integrated net production rate of dissolved manganese from 0.59 to 418 mbsf is 1.2 105 ± 9 106 mol m2 yr1. This negative value is almost certainly a consequence of beginning the Mn reaction rate modeling at 0.59 mbsf; Mn produced at shallower depths migrates downward past this depth to be precipitated below (Fig. 2f). Unfortunately, there are no dissolved Mn data for depths between 0 and 0.59 mbsf (Shipboard Scientific Party, 2003). The depth-integrated net production rate of dissolved methane from 11.9 to 384 mbsf is 3.4 107 ± 5 108 mol m2 yr1. Sulfate reduction rates are greatest in the uppermost 40 m, but the rates do not decrease smoothly with depth as is often assumed. Sulfate reduction occurs in the upper sediment column (above 40 mbsf) at a rate of 2.5 105 mol m3 yr1, in the intervals between 55 and 131 mbsf at a rate of 4.5 106 mol m3 yr1, between 198 and 251 mbsf at a rate of about 5.7 106 mol m3 yr1, and from 275 mbsf to the bottom of the column at a rate of 6.6 107 mol m3 yr1. For the other reaction zones (from 40 to 55 mbsf, from 131 to 198 mbsf, and from 251 to 275 mbsf), the rates are indistinguishable from 0. The depths of two of the sulfate reduction maxima (one near the seafloor around 40 mbsf and the other around 250 mbsf) correspond to the depths of maximum rates determined by radiotracer measurements at this site and approximately correspond to the chemical interfaces that have been suggested to stimulate microbial activities (Parkes et al., 2005). Sulfide, methane, and manganese profiles represent the net difference between production and consumption. Sulfide is produced via sulfate reduction and consumed by precipitation into sulfide minerals. There is net consumption near the seafloor (from 0.21 to 9 mbsf) at a rate of 3.4 105 mol m3 yr1 and from 256 to 286 mbsf at a rate of 2.1 ;106 mol m3 yr1. Net production occurs in the intervals from 9 to 36 mbsf at a rate of 1.3 105 mol m3 yr1, and from 36 to 178 mbsf at a rate about one order of magnitude smaller than the rate in the overlying interval (Fig. 2e). This interval from 36 to 178 mbsf covers some of the depths (between 40 and 55 mbsf and between 131 and 178 mbsf) at which rates of sulfate reduction are too small to be quantified from the sulfate profile, indicating that sulfate reduction still occurs
Quantification of co-occurring reactions
a
3485
c
b
0 Depth (mbsf)
input measurement data 100
modeled rates and concentrations
200
1σ uncertainty envelope
300 400 0.6 0.7 0.8 0.9 1.5 2.5 3.5 4.5 18 Porosity Formation factor
Depth (mbsf)
0
22 26 30 Sulfate (mol m-3)
-2 -1 0 Reaction rate -5 (10 mol m-3 yr -1) -3
e
d
100 150
200
250
300
350 400 0
1 2 Methane -3 (10 mol m-3)
-1
0 1 Reaction rate -8 (10 molm-3 yr -1)
0
0.2 0.4 0.6 Sulfide (mol m-3)
-3
-2 0 2 Reaction rate (10-5 mol m-3 yr -1)
450
0 2 -2 0 0.2 0.4 0.6 Reaction rate Sulfide -6 (mol m-3) (10 mol m-3 yr -1)
f 0
Depth (mbsf)
0
10
100
20 30
200
40
0 2 4 -8 -4 Manganese Reaction rate (10-2 mol m-3) (10-6 mol m-3 yr -1 )
0
300 400 0
2 4 -8 -4 0 Manganese Reaction rate (10-2 mol m-3) (10-6 mol m-3 yr -1)
Fig. 2. Calculated reaction rates based on concentrations and physical properties at ODP Site 1226. (a) Porosity, (b) formation factor, (c) dissolved sulfate, (d) dissolved methane, (e) total dissolved sulfide, and (f) dissolved manganese. Physical properties and dissolved chemical concentrations are from Shipboard Scientific Party (2003).
in these intervals at rates at least the same as the net rate of sulfide production. From 178 to 256 mbsf, net reaction rate of sulfide is indistinguishable from 0. From 286 mbsf to the bottom of the column, dissolved sulfide is below the detection limit (0.17 lM) and we take sulfide concentrations in this interval to be 0 (Fig. 2e). Modeled net sulfide reaction rates are close to 0 in this interval (Fig. 2e). The distributions of the net rates of dissolved sulfate consumption and dissolved sulfide production and consumption suggest that sulfate reduction occurs at most, if not all, of the sediment depths (rates are indistinguishable from 0 in the interval from 178 to 198 mbsf). Dissolved sulfide produced by sulfate reduction is precipitated almost completely throughout the sediment column except in the intervals from 40 to 55 mbsf and from 131 to 178 mbsf, and from 251 to 256 mbsf, where there is no quantitative evidence of net consumption of dissolved sulfide. Based
on the depth-integrated rates, about 96% of the dissolved sulfide produced by sulfate reduction is precipitated in the sediment column. The most likely ultimate sink for dissolved hydrogen sulfide is pyrite (FeS2), which is present throughout the sediments at this site (Shipboard Scientific Party, 2003) and is common in anoxic sediments. This local precipitation of iron sulfide implies that reduced iron [Fe(II)] must be produced in depth intervals where sulfate reduction occurs because there are almost no gradients of dissolved Fe(II) concentrations throughout the sediment column in close proximity to where sulfate reduction occurs (Fig. 3). Considering that concentrations of dissolved Fe(II) remain low in the sediment column, the rates of iron reduction and sulfate reduction seem to be balanced locally. For example, assuming pyrite is the major sink of dissolved Fe(II) and dissolved sulfide, the rate of iron reduction
3486
G. Wang et al. / Geochimica et Cosmochimica Acta 72 (2008) 3479–3488 0
Depth (mbsf)
100
200
300
400
0
0.02 0.04 Iron (mol m-3)
Fig. 3. The concentration profile of dissolved iron in the sediments of Site 1226 (Shipboard Scientific Party, 2003).
would be about half of the rate of sulfate reduction. In this case, the depth-integrated rate of iron reduction is 8.6 104 ± 5 105 mol m2 yr1. In anoxic marine sediments, methane is produced by microbial methanogenesis and consumed by microbially mediated anaerobic methane oxidation using sulfate as the terminal electron acceptor (Devol and Ahmed, 1981) as well as by adsorption onto sediments (Whiticar and Suess, 1990; Knies et al., 2004). At Site 1226, net methane production occurs in the intervals between 50 and 107 mbsf at a rate of 1.4 108 mol m3 yr1 and from 164 to 279 mbsf at a smaller rate of 2.2 109 mol m3 yr1 (Fig. 2d). Methane is consumed and/or adsorbed in the intervals from 3.2 to 50 mbsf at a rate of 4.1 109 mol m3 yr1, from 107 to 164 mbsf at a rate of 6.4 109 mol m3 yr1, and from 280 to 384 mbsf at a rate of 1.6 109 mol m3 yr1. Net dissolved manganese production occurs near the seafloor: from 4.5 to 15 mbsf at a rate of 9.0 107 mol m3 yr1 and from 40 to 78 mbsf at a much smaller rate of 5.7 108 mol m3 yr1 (Fig. 2f). It also occurs deep in the sediment column between 272 and 313 mbsf at a rate of about 3.0 107 mol m3 yr1 and from 391 mbsf to the bottom of the column at a rate of 1.2 107 mol m3 yr1. At other depths, dissolved manganese is consumed (e.g., due to precipitation into MnCO3 or in solid solution with CaCO3) or the rate is indistinguishable from 0. Usually it is inferred that net manganese reduction occurs at a rate at least the same as the rate of net dissolved manganese production (e.g., D’Hondt et al., 2004). In theory, some dissolved manganese may be released from solid phases when they become undersaturated. However, there is no indication from other solutes such as Ca2+ and alkalinity (Shipboard Scientific Party, 2003) that the local production of dissolved Mn is from dissolution. Thus, it is most likely due to manganese reduction. Except in the intervals from 178 to 198 mbsf and from 251 to 256 mbsf, zones of net dissolved methane production correspond to zones where sulfate reduction clearly occurs,
indicating that sulfate reduction does not result in the competitive exclusion of methanogenesis in these intervals. Net dissolved manganese production occurs in intervals near the seafloor and deep in the sediment column, where sulfate reduction clearly occurs. Our results quantitatively confirm the co-occurrence of sulfate reduction, net methanogenesis, and net manganese production at this site and indicate that iron reduction occurs in the depth intervals where sulfate reduction occurs. Besides microbial reduction of iron and manganese, abiotic reduction of iron and manganese by sulfide may co-occur in these sediments as observed in sediment incubations (Canfield, 1989; Koretsky et al., 2003). Our results are in agreement with the qualitative interpretation of D’Hondt et al. (2004). To test the effect of downcore variations of porosity and formation factor on the modeled reaction rates, we calculated sulfate reaction rates at Site 1226 using the downcore average formation factor and porosity. Rates calculated with constant formation factor and porosity do not resolve spatial variations in rates that are identified when measured formation factors and porosities are used. The estimated rate profile is consisted of 2 reaction zones compared with the 5 reaction zones estimated using measured formation factor and porosity. In addition, to examine the effect of compaction on calculated rates, we estimated sulfate reaction rate at this site using measured formation factor and porosity without considering compaction. The resulting rate profile has 6 reaction zones with rates about 16% different from the rates estimated considering compaction. More importantly, the uncertainties of the rates are about twice those when compaction is considered. Results of the application of the approach to these data demonstrate that our approach is an effective tool for quantifying metabolic activities in deeply buried sediments from profiles of dissolved chemicals and sediment physical properties. Our approach is also a useful tool in studying biogeochemical cycles in deep sediments. For example, our Site 1226 results indicate that removal of dissolved sulfide from the pore water into solid phases occurs at most sediment depths and this sediment column is a significant sink of sulfur. In addition, our approach provides a realistic description of metabolic reactions in deeply buried sediments where direct radiotracer measurements are too insensitive or problematic due to difficult sample manipulations (Jørgensen et al., 2001). 4. COMMENTS AND RECOMMENDATIONS ON THE METHOD Our method estimates reaction rates without assuming any a priori relationship of rate to sediment depth and incorporates variations in physical properties that affect pore-water chemical gradients and inferred reaction rates. It also successfully deals with variable spacing of data and changes in sediment physical properties. In comparison, the use of equal length for each zone (e.g., Berg et al., 1998) can overestimate rates in intervals with few concentration data or miss reaction zones in relatively narrow intervals with apparent chemical gradients. Gradients of chemical profiles caused by changes of physical proper-
Quantification of co-occurring reactions
ties are considered by the incorporation of in situ distributions of physical properties. Using the numerical procedures, the approach provides objective judgments on where and how much a dissolved chemical is produced or consumed. The uncertainties in estimated reaction rates that are introduced due to the analytical error associated with chemical measurements are well constrained using a Monte Carlo simulation of concentrations. Although our approach effectively quantifies net rates of respiration reactions in deeply buried sediments, it has limitations. The first limitation is the assumption of constant rates in each reaction zone. Rates within a zone may vary on scales too fine to resolve; the constant rate calculated for a reaction zone is the approximate average of the rates within the zone. The second limitation is, as for all numerical approximations, that the rates at the boundary of a column cannot be predicted without additional inputs. For example, rates cannot be calculated at the seafloor and the sediment/basement interface. Third, the user determines the minimum number of concentration data points in each reaction zone, which is somewhat subjective. This determination may require a trial-and-error process because the user will not want to overestimate rates by using too small a number, nor will the user want to miss apparent reaction zones by using too large a number. Our approach accounts for variations in temperature, diffusivity, formation factor, porosity, viscosity, externally driven flow velocity, and pore water burial rate with depth to reflect the complexity of deep sedimentary environments. However, users have the option of using constant values for these parameters to fit their study environments. If a user wants more smoothing of his/her measurement data than 5-point Gaussian smoothing, locally weighted regression scatter plot smoothing (Lowess) (Chambers et al., 1983) is recommended. Lowess smoothing can be found in a Matlab toolbox and other commercial software. The user can smooth the data profiles using Lowess before loading them into the program, or embed Lowess in NRR1 if the toolbox is available. For very noisy chemical data, the uncertainty of the depth-integrated rate may be 100% of the rate or even higher. In this case, the measured data may not define profiles clean enough for the approach to identify reaction zones. The mass-balance equation that we used makes this approach applicable in environments where advection may be one of the dominant transport processes. Mathematically, externally driven flow is similar to the pore water burial process as can be seen in Eq. (4) and the model handles the external flow rate almost the same way as the pore water burial rate. For radioactive chemicals, the user can incorporate radioactive decay into Eq. (4) and follow the same procedures to estimate reaction rates. The user must keep in mind that the approach applies to dissolved species approaching steady state. 5. CONCLUSIONS Our approach successfully identifies variable reaction rates from sedimentary pore fluid concentration profiles and provides a good approximation of rates with con-
3487
strained uncertainties. It is an effective tool for using profiles of dissolved chemicals and physical properties to study metabolic activities in deeply buried sediments. Application of our approach to data from ODP Site 1226 provides quantitative evidence of the depth distributions of microbial processes deep beneath the seafloor. These results quantitatively demonstrate that sulfate reduction and methanogenesis occur in the same reaction zones deep in a subseafloor sediment column, as do sulfate reduction and dissolved manganese production. This sediment column is a sink for sulfur, where the dissolved sulfide produced from sulfate reduction is precipitated almost completely into solid phases in the depth intervals where sulfate reduction occurs. The existence of this sulfide sink suggests that iron reduction may also occur throughout the sediment column in the reaction zones where sulfate reduction occurs. ACKNOWLEDGMENTS This manuscript was greatly improved by the thoughtful comments of Andrew Dale, Timothy Ferdelman, David Burdige, and an anonymous reviewer. We thank Andrew Dale and Heather Schrum for independently testing the Matlab programs and Heather Schrum for reviewing the manual for the Matlab programs. This study was supported by the NASA Astrobiology Institute.
APPENDIX A. SUPPLEMENTARY DATA Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.gca.2008. 04.024. REFERENCES Arndt S., Brumsack H. and Wirtz K. W. (2006) Cretaceous black shales as active bioreactors: a biogeochemical model for the deep biosphere encountered during ODP Leg 207 (Demerara Rise). Geochim. Cosmochim. Acta. 70, 408–425. Baker P. A., Stout P. M., Kastner M. and Elderfield H. (1991) Largescale lateral advection of seawater through oceanic crust in the central equatorial Pacific. Earth Planet. Sci. Lett. 105, 522–533. Berg P., Risgaard-Petersen N. and Rysgaard S. (1998) Interpretation of measured concentration profiles in sediment pore water. Limnol. Oceanogr. 43, 1500–1510. Berner R. A. (1964) An idealized model of dissolved sulfate distribution in recent sediments. Geochim. Cosmochim. Acta 28, 1497–1503. Berner R. A. (1980) Early Diagenesis. Princeton University Press, Princeton. Bevington P. R. and Robinson D. K. (2003) Data Reduction and Error Analysis for the Physical Sciences, third ed. McGraw Hill. Bockris J. O’M. and Reddy A. K. N. (1970). Modern Electrochemistry, vol. 1. Plenum Press, New York. Bottcher M. E., Khim B. K., Suzuki A., Gehre M., Wortmann U. G. and Brumsack H. J. (2004) Microbial sulfate reduction in deep sediments of the Southwest Pacific (ODP Leg 181, Sites 1119–1125): evidence from stable sulfur isotope fractionation and pore water modeling. Mar. Geol. 205, 249–260. Boudreau B. P. (1997) Diagenetic Models and their Implementation: Modeling Transport and Reactions in Aquatic Sediments. Springer-Verlag, Berlin.
3488
G. Wang et al. / Geochimica et Cosmochimica Acta 72 (2008) 3479–3488
Bralower T. J., Premoli Silva I., Malone M. J., et al. (2002). Proc. ODP, Init. Repts. 198 (Online). Available from: http://wwwodp.tamu.edu/publications/198_IR/198ir.htm, accessed on August 8, 2006. Breitzke M. (2000) Physical properties of marine sediments. In Marine Geochemistry (eds. H. D. Schulz and M. Zabel). Springer-Verlag, Berlin, pp. 29–83. Burdige D. J. (2006) Geochemistry of Marine Sediments. Princeton University Press, Princeton. Canfield D. E. (1989) Reactive iron in marine sediments. Geochim. Cosmochim. Acta 53, 619–632. Chambers J. M., Cleveland W. S., Kleiner B. and Tukey P. A. (1983) Graphing Methods for Data Analysis. Wadsworth International Group, Belmont, California, and Duxbury Press, Boston. D’Hondt S., Jørgensen B. B., Miller D. J., Batzke A., Blake R., Cragg B. A., Cypionka H., Dickens G. R., Ferdelman T., Hinrichs K. U., Holm N. G., Mitterer R., Spivack A., Wang G., Bekins B., Engelen B., Ford K., Gettemy G., Rutherford S. D., Sass H., Skilbeck C. G., Aiello I. W., Gue`rin G., House C., Inagaki F., Meister P., Naehr T., Niitsuma S., Parkes R. J., Schippers A., Smith D. C., Teske A., Wiegel J., Padilla C. N. and Acosta J. L. S. (2004) Distributions of microbial activities in deep subseafloor sediments. Science 306, 2216–2221. D’Hondt S., Jørgensen B. B., Miller D. J., et al. (2003) Proc. ODP Init. Repts. 201 [CD-ROM]. Available from: Ocean Drilling Program, Texas A&M University, College Station, TX 778459547, USA. D’Hondt S., Rutherford S. and Spivack A. J. (2002) Metabolic activity of subsurface life in deep-sea sediments. Science 295, 2067–2070. Devol A. H. and Ahmed S. I. (1981) Are high rates of sulfate reduction associated with anaerobic oxidation of methane? Nature 291, 407–408. Emerson S. and Bender M. (1981) Carbon fluxes at the sedimentwater interface of the deep-sea: calcium carbonate preservation. J. Mar. Res. 39, 139–162. Froelich P. N., Klinkhammer G. P., Bender M. L., Luedtke N. A., Heath G. R., Cullen D., Dauphin P., Hammond D. and Hartman B. (1979) Early oxidation of organic matter in pelagic sediments of the eastern equatorial Atlantic: suboxic diagenesis. Geochim. Cosmochim. Acta 43, 1075–1090. Hales B., Emerson S. and Archer D. (1994) Respiration and dissolution in the sediments of the western North Atlantic: estimates from models of in situ microelectrode measurements of porewater oxygen and pH. Deep-Sea Res. I I41, 695–719. Jørgensen B. B., Weber A. and Zopfi J. (2001) Sulfate reduction and anaerobic methane oxidation in Black Sea sediments. DeepSea Res. I 48, 2097–2120. Kleinbaum D. G. and Kupper L. L. (1978) Applied Regression Analysis and Other Multivariable Methods. Duxbury Press, North Scituate, Mass. Knies J., Damm E., Gutt J., Mann U. and Pinturier L. (2004) Near-surface hydrocarbon anomalies in shelf sediments off Spitsbergen: evidences for past seepages. Geochem. Geophys. Geosyst. 5. doi:10.1029/2003GC000687.
Koretsky C. M., Moore C. M., Lowe K. L., Meile C., Dichristina T. J. and Cappellen P. V. (2003) Seasonal oscillation of microbial iron and sulfate reduction in saltmarsh sediments (Sapelo Island, GA, USA). Biogeochemistry 64, 179–203. Lagarias J. C., Reeds J. A., Wright M. H. and Wright P. E. (1998) Convergence properties of the Nelder-Mead simplex method in low dimensions. SIAM J. Optim. 9, 112–147. Lide D. R. (ed.) (2004) CRC Handbook of Chemistry and Physics, 84th ed. CRC press, 2003–2004. McDuff R. E. and Ellis R. A. (1979) Determining diffusion coefficients in marine sediments: a laboratory study of the validity of resistivity techniques. Am. J. Sci. 279, 666–675. Mitterer R. M., Malone M. J., Goodfriend G. A., Swart P. K., Wortmann U. G., Logan G. A., Feary D. A. and Hine A. C. (2001) Co-Generation of hydrogen sulfide and methane in marine carbonate sediments. Geophys. Res. Lett. 28(20), 3931– 3934. Parkes R. J., Webster G., Cragg B. A., Weightman A. J., Newberry C. J., Ferdelman T. G., Kallmeyer J., Jørgensen B. B., Aiello I. W. and Fry J. C. (2005) Deep sub-seafloor prokaryotes stimulated at interfaces over geological time. Nature 436, 390– 394. Schippers A., Neretin L. N., Kallmeyer J., Ferdelman T. G., Cragg B. A., Parkes R. J. and Jørgensen B. B. (2005) Prokaryotic cells of the deep sub-seafloor biosphere identified as living bacteria. Nature 433, 861–864. Schulz H. D. (2000) Quantification of early diagenesis: dissolved constituents in marine pore water. In Marine Geochemistry (eds. H. D. Schulz and M. Zabel). Springer-Verlag, Berlin, pp. 85– 128. Shipboard Scientific Party (1992) Site 846. In Proc. ODP Init. Repts. 138 (Pt. 1) (eds. L. Mayer, N. Pisias, T. Janecek, et al.). Ocean Drilling Program, College Station, TX, pp. 265–333. Shipboard Scientific Party (2003) Site 1226. In Proc. ODP Init. Repts. 201 [CD-ROM] (eds. S. L. D’Hondt, B. B. Jørgensen, D. J. Miller, et al.). Available from: Ocean Drilling Program, Texas A&M University, College Station, TX 77845-9547, USA. Von Rosenberg D. U. (1969) Methods for the Numerical Solution of Partial Differential Equations. American Elsevier Pub. Co., New York. Wang G. (2006) Metabolic Activities in Deep Subseafloor Sediments. Ph.D. thesis. University of Rhode Island. Westrich J. T. and Berner R. A. (1984) The role of sedimentary organic matter in bacterial sulfate reduction: the G model tested. Limnol. Oceanogr. 29, 236–249. Whiticar M. J. and Suess E. (1990) Characterization of sorbed volatile hydrocarbons from the Peru margin, Leg 112, Site 679, 680/681, 682, 684, and 686/687. Proc. ODP Sci. Res. 112. doi:10.2973/odp.proc.sr.112.165.1990. Wortmann U. G. (2006) A 300 m long depth profile of metabolic activity of sulfate-reducing bacteria in the continental margin sediments of South Australia (ODP Site 1130) derived from inverse reaction-transport modeling. Geochem. Geophys. Geosyst. 7. doi:10.1029/2005GC001143. Associate editor: David J. Burdige
Geochimica et Cosmochimica Acta 72 (2008) 3479–3488 www.elsevier.com/locate/gca
Quantification of co-occurring reaction rates in deep subseafloor sediments Guizhi Wang a,*, Arthur J. Spivack a, Scott Rutherford b, Uri Manor c,1, Steven D’Hondt a a
Graduate School of Oceanography, University of Rhode Island, South Ferry Road, Narragansett, RI 02882, USA b Department of Natural Sciences, Roger Williams University, One Ferry Road, Bristol, RI 02809, USA c Parks College of Engineering, Saint Louis University, Saint Louis, MO 63130, USA Received 10 August 2007; accepted in revised form 24 April 2008
Abstract Net rates of biogeochemical reactions in subseafloor sediments can be quantified from concentration profiles of dissolved reactants or products and physical properties of the sediment. To study net rates of microbial activities in deep sediments, we developed a robust approach that is well suited to use over a broad range of sediment depths. Our approach is based on a finite-difference solution to a continuity equation that considers molecular diffusion, sediment burial, fluid advection, and reaction under the assumption of steady state. Numerical procedures are adopted to identify the maximum number of depth intervals with statistically different reaction rates. The approach explicitly considers downcore variation in physical properties and sample spacing. Uncertainties in the rate estimates are quantified using a Monte Carlo technique. We tested our approach using synthetic concentration profiles generated from analytical solutions to the continuity equation. We then applied the approach to concentration profiles of dissolved sulfate, sulfide, methane, and manganese in the 420-m thick sediment column of eastern equatorial Pacific Ocean Drilling Program Site 1226. Our results indicate that (i) sulfate reduction and iron reduction occur at most sediment depths, (ii) net methane production occurs in discrete depth intervals and (iii) manganese reduction occurs near the seafloor and deep in the sediments. These results provide quantitative evidence that multiple respiration pathways co-exist in the same depth intervals of these deep subseafloor sediments. Ó 2008 Elsevier Ltd. All rights reserved.
1. INTRODUCTION This study has two principal objectives. The first objective is to provide a rigorous technique for quantifying net rates of biogeochemical reactions throughout deep subseafloor sediment columns. The second objective is to quantify the extent to which multiple energy-yielding microbial reactions co-occur in deeply buried marine sediments. Previous studies have used concentration profiles of dissolved chemicals to infer, but not quantify, the co-occurrence of multi-
*
Corresponding author. E-mail address: [email protected] (G. Wang). 1 Present address: Department of Biology, Johns Hopkins University, 3400 N. Charles, Baltimore, MD 21218, USA. 0016-7037/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.gca.2008.04.024
ple metabolic reactions in deep subseafloor sediments, including sulfate reduction and methane production (Mitterer et al., 2001; Bralower et al., 2002; D’Hondt et al., 2002) and sulfate reduction, metal reduction and methane production (D’Hondt et al., 2004). Profiles of dissolved chemical concentrations in sedimentary pore fluids supply fundamental information about sedimentary processes. They are used to identify and quantify net rates of microbial reactions (e.g., Berner, 1964; Froelich et al., 1979; Devol and Ahmed, 1981; Schulz, 2000; D’Hondt et al., 2002, 2004; Burdige, 2006; Wortmann, 2006), net rates of diagenetic reactions, such as CaCO3 precipitation/dissolution caused by metabolic CO2 production (e.g., Emerson and Bender, 1981; Hales et al., 1994), and net rates of fluid advection through sediments and oceanic crust (Baker et al., 1991). This type of data is
3480
G. Wang et al. / Geochimica et Cosmochimica Acta 72 (2008) 3479–3488
particularly important for quantifying metabolic rates in deeply buried subseafloor ecosystems, where microbial activities are otherwise difficult to quantify (e.g., D’Hondt et al., 2004; Parkes et al., 2005). Two principal approaches have been used to estimate net reaction rates from dissolved chemical concentration data. The first approach uses these data to constrain rate parameters in an analytical solution of the continuity equation where reaction rates vary as a simple function of depth (e.g., Berner, 1964; Westrich and Berner, 1984). If reaction rates are not a simple function of depth, the usefulness of this approach is obviated. In deep subseafloor sediments, dissolved chemical profiles provide qualitative evidence that at least some reaction rates are not simple functions of depth. For example, profiles of dissolved iron and dissolved manganese concentrations indicate that net production of these dissolved chemicals occurs in discrete zones deep beneath the seafloor at some Ocean Drilling Program Sites (D’Hondt et al., 2004). The second approach uses finite-difference methods to directly calculate net reaction rates at discrete depths (e.g., Berg et al., 1998). Several studies have used a finitedifference approach to quantify reaction rates in deep subseafloor sediments (Bottcher et al., 2004; Schippers et al., 2005; Arndt et al., 2006; Wortmann, 2006). In this study, we use a finite-difference method. Our approach differs from most previous finite-difference studies of subseafloor reaction rates by explicitly allowing reaction zones of variable thickness and variable sample spacing and by accommodating vertical variations in porosity, tortuosity, pore water burial rate, and diffusion coefficient. It differs from all previous studies by quantifying uncertainties of the calculated rates, using a Monte Carlo technique. We provide a Matlab program for this approach, which users can download and, if necessary, adjust for their purposes. Briefly, our approach finds the maximum number of constant reaction rate depth intervals that is statistically significant. A Monte Carlo method (Bevington and Robinson, 2003) is applied to estimate the uncertainty of the calculated rates due to error associated with sampling and chemical measurements. In our study, the computer program is written in Matlab, but the principles are general. It consists of two sub-programs. Sub-program NRR1, provided in Electronic annex EA-1, estimates a reaction rate profile from measured concentrations and physical properties. Sub-program NRR2, provided in Electronic annex EA-2, quantifies uncertainties of the rates. The manual for the program is provided in Electronic annex EA-3. In this paper, we describe the approach, test its reliability with synthetic data, and apply it to pore fluid profiles from eastern equatorial Pacific Ocean Drilling Program (ODP) Site 1226. 2. PROCEDURES 2.1. Continuity equation For steady state conditions in one dimension, the mass balance of a pore-water chemical species is expressed as
( ) o DðxÞ o½/ðxÞCðxÞ 2 þ ½/ðxÞbðxÞ þ /ðxÞvðxÞCðxÞ ox ox h ðxÞ þ RðxÞ ¼ 0
ð1Þ
where C(x) is concentration in pore water, x is depth below the seafloor, /(x) is porosity, which is the ratio of the volume of pore water to the volume of the whole sediment (e.g., Berner, 1980), D(x) is the molecular diffusion coefficient, h2(x) is tortuosity, b(x) is pore water burial velocity, v(x) is externally driven flow velocity, R(x) is the rate of diagenetic reaction per unit volume of sediment, and (x) signifies that the parameter is a function of sediment depth. Positive values of R(x) refer to net production, while negative values indicate net consumption. To simplify the notation, in later equations, (x) will be omitted in most cases. Based on scaling analysis, the steady state approximation is reasonable for sediments where D/(Ls0) > 1 and R(x) does not vary on time scales shorter than L2/(2D), where L is the length scale of the sediment column and s0 is the sedimentation rate near the seafloor. 2.2. Diffusion coefficients Diffusion coefficients are calculated for in situ temperatures, based on the Stokes–Einstein relation (Bockris and Reddy, 1970) DðxÞ ¼
T ðxÞ gðT r Þ DðT r Þ g½T ðxÞ T r
ð2Þ
where D(Tr) is a solute’s diffusion coefficient at a reference temperature Tr (kelvin), T(x) is the temperature (kelvin) at depth x, g[T(x)] is the viscosity of the pore water at T(x), and g(Tr) is the viscosity of the pore water at Tr. A compilation of diffusion coefficients of dissolved chemicals at specific temperatures is in Schulz (2000). The viscosity of pure water at in situ temperature is calculated according to the empirical relationship in Lide (2004). We assume that the viscosity of seawater is the same as that of pure water. 2.3. Tortuosity and formation factor Tortuosity, h2, is calculated by the relation, h2 ¼ /f
ð3Þ
where f is the formation factor, which is the ratio of the electrical resistivity in sediment to the resistivity in pore water extracted from the sediment (e.g., McDuff and Ellis, 1979). Following the substitution of Eqs. (3) into (1), we have o D oð/CÞ ð/b þ /vÞC þ R ¼ 0 ð4Þ ox /f ox The parameters in this equation can be in any units as long as they are consistent. In our study, the following units are used: for D, m2 yr1; for C, mol m3; for x, meters below seafloor (mbsf); for b, m yr1; for v, m yr1; and for R, mol m3 yr1. When formation factor is not measured, it can be approximated by an empirical relationship with porosity. A commonly used relationship is f = a/n, where a and n are
Quantification of co-occurring reactions
constants. Several studies have derived the two constants for different sediment compositions from measurements of f and / (e.g., Berner, 1980; Boudreau, 1997; Breitzke, 2000). 2.4. Pore water burial rate and externally driven flow If there is negligible compaction, which is often a good approximation for shallowly buried sediments, the pore water burial rate is the same as the sediment burial rate. If, however, compaction is non-negligible, interstitial waters are displaced and pore water and sediment are buried at different rates (Berner, 1980). If steady state compaction and constant sedimentation rate are assumed, then the product, /(x)b(x), is constant (Berner, 1980; Burdige, 2006) and bðxÞ ¼
where 2 Df ð/j bj þ /0 v0 ÞDx jþ1 j AAj ¼ 2 D D ðDxÞ þ 2 f ðDxÞ2 þ ð/jþ1 bjþ1 /j bj ÞðDxÞ3 2 f jþ1 j D D 2 f 2 f þ ð/jþ1 bjþ1 þ /0 v0 ÞDx j jþ1 BBj ¼ 2 Df ðDxÞ2 þ 2 Df ðDxÞ2 þ ð/jþ1 bjþ1 /j bj ÞðDxÞ3 jþ1 j 2 Df 2 Df þ ð/j1 bj1 þ /0 v0 ÞDx j j1 þ 2 D D ðDxÞ2 þ ð/j bj /j1 bj1 ÞðDxÞ3 2 f ðDxÞ þ 2 f 2
D f
j
j1
and
/L ð1 /0 Þ s0 /ðxÞð1 /L Þ
ð5Þ
where /L is the porosity at a depth where compaction becomes negligible, and /0 is the porosity near the seafloor. Multiplying both sides of Eq. (5) by /(x), we have /ðxÞbðxÞ ¼ /L
3481
ð1 /0 Þ s0 ð1 /L Þ
ð6Þ
The right hand side of Eq. (6) is a constant, which indicates that the product of porosity and pore water burial rate is depth-independent although both can vary with depth. Additionally, assuming that pore water behaves as an incompressible fluid, the product, /(x)v(x), is also constant whether or not there is compaction. Thus vðxÞ/ðxÞ ¼ v0 /0
ð7Þ
where v0 is the flow rate near the seafloor. 2.5. Numerical determination of reaction rates Reaction rates are estimated numerically by the following steps: (a) data filtering; (b) numerical approximation of the reaction rate profile; (c) optimization of the number of reaction zones; and (d) estimation of uncertainty. (a) Data filtering: To minimize the effect of outliers, C, f, and / are first smoothed using a Gaussian 5-point filter. (b) Numerical approximation of the reaction rate profile: To solve for R(x), the sediment column is divided into n layers, each with the same length, the boundaries of each corresponding to discrete depth points (x0, x1, . . ., xn). The thickness of each layer is Dx. C, D, f, /, and b are interpolated to these depths. At depth xj (j = 1, . . ., n 1), Eq. (4) is approximated using a centered finite-difference method as h i h i D oð/CÞ D oð/CÞ ð/b þ /vÞC ð/b þ /vÞC /f /f ox ox jþ1=2
j1=2
Dx
j
ð8Þ Following finite-difference approximations of the two terms of the numerator [see Wang (2006) for details] and substituting /(x)v(x) with /0v0 based on Eq. (7), Eq. (8) is converted into a set of equations ð9Þ
j1
Therefore, for the sediment column there are (n 2) equations of C vs. R, which can be expressed in matrix notation: AC ¼ S where
BB1 AA1 AA2 CC 2 BB2 A ¼ CC n2 BBn2 CC n1 R1 CC 1 C 0 R2 S ¼ Rn2 Rn1 AAn1 C n
ð10Þ ; AAn2 BBn1
and
To solve for R, a least-square minimization method is used. R(x) is assumed to be a piecewise constant function. The sediment column is divided into N zones with constant reaction rates. The length of each zone is determined by a specified number of chemical measurements. In other words, the reaction rate profile over the depth range 0 to L is described by N constant rates. From this rate profile, we generate [R1, . . ., Rn1] corresponding to depths [x1, . . ., xn1]. The smoothed measured concentrations (Cm) at the top and bottom of the sediment column are taken as the boundary conditions so that C0 and Cn are known. The concentration estimated by [R1, . . ., Rn1] is ^ ¼ A1 S C
¼ Rj
AAj C jþ1 þ BBj C j þ CC j C j1 ¼ Rj
2 Df þ ð/j bj þ /0 v0 ÞDx j1 j CC j ¼ 2 D D ðDxÞ2 þ ð/j bj /j1 bj1 ÞðDxÞ3 2 f ðDxÞ þ 2 f 2 Df
where C1 C2 ^ ¼ C C n2 C n1
ð11Þ
3482
G. Wang et al. / Geochimica et Cosmochimica Acta 72 (2008) 3479–3488
We used the Thomas algorithm (Von Rosenberg, 1969) ^ in our Matlab program as A is a tri-diagonal to solve for C ^ is then interpolated to the depths of the measurematrix. C ments and the sum of squared errors (SSE) between Cm and ^ is calculated. SSE is then minimized to determine C [R1, . . ., Rn1] (and hence a reaction rate profile with values of N rates) that provides a best-fit concentration profile to the measured concentration profile. We used Fminsearch, a multi-dimensional minimization function in Matlab that uses the simplex search method of Lagarias et al. (1998), for this purpose, though other minimization techniques could be used. The minimization procedure we wrote in Matlab is provided in Electronic annex EA-4. For a total number of concentration data points, M, each at a different depth, the maximum number of reaction zones is equal to (M 2) because at least three points are needed to approximate a reaction rate. Using (M 2), however, most likely will result in relatively noisy estimated rates, of which many are not statistically significant, as is the case when reaction rates are directly calculated using a finite-difference method at discrete depths. It also results in long runtimes of the program. From our tests with synthetic concentration profiles, M/3 is a reasonable compromise, for the maximum number of zones, between runtime and potential number of reaction zones that are statistically significant. That is, at least three concentration data points, at different depths, must be available within each reaction zone to determine a constant rate. In this portion of the model, N is initially set at 1 and then increased to the maximum number of reaction zones, e.g., M/3. The minimization is repeated for each N, and subsequently, a set of rate profiles and corresponding concentration profiles are generated. (c) Optimization of the number of reaction zones: Two steps are used to optimize the reaction rate profile. First the largest number of statistically significant reaction zones is determined from the set of rate profiles generated above, using partial F-tests. Second, adjacent zones with similar rates are combined into zones with constant rates if combining does not decrease the statistical significance. For the first step, the null hypothesis tested here is, given a number of reaction zones, K (i.e., specific number of variables), the addition of more zones (i.e., more variables) does not significantly improve the fit of modeled concentrations to the measured concentrations. The F value for the hypothesis is calculated as (Kleinbaum and Kupper, 1978): F ¼
SSEK SSEL LK SSEL ML1
ð12Þ
where L is a number of reaction zones (L > K). A probability is calculated for each F value and compared with the significance level the user chooses (0.05 is commonly used) to accept or reject the hypothesis. In the second step, adjacent zones with the most similar rates are combined and then the minimization is repeated. The merging of adjacent zones and re-determination of reaction rates are repeated until there is only one zone for the entire data set. The new set of fit profiles is compared again using partial F-tests to determine a final rate profile.
(d) Uncertainty estimates: A Monte Carlo method is used because the uncertainty limits in the estimates of R due to analytical and sampling uncertainty cannot be expressed in a simple formula. The Monte Carlo method we use in NRR2 makes use of randomly determined numbers to simulate processes that involve an element of chance, as there is in the errors associated with analysis and sampling. In our application, for each sample, random chemical measurements (Cr) are assumed to be described by a Gaussian distribution, with a mean of Cm and a relative deviation due to uncertainties associated with chemical measurements and sampling (e). A random measurement can then be expressed as: C r ¼ C m ð1 þ e rÞ
ð13Þ
where r is a random number generated from the standard Gaussian distribution N(0, 1). In our approach, a set of randomized concentration profiles is generated, based on the smoothed original (measured) concentration profile, Cm. A rate profile is estimated for each randomized concentration profile, using the same procedures and the same values of other parameters as in NRR1. A standard deviation is calculated at each depth of the estimated rates from this set of rate profiles, taking the rates calculated from the original concentration profile as the mean. Using a Monte Carlo method to determine the uncertainty of estimated rates, the greater the number of random concentration profiles generated, the closer the calculated uncertainty is to the true uncertainty. However, the more profiles generated, the longer it takes to run NRR2. From our tests, generation of 50 random profiles provides a good compromise between length of runtime and uncertainties. 3. ASSESSMENT 3.1. Reliability test The approach was tested following three steps: (i) creating a concentration profile with an analytical solution and a specified rate profile (hereafter referred to as analytical rates), (ii) determination of rates from the analytical concentration profile using the program, and (iii) comparison of the model-derived rates with the analytical rates. We tested the approach with constant analytical rates, variable analytical rates, and varied sample spacing. In all of these tests (Table 1), the approach reproduced the analytical rates very well (Fig. 1). As an illustrative example, we present Test 1, which is a test with a relatively complicated reaction rate profile where there is an increasing rate of depletion of a dissolved species with depth and then an abrupt decrease in rate to zero at 50 m. This test works well in demonstrating the strength of our approach in identifying both smoothly varying and abrupt changes of rates. We parameterized the consumption rate, under the assumption that it increases exponentially with depth, RðxÞ ¼ Rð0ÞeðaxÞ
ð14Þ
until the concentration reaches 0 at 50 mbsf below which the rate is zero, where R(0) is the consumption rate at
Quantification of co-occurring reactions
3483
Table 1 Tests of the numerical method Parameters R (mol m3 yr1) 2
1
D (m yr ) f / v (m yr1) b (m yr1) BCa: C (mol m3), x (m)
Test 2
Test 3
Test 4
1.5 104
1.5 104
0
0.018 2.2 0.8 1.0 105 2.5 105 C(x = 1) = 30 C(x = 50) = 0
0.018 2.2 0.8 0 0 C(x = 1) = 30 C(x = 100) = 25
0.018 2.2 0.8 1.0 105 2.5 105 C(x = 1) = 30 C(x = 100) = 25
0.018 2.2e0.002x 0.8 0 0 C(x = 1) = 30 C(x = 200) = 25
BC is the abbreviation for boundary conditions.
a
0
Depth (mbsf)
a
Test 1 1:5 104 e0:01x ; 0 6 x 6 50 x > 50 0;
20
b
analytical concentrations modeled rates
40
analytical rates
60
1 σ uncertainty envelope
80 100 0
Depth (mbsf)
c
10 20 30 -3 -2 -1 0 Concentration Reaction rate (mol m -3 ) (10-4 mol m-3 yr -1)
-2 -1 10 20 30 0 Reaction rate Concentration (10-4 mol m-3 yr -1) (mol m-3 )
0
d
0 20
0
50
40 100 60 150
80
200
100 0
10 20 30 0 -2 -1 Concentration Reaction rate (mol m-3 ) (10 -5 mol m-3 yr -1)
0
10 20 30 -2 0 2 Concentration Reaction rate (mol m-3 ) (10-5 mol m-3 yr -1)
Fig. 1. Method assessment results. Rates and uncertainties estimated by the numerical method are compared to the specified analytical rate profiles, which are described in detail in Table 1. (a) Test 1: variable reaction rates with depth, (b) Test 2: constant reaction rate with depth and negligible pore water burial rate and externally driven flow rate (i.e., b, v = 0), (c) Test 3: constant reaction rate with depth, and pore water burial rate and externally driven flow not equal to 0, and (4) Test 4: no reaction (i.e., R = 0), b and v equal to 0, and f changing with depth.
x = 0, and a is a constant. Then Eq. (4) was solved analytically assuming D, f, and / are constant with depth /f ðvþbÞ x D
C ¼P þQe
D
a2
Rð0Þ f eax / f ðv þ bÞ a
ð15Þ
where P and Q are constants. From this analytical solution, we generated a concentration profile for a 100-m thick sediment column. The sediment depth evaluated started at 1 mbsf. The concentration at the top of the sediment column was set to 30 mol m3, D = 1.80 102 m2 yr1, f = 2.2, / = 0.8, v = 1.0 105 m yr1, 5 1 b = 2.5 10 m yr , R(0) = 1.50 104 mol m3 yr1, and a = 0.01 m1. P and Q were solved from the boundary
conditions in Table 1. A variable depth increment was used to generate a concentration profile that has sample spacing similar to real measurements. Along with this analytical concentration profile, constant D, f, /, v, and b were input into NRR1. The chosen minimum number of concentration data points in each reaction zone was 3, since the analytical concentration profile is smooth (Fig. 1a), and a significance level of 0.05 was used. NRR1 was run with a maximum of 12 zones and identified 7 statistically significant reaction zones. In the subsequent run of NRR2, 50 random concentration profiles were generated, in which a relative standard deviation in concentration of 1% was assumed.
3484
G. Wang et al. / Geochimica et Cosmochimica Acta 72 (2008) 3479–3488
The results, compared with the analytical rates, are shown in Fig. 1a. The estimated rates clearly identify and reproduce the exponential increase in magnitude in the analytical rate. Moreover, the estimated rates explicitly reflect the abrupt change of the analytical rates near 50 mbsf. Additionally, the 1r uncertainty envelope closely brackets the analytical rates. In each reaction zone, the estimated constant rate is approximately equal to the mean of the analytical rates. The results of this test indicate that our approach performs well in identifying variable and abruptly changing rates from a concentration profile and in providing an accurate estimate of rates with objective uncertainty limits. 3.2. Application of the approach In most deeply buried sediments, metabolic activities are too low to be quantified with radiotracer methods. Additionally, low spatial resolution in sampling has previously limited the applicability of pore fluid modeling in deep sediments. ODP Site 1226 provides high-resolution sampling of a deep sediment column, which makes it ideal for the application of our method. At this site, co-occurrence of multiple metabolic activities is suggested from qualitative analysis of chemical profiles (D’Hondt et al., 2004). The application of the method described in this paper to data from Site 1226 illustrates the power of the method for quantifying metabolic activity. These results provide quantitative evidence of the co-existence of multiple metabolic activities through a deep sediment column and can be used to test basic aspects of metabolism in this environment. In particular, this work quantifies rates of sulfate, iron, and manganese reduction along with methane generation as a function of depth. Site 1226 is located in the eastern equatorial Pacific at a water depth of 3297 m (Shipboard Scientific Party, 2003). The sediment column at this site is 420-m thick. Seawater flows through the underlying basalt driving sulfate back up toward seawater concentrations at the bottom of the hole (D’Hondt et al., 2003). High-resolution measurements are available for dissolved sulfate, total dissolved sulfide, dissolved methane, dissolved manganese, porosity, and formation factor throughout the sediment column (Fig. 2) (Shipboard Scientific Party, 2003). In situ down-hole temperature (t, °C) is given by: t = 1.74 + 0.0572x (Shipboard Scientific Party, 2003). The sedimentation rate at Site 846 (Shipboard Scientific Party, 1992), which is within 100 m of Site 1226, was taken as the sedimentation rate at Site 1226. The effect of compaction on the pore water burial rate was considered in our calculation. There is no identifiable advection at this site (Shipboard Scientific Party, 2003). A significance level of 0.05 was used in the calculation of the significance of reaction rates. The combined analytical and sampling uncertainty for each chemical species was estimated from the deviations of the measured profiles from their smoothed profiles. These are 1%, 10%, 7%, and 1% for sulfate, sulfide, methane, and manganese, respectively. For each chemical species, 50 random concentration profiles were generated by NRR2 to estimate uncertainty by the Monte Carlo method.
The number of statistically significant reaction intervals varied between chemical species. For dissolved sulfate (Fig. 2c) and methane (Fig. 2d), there were 5 intervals. For dissolved sulfide (Fig. 2e) and manganese (Fig. 2f), there were 6 and 14 intervals, respectively. The confidence limits were 6%, 22%, 15%, and 75% for sulfate, sulfide, methane, and manganese, respectively. The concentrations predicted from the calculated rates fit the observed data very well. Overall, sulfate is the dominant terminal electron acceptor in the sediment column although in some intervals, its reduction rate is below the limit of detection. The depthintegrated net sulfate reduction rate from 3 to 418 mbsf is 1.8 103 ± 1 104 mol m2 yr1. The depth-integrated net rate of dissolved sulfide production from 0.59 to 410 mbsf is 7.6 105 ± 1.7 105 mol m2 yr1. The depth-integrated net production rate of dissolved manganese from 0.59 to 418 mbsf is 1.2 105 ± 9 106 mol m2 yr1. This negative value is almost certainly a consequence of beginning the Mn reaction rate modeling at 0.59 mbsf; Mn produced at shallower depths migrates downward past this depth to be precipitated below (Fig. 2f). Unfortunately, there are no dissolved Mn data for depths between 0 and 0.59 mbsf (Shipboard Scientific Party, 2003). The depth-integrated net production rate of dissolved methane from 11.9 to 384 mbsf is 3.4 107 ± 5 108 mol m2 yr1. Sulfate reduction rates are greatest in the uppermost 40 m, but the rates do not decrease smoothly with depth as is often assumed. Sulfate reduction occurs in the upper sediment column (above 40 mbsf) at a rate of 2.5 105 mol m3 yr1, in the intervals between 55 and 131 mbsf at a rate of 4.5 106 mol m3 yr1, between 198 and 251 mbsf at a rate of about 5.7 106 mol m3 yr1, and from 275 mbsf to the bottom of the column at a rate of 6.6 107 mol m3 yr1. For the other reaction zones (from 40 to 55 mbsf, from 131 to 198 mbsf, and from 251 to 275 mbsf), the rates are indistinguishable from 0. The depths of two of the sulfate reduction maxima (one near the seafloor around 40 mbsf and the other around 250 mbsf) correspond to the depths of maximum rates determined by radiotracer measurements at this site and approximately correspond to the chemical interfaces that have been suggested to stimulate microbial activities (Parkes et al., 2005). Sulfide, methane, and manganese profiles represent the net difference between production and consumption. Sulfide is produced via sulfate reduction and consumed by precipitation into sulfide minerals. There is net consumption near the seafloor (from 0.21 to 9 mbsf) at a rate of 3.4 105 mol m3 yr1 and from 256 to 286 mbsf at a rate of 2.1 ;106 mol m3 yr1. Net production occurs in the intervals from 9 to 36 mbsf at a rate of 1.3 105 mol m3 yr1, and from 36 to 178 mbsf at a rate about one order of magnitude smaller than the rate in the overlying interval (Fig. 2e). This interval from 36 to 178 mbsf covers some of the depths (between 40 and 55 mbsf and between 131 and 178 mbsf) at which rates of sulfate reduction are too small to be quantified from the sulfate profile, indicating that sulfate reduction still occurs
Quantification of co-occurring reactions
a
3485
c
b
0 Depth (mbsf)
input measurement data 100
modeled rates and concentrations
200
1σ uncertainty envelope
300 400 0.6 0.7 0.8 0.9 1.5 2.5 3.5 4.5 18 Porosity Formation factor
Depth (mbsf)
0
22 26 30 Sulfate (mol m-3)
-2 -1 0 Reaction rate -5 (10 mol m-3 yr -1) -3
e
d
100 150
200
250
300
350 400 0
1 2 Methane -3 (10 mol m-3)
-1
0 1 Reaction rate -8 (10 molm-3 yr -1)
0
0.2 0.4 0.6 Sulfide (mol m-3)
-3
-2 0 2 Reaction rate (10-5 mol m-3 yr -1)
450
0 2 -2 0 0.2 0.4 0.6 Reaction rate Sulfide -6 (mol m-3) (10 mol m-3 yr -1)
f 0
Depth (mbsf)
0
10
100
20 30
200
40
0 2 4 -8 -4 Manganese Reaction rate (10-2 mol m-3) (10-6 mol m-3 yr -1 )
0
300 400 0
2 4 -8 -4 0 Manganese Reaction rate (10-2 mol m-3) (10-6 mol m-3 yr -1)
Fig. 2. Calculated reaction rates based on concentrations and physical properties at ODP Site 1226. (a) Porosity, (b) formation factor, (c) dissolved sulfate, (d) dissolved methane, (e) total dissolved sulfide, and (f) dissolved manganese. Physical properties and dissolved chemical concentrations are from Shipboard Scientific Party (2003).
in these intervals at rates at least the same as the net rate of sulfide production. From 178 to 256 mbsf, net reaction rate of sulfide is indistinguishable from 0. From 286 mbsf to the bottom of the column, dissolved sulfide is below the detection limit (0.17 lM) and we take sulfide concentrations in this interval to be 0 (Fig. 2e). Modeled net sulfide reaction rates are close to 0 in this interval (Fig. 2e). The distributions of the net rates of dissolved sulfate consumption and dissolved sulfide production and consumption suggest that sulfate reduction occurs at most, if not all, of the sediment depths (rates are indistinguishable from 0 in the interval from 178 to 198 mbsf). Dissolved sulfide produced by sulfate reduction is precipitated almost completely throughout the sediment column except in the intervals from 40 to 55 mbsf and from 131 to 178 mbsf, and from 251 to 256 mbsf, where there is no quantitative evidence of net consumption of dissolved sulfide. Based
on the depth-integrated rates, about 96% of the dissolved sulfide produced by sulfate reduction is precipitated in the sediment column. The most likely ultimate sink for dissolved hydrogen sulfide is pyrite (FeS2), which is present throughout the sediments at this site (Shipboard Scientific Party, 2003) and is common in anoxic sediments. This local precipitation of iron sulfide implies that reduced iron [Fe(II)] must be produced in depth intervals where sulfate reduction occurs because there are almost no gradients of dissolved Fe(II) concentrations throughout the sediment column in close proximity to where sulfate reduction occurs (Fig. 3). Considering that concentrations of dissolved Fe(II) remain low in the sediment column, the rates of iron reduction and sulfate reduction seem to be balanced locally. For example, assuming pyrite is the major sink of dissolved Fe(II) and dissolved sulfide, the rate of iron reduction
3486
G. Wang et al. / Geochimica et Cosmochimica Acta 72 (2008) 3479–3488 0
Depth (mbsf)
100
200
300
400
0
0.02 0.04 Iron (mol m-3)
Fig. 3. The concentration profile of dissolved iron in the sediments of Site 1226 (Shipboard Scientific Party, 2003).
would be about half of the rate of sulfate reduction. In this case, the depth-integrated rate of iron reduction is 8.6 104 ± 5 105 mol m2 yr1. In anoxic marine sediments, methane is produced by microbial methanogenesis and consumed by microbially mediated anaerobic methane oxidation using sulfate as the terminal electron acceptor (Devol and Ahmed, 1981) as well as by adsorption onto sediments (Whiticar and Suess, 1990; Knies et al., 2004). At Site 1226, net methane production occurs in the intervals between 50 and 107 mbsf at a rate of 1.4 108 mol m3 yr1 and from 164 to 279 mbsf at a smaller rate of 2.2 109 mol m3 yr1 (Fig. 2d). Methane is consumed and/or adsorbed in the intervals from 3.2 to 50 mbsf at a rate of 4.1 109 mol m3 yr1, from 107 to 164 mbsf at a rate of 6.4 109 mol m3 yr1, and from 280 to 384 mbsf at a rate of 1.6 109 mol m3 yr1. Net dissolved manganese production occurs near the seafloor: from 4.5 to 15 mbsf at a rate of 9.0 107 mol m3 yr1 and from 40 to 78 mbsf at a much smaller rate of 5.7 108 mol m3 yr1 (Fig. 2f). It also occurs deep in the sediment column between 272 and 313 mbsf at a rate of about 3.0 107 mol m3 yr1 and from 391 mbsf to the bottom of the column at a rate of 1.2 107 mol m3 yr1. At other depths, dissolved manganese is consumed (e.g., due to precipitation into MnCO3 or in solid solution with CaCO3) or the rate is indistinguishable from 0. Usually it is inferred that net manganese reduction occurs at a rate at least the same as the rate of net dissolved manganese production (e.g., D’Hondt et al., 2004). In theory, some dissolved manganese may be released from solid phases when they become undersaturated. However, there is no indication from other solutes such as Ca2+ and alkalinity (Shipboard Scientific Party, 2003) that the local production of dissolved Mn is from dissolution. Thus, it is most likely due to manganese reduction. Except in the intervals from 178 to 198 mbsf and from 251 to 256 mbsf, zones of net dissolved methane production correspond to zones where sulfate reduction clearly occurs,
indicating that sulfate reduction does not result in the competitive exclusion of methanogenesis in these intervals. Net dissolved manganese production occurs in intervals near the seafloor and deep in the sediment column, where sulfate reduction clearly occurs. Our results quantitatively confirm the co-occurrence of sulfate reduction, net methanogenesis, and net manganese production at this site and indicate that iron reduction occurs in the depth intervals where sulfate reduction occurs. Besides microbial reduction of iron and manganese, abiotic reduction of iron and manganese by sulfide may co-occur in these sediments as observed in sediment incubations (Canfield, 1989; Koretsky et al., 2003). Our results are in agreement with the qualitative interpretation of D’Hondt et al. (2004). To test the effect of downcore variations of porosity and formation factor on the modeled reaction rates, we calculated sulfate reaction rates at Site 1226 using the downcore average formation factor and porosity. Rates calculated with constant formation factor and porosity do not resolve spatial variations in rates that are identified when measured formation factors and porosities are used. The estimated rate profile is consisted of 2 reaction zones compared with the 5 reaction zones estimated using measured formation factor and porosity. In addition, to examine the effect of compaction on calculated rates, we estimated sulfate reaction rate at this site using measured formation factor and porosity without considering compaction. The resulting rate profile has 6 reaction zones with rates about 16% different from the rates estimated considering compaction. More importantly, the uncertainties of the rates are about twice those when compaction is considered. Results of the application of the approach to these data demonstrate that our approach is an effective tool for quantifying metabolic activities in deeply buried sediments from profiles of dissolved chemicals and sediment physical properties. Our approach is also a useful tool in studying biogeochemical cycles in deep sediments. For example, our Site 1226 results indicate that removal of dissolved sulfide from the pore water into solid phases occurs at most sediment depths and this sediment column is a significant sink of sulfur. In addition, our approach provides a realistic description of metabolic reactions in deeply buried sediments where direct radiotracer measurements are too insensitive or problematic due to difficult sample manipulations (Jørgensen et al., 2001). 4. COMMENTS AND RECOMMENDATIONS ON THE METHOD Our method estimates reaction rates without assuming any a priori relationship of rate to sediment depth and incorporates variations in physical properties that affect pore-water chemical gradients and inferred reaction rates. It also successfully deals with variable spacing of data and changes in sediment physical properties. In comparison, the use of equal length for each zone (e.g., Berg et al., 1998) can overestimate rates in intervals with few concentration data or miss reaction zones in relatively narrow intervals with apparent chemical gradients. Gradients of chemical profiles caused by changes of physical proper-
Quantification of co-occurring reactions
ties are considered by the incorporation of in situ distributions of physical properties. Using the numerical procedures, the approach provides objective judgments on where and how much a dissolved chemical is produced or consumed. The uncertainties in estimated reaction rates that are introduced due to the analytical error associated with chemical measurements are well constrained using a Monte Carlo simulation of concentrations. Although our approach effectively quantifies net rates of respiration reactions in deeply buried sediments, it has limitations. The first limitation is the assumption of constant rates in each reaction zone. Rates within a zone may vary on scales too fine to resolve; the constant rate calculated for a reaction zone is the approximate average of the rates within the zone. The second limitation is, as for all numerical approximations, that the rates at the boundary of a column cannot be predicted without additional inputs. For example, rates cannot be calculated at the seafloor and the sediment/basement interface. Third, the user determines the minimum number of concentration data points in each reaction zone, which is somewhat subjective. This determination may require a trial-and-error process because the user will not want to overestimate rates by using too small a number, nor will the user want to miss apparent reaction zones by using too large a number. Our approach accounts for variations in temperature, diffusivity, formation factor, porosity, viscosity, externally driven flow velocity, and pore water burial rate with depth to reflect the complexity of deep sedimentary environments. However, users have the option of using constant values for these parameters to fit their study environments. If a user wants more smoothing of his/her measurement data than 5-point Gaussian smoothing, locally weighted regression scatter plot smoothing (Lowess) (Chambers et al., 1983) is recommended. Lowess smoothing can be found in a Matlab toolbox and other commercial software. The user can smooth the data profiles using Lowess before loading them into the program, or embed Lowess in NRR1 if the toolbox is available. For very noisy chemical data, the uncertainty of the depth-integrated rate may be 100% of the rate or even higher. In this case, the measured data may not define profiles clean enough for the approach to identify reaction zones. The mass-balance equation that we used makes this approach applicable in environments where advection may be one of the dominant transport processes. Mathematically, externally driven flow is similar to the pore water burial process as can be seen in Eq. (4) and the model handles the external flow rate almost the same way as the pore water burial rate. For radioactive chemicals, the user can incorporate radioactive decay into Eq. (4) and follow the same procedures to estimate reaction rates. The user must keep in mind that the approach applies to dissolved species approaching steady state. 5. CONCLUSIONS Our approach successfully identifies variable reaction rates from sedimentary pore fluid concentration profiles and provides a good approximation of rates with con-
3487
strained uncertainties. It is an effective tool for using profiles of dissolved chemicals and physical properties to study metabolic activities in deeply buried sediments. Application of our approach to data from ODP Site 1226 provides quantitative evidence of the depth distributions of microbial processes deep beneath the seafloor. These results quantitatively demonstrate that sulfate reduction and methanogenesis occur in the same reaction zones deep in a subseafloor sediment column, as do sulfate reduction and dissolved manganese production. This sediment column is a sink for sulfur, where the dissolved sulfide produced from sulfate reduction is precipitated almost completely into solid phases in the depth intervals where sulfate reduction occurs. The existence of this sulfide sink suggests that iron reduction may also occur throughout the sediment column in the reaction zones where sulfate reduction occurs. ACKNOWLEDGMENTS This manuscript was greatly improved by the thoughtful comments of Andrew Dale, Timothy Ferdelman, David Burdige, and an anonymous reviewer. We thank Andrew Dale and Heather Schrum for independently testing the Matlab programs and Heather Schrum for reviewing the manual for the Matlab programs. This study was supported by the NASA Astrobiology Institute.
APPENDIX A. SUPPLEMENTARY DATA Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.gca.2008. 04.024. REFERENCES Arndt S., Brumsack H. and Wirtz K. W. (2006) Cretaceous black shales as active bioreactors: a biogeochemical model for the deep biosphere encountered during ODP Leg 207 (Demerara Rise). Geochim. Cosmochim. Acta. 70, 408–425. Baker P. A., Stout P. M., Kastner M. and Elderfield H. (1991) Largescale lateral advection of seawater through oceanic crust in the central equatorial Pacific. Earth Planet. Sci. Lett. 105, 522–533. Berg P., Risgaard-Petersen N. and Rysgaard S. (1998) Interpretation of measured concentration profiles in sediment pore water. Limnol. Oceanogr. 43, 1500–1510. Berner R. A. (1964) An idealized model of dissolved sulfate distribution in recent sediments. Geochim. Cosmochim. Acta 28, 1497–1503. Berner R. A. (1980) Early Diagenesis. Princeton University Press, Princeton. Bevington P. R. and Robinson D. K. (2003) Data Reduction and Error Analysis for the Physical Sciences, third ed. McGraw Hill. Bockris J. O’M. and Reddy A. K. N. (1970). Modern Electrochemistry, vol. 1. Plenum Press, New York. Bottcher M. E., Khim B. K., Suzuki A., Gehre M., Wortmann U. G. and Brumsack H. J. (2004) Microbial sulfate reduction in deep sediments of the Southwest Pacific (ODP Leg 181, Sites 1119–1125): evidence from stable sulfur isotope fractionation and pore water modeling. Mar. Geol. 205, 249–260. Boudreau B. P. (1997) Diagenetic Models and their Implementation: Modeling Transport and Reactions in Aquatic Sediments. Springer-Verlag, Berlin.
3488
G. Wang et al. / Geochimica et Cosmochimica Acta 72 (2008) 3479–3488
Bralower T. J., Premoli Silva I., Malone M. J., et al. (2002). Proc. ODP, Init. Repts. 198 (Online). Available from: http://wwwodp.tamu.edu/publications/198_IR/198ir.htm, accessed on August 8, 2006. Breitzke M. (2000) Physical properties of marine sediments. In Marine Geochemistry (eds. H. D. Schulz and M. Zabel). Springer-Verlag, Berlin, pp. 29–83. Burdige D. J. (2006) Geochemistry of Marine Sediments. Princeton University Press, Princeton. Canfield D. E. (1989) Reactive iron in marine sediments. Geochim. Cosmochim. Acta 53, 619–632. Chambers J. M., Cleveland W. S., Kleiner B. and Tukey P. A. (1983) Graphing Methods for Data Analysis. Wadsworth International Group, Belmont, California, and Duxbury Press, Boston. D’Hondt S., Jørgensen B. B., Miller D. J., Batzke A., Blake R., Cragg B. A., Cypionka H., Dickens G. R., Ferdelman T., Hinrichs K. U., Holm N. G., Mitterer R., Spivack A., Wang G., Bekins B., Engelen B., Ford K., Gettemy G., Rutherford S. D., Sass H., Skilbeck C. G., Aiello I. W., Gue`rin G., House C., Inagaki F., Meister P., Naehr T., Niitsuma S., Parkes R. J., Schippers A., Smith D. C., Teske A., Wiegel J., Padilla C. N. and Acosta J. L. S. (2004) Distributions of microbial activities in deep subseafloor sediments. Science 306, 2216–2221. D’Hondt S., Jørgensen B. B., Miller D. J., et al. (2003) Proc. ODP Init. Repts. 201 [CD-ROM]. Available from: Ocean Drilling Program, Texas A&M University, College Station, TX 778459547, USA. D’Hondt S., Rutherford S. and Spivack A. J. (2002) Metabolic activity of subsurface life in deep-sea sediments. Science 295, 2067–2070. Devol A. H. and Ahmed S. I. (1981) Are high rates of sulfate reduction associated with anaerobic oxidation of methane? Nature 291, 407–408. Emerson S. and Bender M. (1981) Carbon fluxes at the sedimentwater interface of the deep-sea: calcium carbonate preservation. J. Mar. Res. 39, 139–162. Froelich P. N., Klinkhammer G. P., Bender M. L., Luedtke N. A., Heath G. R., Cullen D., Dauphin P., Hammond D. and Hartman B. (1979) Early oxidation of organic matter in pelagic sediments of the eastern equatorial Atlantic: suboxic diagenesis. Geochim. Cosmochim. Acta 43, 1075–1090. Hales B., Emerson S. and Archer D. (1994) Respiration and dissolution in the sediments of the western North Atlantic: estimates from models of in situ microelectrode measurements of porewater oxygen and pH. Deep-Sea Res. I I41, 695–719. Jørgensen B. B., Weber A. and Zopfi J. (2001) Sulfate reduction and anaerobic methane oxidation in Black Sea sediments. DeepSea Res. I 48, 2097–2120. Kleinbaum D. G. and Kupper L. L. (1978) Applied Regression Analysis and Other Multivariable Methods. Duxbury Press, North Scituate, Mass. Knies J., Damm E., Gutt J., Mann U. and Pinturier L. (2004) Near-surface hydrocarbon anomalies in shelf sediments off Spitsbergen: evidences for past seepages. Geochem. Geophys. Geosyst. 5. doi:10.1029/2003GC000687.
Koretsky C. M., Moore C. M., Lowe K. L., Meile C., Dichristina T. J. and Cappellen P. V. (2003) Seasonal oscillation of microbial iron and sulfate reduction in saltmarsh sediments (Sapelo Island, GA, USA). Biogeochemistry 64, 179–203. Lagarias J. C., Reeds J. A., Wright M. H. and Wright P. E. (1998) Convergence properties of the Nelder-Mead simplex method in low dimensions. SIAM J. Optim. 9, 112–147. Lide D. R. (ed.) (2004) CRC Handbook of Chemistry and Physics, 84th ed. CRC press, 2003–2004. McDuff R. E. and Ellis R. A. (1979) Determining diffusion coefficients in marine sediments: a laboratory study of the validity of resistivity techniques. Am. J. Sci. 279, 666–675. Mitterer R. M., Malone M. J., Goodfriend G. A., Swart P. K., Wortmann U. G., Logan G. A., Feary D. A. and Hine A. C. (2001) Co-Generation of hydrogen sulfide and methane in marine carbonate sediments. Geophys. Res. Lett. 28(20), 3931– 3934. Parkes R. J., Webster G., Cragg B. A., Weightman A. J., Newberry C. J., Ferdelman T. G., Kallmeyer J., Jørgensen B. B., Aiello I. W. and Fry J. C. (2005) Deep sub-seafloor prokaryotes stimulated at interfaces over geological time. Nature 436, 390– 394. Schippers A., Neretin L. N., Kallmeyer J., Ferdelman T. G., Cragg B. A., Parkes R. J. and Jørgensen B. B. (2005) Prokaryotic cells of the deep sub-seafloor biosphere identified as living bacteria. Nature 433, 861–864. Schulz H. D. (2000) Quantification of early diagenesis: dissolved constituents in marine pore water. In Marine Geochemistry (eds. H. D. Schulz and M. Zabel). Springer-Verlag, Berlin, pp. 85– 128. Shipboard Scientific Party (1992) Site 846. In Proc. ODP Init. Repts. 138 (Pt. 1) (eds. L. Mayer, N. Pisias, T. Janecek, et al.). Ocean Drilling Program, College Station, TX, pp. 265–333. Shipboard Scientific Party (2003) Site 1226. In Proc. ODP Init. Repts. 201 [CD-ROM] (eds. S. L. D’Hondt, B. B. Jørgensen, D. J. Miller, et al.). Available from: Ocean Drilling Program, Texas A&M University, College Station, TX 77845-9547, USA. Von Rosenberg D. U. (1969) Methods for the Numerical Solution of Partial Differential Equations. American Elsevier Pub. Co., New York. Wang G. (2006) Metabolic Activities in Deep Subseafloor Sediments. Ph.D. thesis. University of Rhode Island. Westrich J. T. and Berner R. A. (1984) The role of sedimentary organic matter in bacterial sulfate reduction: the G model tested. Limnol. Oceanogr. 29, 236–249. Whiticar M. J. and Suess E. (1990) Characterization of sorbed volatile hydrocarbons from the Peru margin, Leg 112, Site 679, 680/681, 682, 684, and 686/687. Proc. ODP Sci. Res. 112. doi:10.2973/odp.proc.sr.112.165.1990. Wortmann U. G. (2006) A 300 m long depth profile of metabolic activity of sulfate-reducing bacteria in the continental margin sediments of South Australia (ODP Site 1130) derived from inverse reaction-transport modeling. Geochem. Geophys. Geosyst. 7. doi:10.1029/2005GC001143. Associate editor: David J. Burdige