The influence of porosity gradients on mixing coefficients in sediments

Geochimica et Cosmochimica Acta 71 (2007) 961–973 www.elsevier.com/locate/gca

The influence of porosity gradients on mixing coefficients in sediments Filip J.R. Meysman a

a,*

, Volodymyr S. Malyuga a, Bernard P. Boudreau b, Jack J. Middelburg a

Centre for Estuarine and Marine Ecology, The Netherlands Institute of Ecology (NIOO-KNAW), Korringaweg 7, 4401 NT Yerseke, The Netherlands b Department of Oceanography, Dalhousie University, Halifax , NS, Canada NS B3H 4J1 Received 22 June 2006; accepted in revised form 1 November 2006

Abstract The surface layer of aquatic sediments is a zone characterized by both porosity gradients and intensive mixing. In the standard approach, porosity gradients are ignored when estimating mixing intensity. Here, model formulations with both constant and varying porosity are contrasted to estimate mixing coefficients Db from tracer depth profiles. Complementing the well-known exponential solution of the constant-porosity model, we present a general solution to the variable-porosity model in terms of hypergeometric functions. When using these models in a forward way, the tracer activities predicted by the variable-porosity model are higher than those generated by the constant-porosity model. Similarly, when inverse modelling, Db values estimated by the variable-porosity model are systematically higher than those derived from the constant-porosity model. Still, differences in Db values remain relatively small. When applying both mixing models to excess 210Pb data profiles from slope sediments, a maximal difference of 30% is obtained between Db values, the average deviation being 16%. A systematic exploration of parameter space predicts a maximal underestimation of 60% when deriving Db values from the constant-porosity mixing model. Given the uncertainty imposed by other model assumptions underlying the diffusive mixing model, the influence of porosity gradients on Db values must be classified as rather modest. Hence, the current mixing coefficient database is not biased by the constant porosity approximation.  2006 Elsevier Inc. All rights reserved.

1. Introduction The top layer (upper decimetres) of aquatic sediments is continuously reworked by the interplay of biological activity (feeding, moving, defecation. . .) and physical processes (erosion, resuspension, waves. . .). This mixing strongly affects the spatial distribution of various types of ‘‘particles’’, ranging from truly chemical particles (organic matter, metal oxides, adsorbed contaminants) to biological ones (viruses, bacteria, macrofaunal eggs) (Meysman et al., 2006). The intensity of this particle transport can be determined from the depth distribution of certain naturally occurring or deliberately released tracers (e.g., Li et al., 1985; Legeleux et al., 1994; Thomson et al., 2000). Typical

*

Corresponding author. Fax: +31 113 573 616. E-mail address: [email protected] (F.J.R. Meysman).

0016-7037/$ - see front matter  2006 Elsevier Inc. All rights reserved. doi:10.1016/j.gca.2006.11.024

tracers are either inert particles (e.g., luminophores, glass beads) or particle-associated chemicals that have a wellcharacterized reactivity (e.g., radio-tracers, photosynthetic pigments). To interpret such tracer profiles and estimate a mixing intensity one needs to apply an appropriate reactive-transport model (Officer and Lynch, 1983). A variety of mixing models are available, depending on the temporal and spatial length scale of the particle displacements (Meysman et al., 2003). However, the most widely used model to interpret downward tracer profiles is the diffusive mixing model, originally proposed by Goldberg and Koide (1962) and popularized by Guinasso and Schink (1975). Most applications of this diffusive model are based on the analytical solution derived for the case of constant porosity. The significance of this assumption is the focus of the present study. The diffusive mixing model assumes small-scale, random particle displacement, and implements Fick’s law to

962

F.J.R. Meysman et al. 71 (2007) 961–973

describe the particle flux (Boudreau, 1986; Wheatcroft et al., 1990; Meysman et al., 2003). Based on the diffusion analogy, the mass conservation equation for a radio-tracer typically reduces to the classical diffusion-reaction equation oC si o2 C s ¼ Db 2i  kC si ot ox

ð1Þ

where C si denotes the tracer activity (measured per g of dry sediment), Db is the mixing coefficient (cm2 yr1), and k the decay constant (yr1). The diffusive mixing model (1) has two clear advantages, which undoubtedly contribute to its popularity. Firstly, mixing is conveniently characterized by a single parameter, the mixing coefficient Db, which enables a straightforward comparison of the mixing intensity between different sediments. Secondly, the diffusion equation emerges in many disciplines of the natural sciences, and hence, its solution strategy is well documented (Carslaw and Jaeger, 1959; Crank, 1975). As a rule, analytical solutions are available for those conditions that are relevant to sediment mixing, e.g., steady-state profiles under a constant input flux, or time-dependent mixing of a pulse input (Duursma and Hoede, 1967; Guinasso and Schink, 1975; Officer, 1982; Boudreau, 1997). In the standard analysis of radionuclide depth profiles, these analytical solutions are applied in an inverse manner to estimate the values of the mixing coefficient (e.g., Robbins et al., 1977; Nozaki et al., 1977; DeMaster and Cochran, 1982; Cochran, 1985; Pope et al., 1996; Fornes et al., 1999; Henderson et al., 1999). As a result, a considerable database of Db coefficients has been built up over the years, documenting the mixing intensity in a variety of marine, estuarine and freshwater sediment locations. This database is then used to relate sediment mixing to other variables that control sediment biogeochemistry (Boudreau, 1994; Middelburg et al., 1997). These Db values have all been estimated based on the diffusive mixing Eq. (1), and accordingly, these numbers are all dependent on some important assumptions that underlie this model: (a) the mixing phenomenon can be described as a diffusive process; (b) the mixing intensity is constant over the mixed layer; (c) the radio-tracer is irreversibly adsorbed on the solid particles; (d) the solid volume fraction /s = 1  /f, where /f is the porosity, is assumed constant over the mixed layer. Here, our aim is to assess the effect of the last idealization in the classical procedure of quantifying sediment mixing. The surface layer of aquatic sediments exhibits appreciable gradients in the solid volume fraction due to compaction (Mulsow et al., 1998). Porosity typically decreases from 0.9 to 0.7 within the mixed zone of muddy sediments (Archer et al., 1989; Iversen and Jorgensen, 1993), implying a threefold increase of the solid volume fraction. The constant-porosity assumption effectively implies that this effect of compaction is ignored. An important motivation for our study is to retain the internal consistency of the current database of Db values. Due to increased computing power and improved access to numerical software, recent studies have started to explore

biodiffusion models that account for variable porosity (Berg et al., 2001; Nie et al., 2001; Alperin et al., 2002). When applying different mixing models, i.e., a variable porosity versus constant porosity, one is bound to obtain different mixing coefficients. To allow comparison, one should investigate the influence of porosity gradients on the estimation of Db values. 2. Model formulation 2.1. Conservation equations We consider a sediment environment of variable porosity, where the solid particles are subject to both compaction and diffusive mixing. To describe the reactive transport of a solid tracer in this sediment, one needs two separate conservation statements: one for the solid phase as a whole, and one for the specific tracer. Under the assumption of a constant density qs for the solid phase, these two equations are respectively given by (see Meysman et al., 2005 for a detailed derivation) o s 1 o s o ½/  ¼  s ½J  ¼  ½/s vs  ot q ox ox o  s s 1 o  s / Ci ¼  s J  k/s C si ot q ox i   o  s s s o oC si s v / Ci þ / Db ¼  k/s C si ox ox ox

ð2Þ

ð3Þ

where /s is the solid volume fraction (SVF), i.e., one minus the porosity, Js is the total solid flux (g cm2 yr1), J si is the total tracer flux (mole, Bq or dpm cm2 yr1), and ms is the mass-averaged velocity of the solid phase (cm yr1). The x-coordinate represents depth and the corresponding x-axis points downwards from the sediment–water interface. Note that in the past there has been some discussion on the proper formulation of component mass balance (3). The issue was whether the solid volume fraction /s should be placed inside or outside the oC/ox term in (3). Recently, it has been shown that /s should always be placed outside when using the mass-averaged velocity vs as the reference velocity (Meysman et al., 2005). The solution of the partial differential Eqs. (2) and (3) requires the specification of proper boundary conditions. To this end, we assume that the total flux Fs(t) of solid matter arriving at the sediment–water interface, also referred to as the sedimentation rate, is a known quantity, e.g., as determined from sediment traps. From the sedimentation rate (g cm2 yr1), one can directly calculate the advective velocity at the sediment–water interface (cm yr1), which is termed the sedimentation velocity x (Meysman et al., 2005), via x¼

F s ðtÞ qs /s0 ðtÞ

ð4Þ

where /s0 represents the SVF at the sediment–water interface (SWI). The parameter /s0 is typically referred to as

Influence of porosity gradients on sediment mixing

963

the ‘‘stress-free’’ solid volume fraction, i.e., the SVF that the sediment attains in the absence of any overlying weight. The flux of material arriving from the overlying water column should equal the solid flux transported downwards into the sediment. Formally, this leads to the boundary condition  F s ðtÞ ¼ qs /s0 x ¼ J s x¼0 ¼ ðqs /s vs Þjx¼0 ð5Þ

both Fs and /s0 are constant, the sedimentation velocity x also becomes invariant with time; see Eq. (4). When implementing condition (9) into the solid phase conservation Eq. (2), one observes that the total solid mass flux becomes invariant with depth

If C si;input denotes the tracer activity of the incoming material (measured per g of dry sediment), the tracer flux F si can be written in terms of the total flux Fs(t) as

From this, it directly follows that at any given point in the sediment

F si ðtÞ ¼ C si;input F s ðtÞ ¼ qs /s0 xC si;input

ð6Þ

o s o ½J  ¼ ½qs /s vs  ¼ 0 ox ox

vs ¼

F s /s0 ¼ x qs /s

ð10Þ

ð11Þ

The tracer flux (6) should equal the flux transported down from the sediment–water interface into the sediment. Formally, this leads to the second boundary condition at the sediment surface   F si ðtÞ ¼ qs /s0 xC si;input ¼ J si  x¼0   oC s  s s s s s s ¼ q v / C i  q / Db i  ð7Þ ox x¼0

Expression (11) is the sought-after expression for the massaveraged velocity, which completes our model statement. The adoption of steady-state compaction significantly simplifies the conservation equation for a solid component. Implementing (9) and (11), the tracer conservation Eq. (3) simplifies to   s o oC si o  s s oC i s / Db C  k/s C si / ¼ ð12Þ  x/s0 ox ox i ot ox

Finally, as Eq. (3) comprises a second-order differential equation, one also requires a boundary condition at the lower boundary Lmix of the model domain. In the conventional procedure of estimating mixing coefficients, the data are fitted by a tracer profile that is described by a single exponential function. The assumptions underlying this solution are that the sediment has a constant porosity and that it is mixed homogeneously over the semi-infinite domain (Lmix = 1). In reality however, sediments are only mixed to a finite depth, typically down to 10 cm (Boudreau, 1998). Therefore, we assume that the sediment is mixed homogeneously down to a finite depth Lmix. At this depth, diffusive transport should vanish, thus leading to the Neumann-type boundary condition  s  oC i  ¼0 ð8Þ ox x¼Lmix

Similarly, the boundary condition (7) at the sediment water–interface reduces to    s  /s Db oC si   C i x¼0  s ¼ C si;input ð13Þ /0 x ox x¼0

The equation set (2)–(8) forms a closed mathematical problem statement, provided that the mass-averaged velocity vs is specified. When modelling radio-tracer profiles, vs is typically constrained via the assumption of steady-state compaction. 2.2. Steady-state compaction

The actual solid volume fraction /s that should be substituted in Eqs. (12) and (13) depends on the specific model formulation (the depth-dependent SVF profile in the variable-porosity model; the average SVF in constant-porosity model; see below). 2.3. Solid volume fraction profile When adopting steady-state compaction, the solid volume fraction is no longer treated as an unknown variable. Rather, a predefined depth profile is used as a model forcing function. SVF profiles typically exhibit a smooth and regular shape in the upper few decimeters of marine sediments (Mulsow et al., 1998; Boudreau and Bennett, 1999). In most cases, the depth dependency of the SVF is well represented by the exponential relation    s  x s s s / ðxÞ ¼ /1  /1  /0 exp  ð14Þ Lpor

ð9Þ

where /s0 denotes the stress-free SVF introduced previously, /s1 is the asymptotic SVF at depth, and Lpor represents the attenuation length. Mulsow et al. (1998) showed that Eq. (14) could represent >90% of the variance in SVF profiles obtained from the Nova Scotian shelf. Based on (14), the average SVF within the mixing zone can be calculated as    s

Lmix s s s Lpor /AVG ¼ /1  /1  /0 1  exp  ð15Þ Lmix Lpor

The constraint set by Eq. (9) should particularly hold for the solid volume fraction /s0 at the sediment surface. When

The attenuation length Lpor represents the typical length scale over which the SVF increases with depth, or for that

The assumption of steady-state compaction involves two conditions (Meysman et al., 2005). Firstly, the flux of solid material to the sediment, i.e., the sedimentation rate Fs, should remain invariant with time. Secondly, the SVF profile also should not change with time, i.e., o/s ¼0 ot

964

F.J.R. Meysman et al. 71 (2007) 961–973

matter, the porosity decreases with depth. In order to evaluate mixing and compaction regimes at different sites, this length scale Lpor should be compared to the depth of the mixed zone Lmix. This can be done via the scaling ratio b¼

Lpor Lmix

ð16Þ

Porosity gradients are only important within the mixed layer when 0 < b < 1. When b P 1, the porosity still decreases to the same asymptotic value /f1 ¼ 1  /s1 , but this change then takes place over a larger horizon, and so, most of porosity decrease occurs outside the mixed zone. In a similar fashion, we can characterize the relative increase in the SVF due to compaction by a second dimensionless quantity v¼

/s1  /s0 /s1

ð17Þ

When v fi 0 there is no prominent change in the SVF with depth; v fi 1 occurs when the porosity at the SWI is very high. Using both dimensionless quantities (16) and (17), the SVF depth relation (14) can be re-arranged into the scaled form  s ðxÞ ¼ 1  v expðx=bÞ /

ð18Þ

where the SVF itself has been rescaled to its asymptotic val s ¼ /s =/s . Equally, the absolute sediment depth ue, i.e., / 1 x has been replaced by the relative position within the mixed zone, i.e., x ¼ x=Lmix . Eq. (18) confirms that the shape of SVF profiles within the mixed zone is determined by the two dimensionless parameters v and b.

constant-porosity mixing model (CP). In order to perform a sensitivity analysis, it is useful to bring this equation set into a dimensionless form. To this end, we introduce the following dimensionless parameters Dadif ¼

kL2mix Db

ð21Þ

Daadv ¼

kLmix x

ð22Þ

These quantities are usually referred to as the diffusive and advective Damkohler numbers (Boudreau, 1997). They respectively provide the ratio of the time scale of diffusive and advective transport to the reaction time scale of the radio-tracer. The Peclet number is obtained as the ratio of both Damkohler numbers Pe ¼

xLmix Dadif ¼ Db Daadv

ð23Þ

Using these dimensionless quantities, and scaling the tracer  ¼ C s =C s activity as C i i;input , the conservation Eq. (19) can be re-arranged into the generic form    o oC oC  ¼0 a2  a0 C ð24Þ  a1 ox ox ox where the coefficients are defined through (‘‘CP’’ denotes constant porosity) 9 /sAVG 1bv½1exp ð1=bÞ Db > aCP ¼ ¼ s 2 > dif 2 /0 kfLmix g Da ð1vÞ > = CP x 1 a1 ¼ Lmix k ¼ Daadv ð25Þ > > > /sAVG 1bv½1exp ð1=bÞ CP ; a0 ¼ /s ¼ 1v 0

3. Two mixing models 3.1. Constant-porosity model In a first approach, we disregard all porosity variations within the mixed zone. To this end, we can rewrite all equations in terms of the average solid volume fraction /sAVG over the mixed zone. The steady-state form of the conservation Eq. (12) thus becomes   o oC si o  s s /AVG Db C  k/sAVG C si ¼ 0 ð19Þ  x/s0 ox ox i ox Note that this equation contains both the stress-free SVF /s0 and the average value /sAVG over the mixed layer. These two different solid volume fractions are present because x is defined as the advective velocity at the sediment–water interface, and not as the average advective velocity within the mixed layer. In a similar fashion as (19), the upper boundary condition (13) simplifies to  s   s  oC i  s s  x/0 C i x¼0  Db /AVG ¼ x/s0 C si;input ð20Þ ox x¼0 while the lower boundary condition (8) remains unchanged. The conservation Eq. (19) subject to the boundary conditions (20) and (8) forms the statement of the

In a similar fashion, we can transform the boundary conditions (20) and (8) to the dimensionless form   s   s  oC i    ¼ a1 a2 þ a1 C ð26Þ i x¼0  ox x¼0   s  oC i  ¼0 ð27Þ ox x¼1 In the case of constant porosity, the coefficients ai in the differential Eq. (24) and the boundary conditions (26) and (27) have constant values. The solution to the equation system (24)–(27) is well-known (Boudreau, 1997), and constitutes a sum of two exponential terms  s ðxÞ ¼ M 1 exp ðd1xÞ þ M 2 exp ðd2xÞ C i

ð28Þ

The explicit expressions for the attenuation constants dj and the coefficients Mj are provided in Appendix A. 3.2. Variable-porosity model In a similar way as was done for the constant-porosity model, we can now derive the dimensionless form of the tracer equation when accounting for a variable, depth-dependent porosity (VP model). Upon substitution of the explicit depth relation for the solid volume fraction (18),

Influence of porosity gradients on sediment mixing

965

the steady-state version of the tracer conservation Eq. (12) becomes   o oC oC s / ðxÞDb  x/s0  /s ðxÞkC ¼ 0 ð29Þ ox ox ox

Table 1 Parameter set representing a typical sediment environment with a strong porosity gradient within the mixed layer

By proper scaling, this differential equation can be re-arranged to the exactly same generic form (24) as for the case of the constant-porosity model. This scaling of the boundary conditions (21) and (8) provides expressions of the same form as (26) and (27). The difference between the variable and constant-porosity models lies in the expressions for the coefficients ai, which are now depth-dependent (‘‘VP’’ denotes variable porosity) 9 s 1v expðx=bÞ Db > aVP xÞ ¼ //ðxÞ s 2 ¼ dif ð1vÞ > 2 ð Da kðLmix Þ > 0 = x 1 aVP ð x Þ ¼ ¼ ð30Þ adv 1 Lmix k Da > > s > / ðxÞ 1v expð x =bÞ ; aVP xÞ ¼ / s ¼ 0 ð 1v The complete model statement of the variable-porosity model now forms a boundary value problem, composed of the tracer Eq. (24), the scaled boundary conditions (26) and (27), and the explicit expressions for the ai coefficients (30). Boundary value problems with non-constant coefficients are more intricate to solve, and frequently, one has to resort to numerical integration routines to obtain an approximate solution. Here however, we were very fortunate to be able to develop an analytical solution to this problem, as documented in Appendix B. The solution is of the form  s ðxÞ ¼ M 1 expðd1xÞF ðh11 ; h12 ; h13 ; v expðx=bÞÞ C i þ M 2 expðd2xÞF ðh11 ; h12 ; h13 ; v expðx=bÞÞ

ð31Þ

where F(a,b;c;y) denotes the hypergeometric function, which is a special, multi-parameter function that can be defined in the form of a hypergeometric series, and which is available in many software packages, such as Excel, R, MAPLE and MATLAB. The derivation of (31), as well as explicit expressions for the coefficients hjk in terms of the ai-coefficients (30) are provided in Appendix B. The hypergeometric function is not novel in early diagenetic modelling theory. Boudreau (1986) obtained it as the solution of the tracer conservation equation with depth-dependent mixing coefficients. 4. Results and discussion

Symbol

Value

Units

Stress-free porosity

/f0

0.9

—

Porosity at infinity Porosity profile attenuation

/f1 Lpor

0.5 5

— cm

Average porosity Mixing depth Sedimentation velocity Mixing coefficient Input tracer activity Tracer decay constant Solid phase density

/fAVG Lmix x Db C si;input k qs

0.69 10 0.1 1 1000 0.031 2.5

— cm cm yr1 cm2 yr1 Bq kg1 yr1 g cm3

The sedimentation velocity and mixing coefficient are typical values for marine continental shelf sediments. The tracer decay constant is that of 210 Pbxs.

average mixing intensity over a timescale of 100 yr (Alperin et al., 2002). The depth of the mixed layer was fixed at Lmix = 10 cm, which is consistent with the average mixing depth observed in marine sediments (Boudreau, 1998). The porosity profile has the values /f0 ¼ 0:95 at the surface and /f1 ¼ 0:5 at depth, while the attenuation depth is Lpor = 5 cm, so that a clear porosity gradient is present within the mixed layer (resulting in a length scale ratio b = 0.5). Using expression (15), the average porosity within the mixed layer is /fAVG ¼ 0:69. Adopting a sedimentation velocity x = 0.1 cm yr1, this corresponds to a sedimentation rate F s ¼ ð1  /f0 Þxqs ¼ 120 g m2 yr1 , using an average solid density qs = 2.5 g cm3. The input tracer activity is set to the representative value C si;input ¼ 1000 Bq kg1 . Note that the choice of input tracer activity is arbitrary, since the solution C si ðxÞ of both the CP and VP models will simply scale with C si;input . Fig. 1 compares the simulated tracer depth profiles obtained by the CP and VP models. A first observation is that

0

0

Excess 210Pb (Bq kg–1) 100 200 300

400

2 depth (cm)

0

Parameter

4

6

4.1. Comparison of CP and VP models 8

To illustrate the intrinsic differences between the CP and VP models, they were solved using the same set of parameters (Table 1), exemplifying a typical marine sediment that has a strong porosity gradient within the mixed layer. The radio-tracer that is modelled is excess 210Pb (denoted 210 Pbxs), which is a natural product of the 238U decay series. 210Pbxs has a decay constant k = 0.031 yr1 (half-life 22.3 yr), and is the standard radio-tracer to estimate the

10

Fig. 1. Steady-state activity profiles for excess 210Pb generated by the constant-porosity (CP) model (solid line), and the variable-porosity (VP) model (dashed line) for a typical sediment environment with a strong porosity gradient within the mixed layer. Both profiles are generated with exactly the same porosity profile and tracer input (parameter set as summarized in Table 1).

966

F.J.R. Meysman et al. 71 (2007) 961–973

both solutions have a very similar shape, as both models generate a smoothly decreasing ‘‘exponential-like’’ profile. Note however, that only the CP model truly generates exponential solutions, while the VP model solution is a hypergeometric function (Appendix B). Given the error bars on actual tracer data (see below), it would be very difficult, if not impossible, to actually determine which of the two curves (exponential or hypergeometric) provides a better fit to the data. Moreover, the similar shapes of the CP and VP model solutions prohibits the use of tracer data (that are expressed per mass of sediment) to assess the presence and importance of porosity effects. From a given tracer profile alone, one cannot derive whether significant porosity gradients are present, let alone the degree to which such gradients will affect the estimation of the mixing intensity. A second striking feature is that the VP model predicts larger tracer activities than the CP model. Nonetheless, the input of tracer to the sediment F si is exactly the same for both models. To explain the difference, one can look at the average tracer activity within the mixed layer Z Lmix 1 s  Ci ¼ C si ðxÞdx ð32Þ Lmix 0 and its relation to the total tracer inventory. For a tracer with a fixed decay constant, the tracer flux should always equal the product of the decay constant and tracer inventory Z 1  F si ¼ kI si  k qs /s ðxÞC si ðxÞdx ð33Þ 0

Eq. (33) emphasizes that tracer inventory I si is not only governed by the tracer profile C si ðxÞ, but also by the depth dependence of the solid volume fraction, and possibly the depth dependence of the solid density (which is here, however, assumed constant). Moreover, the total inventory I si results partly from tracer within the mixed layer, and partly from tracer buried in deeper layers, I si ¼ I si;mix þ I si;deep Z Lmix Z ¼ qs /s ðxÞC si ðxÞdx þ 0

1

qs /s ðxÞC si ðxÞdx

ð34Þ

Lmix

Eqs. (33) and (34) can now be applied separately to the CP and VP models. One problem is that both the CP and VP models only describe the mixed zone, and thus, I si;deep cannot be calculated. However, we know that the deep contribution I si;deep should match the flux to the deeper sediment, which is given by x/s0 C si ðx ¼ Lmix Þ in both models. Furthermore, one observes that the VP and CP solutions approach each other at the bottom of the mixed layer (see Fig. 1). From these two premises, one can conclude that the deep contribution I si;deep will be similar in both the CP and VP models. Adopting this approximation and knowing that the flux F si imposed on both CP and VP models is the same, Eqs. (33) and (34) directly provide the following relation

Z Lmix  CP  1 C CP C i ðxÞdx i Lmix 0 Z Lmix  s  1 / ðxÞ=/sAVG C VP ¼ i ðxÞdx Lmix 0 Z Lmix  VP  1
ð35Þ

The inequality (35) demonstrates that the VP model predicts higher tracer activities than the CP model for the same tracer input. In essence, the average tracer activities  CP and C  VP are not the same because by the average of C i i the product of solid volume fraction and tracer activity does not necessarily equal the product of the respective  CP is always smaller than averages. More specifically, C i VP  , because the weighing factor ½/s ðxÞ=/s  is negatively C i AVG correlated with the tracer profile C VP i ðxÞ. In the upper sediment layers, tracer activities are higher than average, but solid volume fractions are lower than average. In deeper sediment layers, the situation is the other way round. Accordingly, one can expect the greatest deviation between VP and CP models, when the porosity and the tracer activity profiles show about the same depth attenuation (generating the highest correlation between the profiles). 4.2. Inverse modelling of a

210

Pbxs—porosity data set

In a second step, differences between the CP and VP models are examined using actual field data. To this end, both models are applied to combined porosity and excess 210 Pb (further denoted as 210Pbxs) data published by Alperin et al. (2002). This dataset was obtained in marine sediments on the continental slope near Cape Hatteras, North Carolina (USA), and is summarized in Table 2. Sediment cores were collected from cross-slope transects, resulting in a total of 12 stations ranging from 200 to 1000 m depth. This dataset is well suited for the present purpose, as it combines substantial porosity gradients within the mixed layer with sizeable mixing activities (Db values ranging from 0.2 to 5 cm2 yr1). Details on the study site, sediment collection, porosity determination and radiochemical analysis can be found in Alperin et al. (2002). Alperin et al. (2002) have already examined this dataset with a sophisticated model that accounts for variable porosity. Their mixing model includes two mixed layers with different mixing intensities, and hence, it is more complex than the VP model presented here. In total, the Alperin et al. (2002) model contained four adjustable parameters (the mixing intensity in the deep mixed layer, the sedimentation velocity, and two mixing depths), which were optimized simultaneously to provide the best fit to the combined tracer profiles of 210Pbxs and 239,240Pu (the latter tracer constrains the sedimentation velocity). For the present purpose, such model complexity is unnecessary. Our prime goal is to compare the adjustment of the Db coefficient between CP and VP models. Therefore, we simplify the Alperin et al. (2002) mixing model by ignoring the

Influence of porosity gradients on sediment mixing

967

Table 2 Porosity and mixing data from marine sediments on the continental slope near Cape Hatteras, North Carolina (USA) Station

Depth (m)

/f0

/f1

Lpor (cm)

Lmix (cm)

x (cm kyr1)

Inventory (mBq cm2)

1 2 3 4 5 6 7 8 9 10 11 12

212 296 325 455 490 532 544 607 746 752 850 1004

0.598 0.722 0.674 0.781 0.661 0.770 0.827 0.891 0.815 0.860 0.916 0.904

0.461 0.534 0.463 0.617 0.595 0.664 0.628 0.748 0.630 0.751 0.747 0.769

4.8 8.5 4.2 8.6 14.3 2.7 4.7 3.6 6.1 5.6 3.4 3.2

11 20 13 8.8 9.8 9.8 13 22 20 27 21 20

47 159 49 61 131 146 269 301 850 205 181 24

1040 2610 1540 1130 2100 2050 3380 3570 10,900 4210 4340 2440

Values are taken from Tables 2 and 3 in Alperin et al. (2002). Sediment cores were collected from cross-slope transects in 12 stations from 200 to 1000 m depth. The values for the sedimentation velocity x and the mixing depth Lmix were originally reported as x1 and L2 in Table 3 in Alperin et al. (2002).

shallow zone near the surface that is intensively mixed. This way, we only retain a single mixed zone. Furthermore, we do not constrain the sedimentation velocity x and the mixing depth Lmix by inverse modelling, but simply use the fitted values obtained by Alperin et al. (2002) as input parameters (Table 2). As a result, we only retain one adjustable parameter in each mixing model, i.e., the mixing VP intensity DCP b in the CP model, and the mixing intensity Db in the VP model. These parameters have been optimized by fitting the respective model solutions (28) and (31) to the tracer data. The goodness-of-fit for the modelled profiles is determined by maximizing the function 2 P  data  C model j j Cj R2 ¼ 1  P  ð36Þ 2 data  j data C  C j j where C data and C model are respectively the measured and j j model-derived tracer activities at each sampled depth j,  data is the mean of all measured tracer activities. and C j Expression (36) is the same as that of Alperin et al. (2002), thus allowing a straightforward comparison. The results of the above inverse modelling procedure are summarized in Fig. 2 and Table 3. The average difference between the Db coefficients determined by the CP and VP models is 16%. The minimal deviation is 2%, while the maximal difference is 31%. In two thirds of the cases (8 out of 12 cores), the mixing coefficient estimated by the VP model is larger than that determined by the CP model. This is consistent with the results of the previous section, where we found that when keeping the tracer input and Db values the same, the VP model typically predicts larger tracer activities than the CP model. This is also confirmed by the modelled tracer profiles in Fig. 2. Accordingly, for the same tracer input to the sediment, but keeping the tracer profiles the same, the VP model will in general need a larger mixing coefficient to match the data as compared to the CP model. This trend of the VP model predicting higher mixing coefficients is not absolute. In four cores, the CP model

value. These deviations from generated the highest DCP b the general rule are probably due to the interplay of strong model idealizations and ‘‘noise’’ within the experimental data. As noted in the introduction, there are other strong assumptions underlying the VP and CP models besides the treatment of porosity. The most important assumptions are that mixing is a truly diffusive mechanism, and that the mixing intensity remains constant with depth. The actual mixing process will violate these assumptions to some extent. For comparison, we have included in Table 3 the mixing coefficients that were estimated by the Alperin et al. (2002) model for the deep mixing zone. These are typically lower than the Db values estimated by the CP and VP models presented here. The reason for this is the inclusion of the extra shallow mixing zone in the Alperin et al. (2002) model. The extent of this zone is indicated in Fig. 2 by the horizontal dotted lines. This zone is intensively mixed, resulting in lower tracer activity values near the surface. Since the Alperin et al. (2002) model considers two zones, these shallow data points are not accounted for in the estimation of the ‘‘deep’’ Db value. In our simulations however, only a single mixed zone is modelled. Hence, the ‘‘wellmixed’’ data points near the SWI are ‘‘weighed’’ by the inverse procedure, thus increasing the Db estimates in our simulations. Table 3 also compares the goodness-of-fit coefficients for the Alperin, the CP and the VP models. The R2 coefficient in the VP model seems to be systematically higher than that of the CP model (11 out of 12 cores). This is not surprising, given the fact that the CP model only uses one parameter to describe porosity, while the VP model uses three. Given this increase in model complexity, the gain in terms of goodness-of-fit seems actually quite modest. In a similar fashion, the goodness of fit obtained by the Alperin et al. (2002) model is mostly higher than the one obtained by the VP model (7 out of 12 cores). Again, this is straightforwardly explained by the increased complexity of the Alperin model. One should also note that the R2 values reported in Alperin et al. (2002) are calculated for the complete data

968

F.J.R. Meysman et al. 71 (2007) 961–973 Excess 0 50 0

210

–1

Excess 210Pb (Bq kg–1) 0 200 400 0

Pb (Bq kg ) 100 150

Excess 0 0

210

–1

210

Pb (Bq kg ) 100 200

1

2

2

2

5

10

Depth (cm)

6

Depth (cm)

Depth (cm)

Depth (cm)

4 4

–1

Excess Pb (Bq kg ) 0 200 400 600 0

6 8

3 4 5 6

8

15

10

10

ST 1

7

12

ST 2

8

ST 3

ST 4

20 210

Excess 0 100 0

–1

Excess 210Pb (Bq kg–1) 0 200 400 0

Pb (Bq kg ) 200 300

Excess 0 0

210

–1

210

Pb (Bq kg ) 200 400

–1

Excess Pb (Bq kg ) 0 500 1000 0

2 2

2

5

6

4

6

Depth (cm)

4

Depth (cm)

Depth (cm)

Depth (cm)

4 6 8

10

15 10 8

8 12

ST 6

Excess 210Pb (Bq kg–1) 0 200 400 600 0

Excess 0 0

210

–1

Excess

Pb (Bq kg ) 500 1000

0 0

5

Depth (cm)

Depth (cm)

Depth (cm)

10 15

ST 9

25

–1

ST 10

20

5

10

20

ST 8

Excess 210Pb (Bq kg–1) 0 500 1000 0

Pb (Bq kg ) 500 1000

15

20

15

210

5

5

10

20

ST 7

Depth (cm)

ST 5

10

15 ST 11

ST 12 20

Fig. 2. Model fits to field data for 12 stations off Cape Hatteras. Excess 210Pb data (markers) are from Alperin et al. (2002); descriptive parameters for each station are provided in Table 2. The best fit of the variable-porosity (VP) model to the data is shown (dashed line). This then provides a best-fit estimate DVP b for the mixing coefficient (values are given in column 5 of Table 3). For comparison, the solution of the constant-porosity (CP) model is shown (solid line) that implements the same DVP b value obtained from the VP model. The procedure can also be reversed: now the CP model is fitted to the data, providing a different estimate DCP b for the mixing coefficient (values given in column 4 of Table 3; fits are not shown as they provide equivalent results). For reference, the horizontal dotted lines differentiate the two mixed layers that were included in the model of Alperin et al. (2002). As explained in the text, the CP and VP models implemented here only incorporate a single mixed layer.

profile, while the models presented here consider only data within the mixed layer. Similar to the approach taken here for the Cape Haterras dataset, Mulsow et al. (1998) applied various mixing models to combined porosity and 210Pbxs data for sediments from the eastern Canadian Margin. These authors actually distinguished three different mixing models. The first model assumed constant porosity mixing, and is identical to our CP model. The two other mixing models

accounted for depth-dependent porosity, and were respectively termed the interphase mixing model (IEM) and intraphase mixing model (IAM). The difference between these latter two models was related to the question of whether or not the mixing process affects the porosity profile. Mathematically, this difference concerns the location of the solid volume fraction in the diffusive term of Eq. (3): IAM incorporates the diffusive flux as J si ¼ /s Db oC=ox (mixing coefficient outside the differential), while IEM states the

Influence of porosity gradients on sediment mixing

969

Table 3 2 CP CP Bioturbation coefficients Db (cm2 yr1), deviation jDVP b  Db j=Db (%), and goodness of fit coefficients R obtained from the application of the constantporosity (CP) and variable-porosity (VP) models to the 12 continental shelf sites near Cape Hatteras (as described in Table 2) Station

Depth (m)

DAlperin b

DCP b

DVP b

CP CP jDVP b  Db j=Db ð%Þ

R2Alperin

R2CP

R2VP

1 2 3 4 5 6 7 8 9 10 11 12

212 296 325 455 490 532 544 607 746 752 850 1004

0.45 0.40 2.00 0.30 0.65 0.60 1.20 0.40 5.20 1.30 2.80 0.55

1.04 0.776 1.75 0.178 0.950 0.427 3.08 0.592 3.82 1.49 2.61 0.886

1.00 1.02 1.68 0.233 0.971 0.412 3.32 0.422 4.73 1.82 3.11 1.01

4 31 4 31 2 4 8 29 24 22 19 14

0.966 0.904 0.962 0.995 0.961 0.961 0.453 0.988 0.821 0.961 0.939 0.943

0.863 0.987 0.893 0.978 0.960 0.933 0.655 0.846 0.970 0.929 0.950 0.961

0.895 0.995 0.939 0.987 0.960 0.950 0.660 0.935 0.965 0.953 0.977 0.990

The values with superscript ‘‘Alperin’’ refer to the model results reported in Table 3 of Alperin et al. (2002).

diffusive flux as J si ¼ /s o½Db C=ox (mixing coefficient inside the differential). Recently, Meysman et al. (2005) have shown that there are indeed two different mixing regimes, i.e., interphase and intraphase mixing, but that in the case of steady-state compaction, both these regimes must be represented by the same model Eq. (3). In other words, the mixing coefficient Db should be always placed inside the differential, and the effect of interphase versus intraphase mixing should be incorporated through different porosity profiles, rather than via different model equations. So for the Mulsow et al. (1998) analysis, this implies that the reported Db values obtained from the IEM model are not valid. However, the IAM model presented by Mulsow et al. (1998) is a valid model, and is identical to our VP model. Table 4 shows the resulting Db values obtained from application of the IAM/VP and the corresponding CP model to the Mulsow et al. (1998) data. These results are very similar as those observed above. The average difference between the Db coefficients determined by the CP and VP models is 27%, slightly larger than the 16% above. The minimal deviation is 5%, the maximal difference is 50%. Again, in 6 out of 8 cases, the mixing coefficient estimated by the VP model is larger than the one determined by the CP model, confirming the trend that was observed for the Cape Haterras sediments. Table 4 CP CP Bioturbation coefficients Db (cm2 yr1) and deviation jDVP b  Db j=Db (%) obtained from the application of the constant-porosity (CP) and variableporosity (VP) models by Mulsow et al. (1998) to eight sediments from the eastern Canadian Margin Station

Depth (m)

DCP b

DVP b

CP CP jDVP b  Db j=Db ð%Þ

1 2 3 4 4 4 5 5

331 262 494 232 230 230 816 830

0.13 0.28 0.08 0.47 0.72 0.19 0.81 0.37

0.15 0.37 0.12 0.67 0.89 0.24 0.63 0.35

15 32 50 43 24 26 22 5

(Winter 93) (Winter 93) (Winter 93) (Summer 93) (Winter 93) (Summer 94) (Summer 93) (Summer 94)

Results reported in Tables 1, 3 and 5 of Mulsow et al. (1998).

4.3. Sensitivity analysis: deep sea, shelf, and coastal marine sediments In the previous section, we quantified the difference beVP tween DCP b and Db for sediments from the eastern Canadian margin and the continental slope off Cape Hatteras. Obviously, these only comprise two particular examples of sediment environments subject to compaction and mixing. To generalize our results, we performed a sensitivity analysis for a set of ‘‘virtual sediments’’ covering a range of conditions found in marine sediments. To this end, we first defined three characteristic settings, denoted as ‘‘deep sea’’, ‘‘shelf’’ and ‘‘coast’’, which are each characterized by a specific value for the sedimentation velocity x and the mixing coefficient Db (see the caption of Fig. 3). At each location, we adopted a fixed mixed depth of 10 cm, and implemented the same set of porosity profiles. Rather than being realistic for a particular location, these porosity profiles were constructed to investigate the maximal effect of porosity gradients on mixing estimates. Each porosity profile is characterized by the same strong decrease of porosity with depth, determined by /f0 ¼ 0:95 and /f1 ¼ 0:5. This fixed porosity decrease D/f ¼ /f0  /f1 ¼ 0:45 can be regarded as an upper bound on the actual porosity decreases found in surface sediments. Although the porosity decrease D/f is fixed, the porosity profile can still be adapted by changing the attenuation depth Lpor. To investigate different profile shapes, we varied the length ratio b continuously over the range from 0.1 to 10 at each of the three settings, thus mimicking a transition from sharply to gradually declining porosity profiles. We followed a two-step procedure to create and analyze this virtual data set. In a first step, we used the VP model to 210 generate ‘‘synthetic’’ data profiles C VP Pbxs i ðxÞ for VP employing a mixing coefficient Db that is characteristic of deep sea, shelf or coastal locations. In each case, the tracer input activity was fixed at C VP i;input  1 (as noted before, the problem simply scales with C VP i;input ). In a second step, we fitted the CP model to the artificial C VP i ðxÞ profiles,

970

F.J.R. Meysman et al. 71 (2007) 961–973

A

Deep sea setting 0

0

45

.2 .4 .6 .8

1

0

.2 .4 .6 .8

1

0

.2 .4 .6 .8

1

40 deviation (%)

depth (cm)

2 4 6

35

30 8

β=0.1

β=1

β=0.5

25

10

B

Continental shelf 0

0

.2 .4 .6 .8

1

0

.2 .4 .6 .8

1

0

.2 .4 .6 .8

.2

.4 .6 .8 length scaling ratio β

1

60

1

deviation (%)

50

depth (cm)

2 4

30

6

20

8

β=5

β=0.5

β=0.1

10 .1

10

C

2

4

6

8

10

Coastal setting 0

0

.2 .4 .6 .8

1

0

.2 .4 .6 .8

1

0

.2 .4 .6 .8

1

60 50 deviation (%)

2 depth (cm)

40

4 6

40 30 20

8

β=0.1

β=0.5

10

β=5 10 .1

2 4 6 8 length scaling ratio β

10

Fig. 3. Sensitivity analysis of the deviation between variable-porosity (VP) and constant-porosity (CP) models for three typical marine sediment 1 1 1 VP VP 2 1 2 1 2 1 environments. Deep-sea [DVP b ¼ 0:1 cm yr ; x = 0.01 cm yr ], Shelf [Db ¼ 1 cm yr ; x = 0.1 cm yr ] and Coast [Db ¼ 10 cm yr ; x = 1 cm yr ]. VP (Left) A ‘‘synthetic’’ dataset is first created from the output of the VP model with a given ‘‘input’’ Db value (diamonds). Next, the CP model is fitted to the synthetic data (bold line) and a corresponding ‘‘output’’ DCP b value is obtained. Three different porosity profiles (dashed line) are depicted for each setting, with shallow (b = 0.1), moderate (b = 0.5) or deep (b = 1 or 5) gradients in porosity. (Right) The deviation in the mixing intensity between CP and VP models (formula (37)) is plotted as a function of the length scale ratio b.

thus estimating a new mixing coefficient DCP b . This was done numerically by optimizing the goodness-of-fit funcmodel tion (36) where C data ¼ C VP ðxÞ ¼ C i CP ðxÞ. i i ðxÞ and C i The deviation of the ‘‘fitted’’ CP mixing coefficient from the ‘‘true’’ VP mixing coefficient is quantified via the ratio   VP CP  ð37Þ Db h ¼ 100  DVP b  Db The results are depicted in Fig. 3. For the deep sea setting with low mixing intensity (0.1 cm2 yr1), the deviation h ranges from 35% for a sharply declining porosity profile

(b = 0.1), reaches a maximum of 40% around a value b = 0.3, and then decreases for gradually declining porosity profiles. As anticipated above (see discussion in Section 4.1), the maximum deviation is obtained when the attenuation depth of the tracer profile coincides with the attenuation depth Lpor of the porosity profile. At this point, the correlation between the tracer and porosity profile is maximal, leading to the largest deviation. Both the shelf (1 cm2 yr1) and coastal (10 cm2 yr1) settings show a similar response, where the maximal deviation amounts

Influence of porosity gradients on sediment mixing

to 60% and is attained at the value b  3. Because of increased mixing, the attenuation depth of the tracer profile is greater in shelf and coastal settings, thus explaining the higher value of b where the maximal deviation is reached. Overall, the pattern is very similar in all three settings: the deviation h goes through a maximum when the attenuation depth Lpor of the porosity profile increases. This is as expected: for both Lpor fi 0 and Lpor fi 1, there are no porosity gradients within the mixed zone, and hence, the VP and CP models should converge to the same mixing estimate. 5. Summary and conclusions The mixing intensity in sediments is typically determined from the inverse fitting of the standard diffusion equation to tracer depth profiles, such as excess 210Pb. A crucial assumption underlying this approach is that the mixed layer has a constant porosity. Because surface sediments can exhibit appreciable decreases in porosity in the mixed layer, we have assessed the influence of a depth-dependent porosity on the mixing coefficient Db. To this end we first derived a general solution to the reactive transport equation for a radionuclide in sediments where porosity decreases exponentially with depth. This solution takes the form of a hypergeometric function, which generates tracer depth profiles that are very similar to the exponential curves of the classical constant-porosity model. When applying both constant-porosity (CP) and variable-porosity (VP) model solutions to the same tracer profile, our theoretical analysis of the governing tracer equations shows that the VP model will typically generate a larger mixing coefficient.

971

Analysis of 210Pbxs profiles from sediments off the coast of Cape Haterras (Alperin et al., 2002) and the eastern Canadian Margin (Mulsow et al., 1998) demonstrates that both the CP and VP models are equally capable of explaining the data. The resulting mixing coefficients estimated from the CP and VP models differ on average by 20%, with a maximum deviation up to 50%. A subsequent sensitivity analysis, which covered the spectrum of possible interactions between mixing and porosity gradients in marine sediments, provides an upper limit to the deviation of around 60%. Both data sets demonstrate the systematic trend towards higher mixing coefficients estimated by the VP model. Still, both field data and theoretical model analysis show that the systematic deviation from the one-to-one relation between VP and CP mixing coefficients can be considered as weak (Fig. 4). This conclusion is reinforced if we compare the deviation resulting from porosity gradients as defined in (37) with other sources of uncertainty (expressed via the coefficient of variation (COV), which is defined as the standard deviation multiplied by one hundred and divided by the mean). A first type of uncertainty is related to the measurement of tracer activities: typical COV values for 210Pbxs data are in the range between 1% and 10% (Legeleux et al., 1994; Santschi et al., 2001). Secondly, one has the uncertainty arising from the least squares fit to the data. Although seldom reported in tracer studies, COV values for Db’s obtained from fitting 210Pbxs profiles can be rather large, ranging beyond 100% (e.g., COV of Db’s in Anderson et al. (1988): range 6–590%, mean 74%; COV of Db’s in Mulsow et al. (1998): range 2–52%, mean 25%). Finally, one should account for the other model assumptions underlying the classical diffusive mixing model

Fig. 4. Comparison of the mixing coefficients estimated from field data by the variable-porosity (VP) and the constant-porosity (CP) models. Two sets of biodiffusion coefficients are plotted: our model estimates based on the data of Alperin et al. (2002) as given in Table 3, and the model results of Mulsow et al. (1998) as given in Table 4. No systematic deviation from the one-to-one line is observed.

972

F.J.R. Meysman et al. 71 (2007) 961–973

(infinitely small and frequent steps, constant mixing intensity, irreversible tracer adsorption; Wheatcroft et al., 1990; Meysman et al., 2003, 2005). Given these various types uncertainties, the influence of porosity gradients on Db values can be classified as rather modest, especially when compared to the COV of the least squares Db estimate. In other words, the constant-porosity mixing model does a respectable job in explaining the data, and hence, the idealization of a constant average porosity may be adopted without any particular sacrifice of accuracy in the calculation of mixing coefficients.

Using these substitution, Eq. (24) is transformed into

o2 W VP VP oW þ 2a2 k þ aVP 2  a1 2 ox ox k v 1 þ dif ðb  2lð1  vÞÞ expðx=bÞW ¼ 0 1v Da

aVP 2

ðA:7Þ

In a next step, we can also rescale the depth coordinate as z ¼ v expðx=bÞ

ðA:8Þ

Acknowledgments Many thanks to Marc Alperin and Larry Benninger for providing the excellent dataset. The authors acknowledge the financial support of the US Office of Naval Research (grants N00014-02-1-0107 and N00014-99-1-0063) and the program manager, Dr. J. Eckman. This research received additional support through a PIONIER grant to J.J.M. from the Netherlands Organization for Scientific Research (NWO, 833.02.2002). This is publication 3906 of the NIOO-KNAW (Netherlands Institute of Ecology).

and as a result, Eq. (A.7) becomes o2 W oW þ ½ð1  2b½k  lð1  vÞÞ  2zð1  kbÞ oz2 oz 2  ½2lb ð1  vÞ  bkW ¼ 0 ðA:9Þ

zð1  zÞ

Eq. (A.9) is the hypergeometric equation. Introducing following notations 9 2kbð2lbð1  vÞ  1Þ ffi > qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > > 2 dif > > 1  2kb þ 1 þ 4b Da =  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > b ¼ 12 1  2kb þ 1 þ 4b2 Dadif > > > > > ; c ¼ 1  2bðk  lð1  vÞÞ

a¼

Associate editor: David J. Burdige

Appendix A. Solution of constant-porosity (CP) model The attenuation constants dj in the solution to the constant porosity model (28) are given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi j dj ¼ c þ ð1Þ c2 þ Dadif j ¼ 1; 2 ðA:1Þ with the auxiliary parameter c = 0.5(a1/a2), the Damkohler number Dadif as defined in (21), and a1 and a2 as given by expressions (25). Similarly, the coefficients Mj in Eq. (28) are determined from the boundary conditions (26) and (27) as  2cdjþð1Þjþ1 exp djþð1Þjþ1 jþ1 j ¼ 1; 2 ðA:2Þ M j ¼ ð1Þ D where the auxiliary parameter D is introduced as D ¼ d1 ðd2  2cÞ expðd1 Þ  d2 ðd1  2cÞ expðd2 Þ

one can once more rewrite Eq. (A.9) as zð1  zÞ

þ N v1c expðð1  cÞx=bÞF ða þ 1  c; b þ 1  c; 2  c; v expðx=bÞÞ

ðA:4Þ

where k is a root of the quadratic equation l¼

Dadif 2Daadv

and the proper root is selected as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k ¼ lð1  vÞ  l2 ð1  vÞ2 þ Dadif

ðA:11Þ

Cðx; tÞ ¼ expðkxÞ½MF ða; b; c; v expðx=bÞÞ ðA:3Þ

We use the substitution

k 2  2lð1  vÞk  Dadif ¼ 0;

o2 W oW  abW ¼ 0 þ ½c  ða þ b þ 1Þz oz2 oz

Deriving the general solution of this equation and taking into account substitutions (A.4) and (A.8), we obtain following solution for our original Eq. (24)

Appendix B. Solution of variable-porosity (VP) model

CðxÞ ¼ W ðxÞ expðkxÞ

ðA:10Þ

where M and N are arbitrary coefficients, which are to be found from the boundary conditions. The substitution of (A.12) into Eqs. (26) and (27) leads to a system of linear algebraic equations, from which we find that

ðA:5Þ M¼

ðA:6Þ

ðA:12Þ

G22 ; D

N ¼

G21 ; D

D ¼ G11 G22  G12 G21

The G’s are respectively defined as

ðA:13Þ

Influence of porosity gradients on sediment mixing



 k v 0 F ða; b; c; vÞ 1 F ða; b; c; vÞ þ 2l 2lb   kb  ð1  cÞ G12 ¼ v1c 1  2lb

G11 ¼

 F ða þ 1  c; b þ 1  c; 2  c; vÞ þ

v 0 F ða þ 1  c; b þ 1  c; 2  c; vÞ 2lb



G21 ¼ kF ða; b; c; v expðb1 ÞÞ  b1 v expðb1 ÞF 0 ða; b; c; v expðb1 ÞÞ G22 ¼ v1c expðb1 ð1  cÞÞ  ½ðk  b1 ð1  cÞÞ  F ða þ 1  c; b þ 1  c; 2  c; v expðb1 ÞÞ  vb1 expðb1 Þ  F 0 ða þ 1  c; b þ 1  c; 2  c; v expðb1 ÞÞ ðA:14Þ where F is the hypergeometric function and F0 represents its first derivative. References Alperin, M.J., Suayah, I.B., Benninger, L.K., Martens, C.S., 2002. Modern organic carbon burial fluxes, recent sedimentation rates, and particle mixing rates from the Upper Continental Slope Near Cape Hatteras, North Carolina (USA). Deep-Sea Res. Part II-Top. Stud. Oceanogr. 49, 4645–4665. Anderson, R.F., Bopp, R.F., Buesseler, K.O., Biscaye, P.E., 1988. Mixing of particles and organic constituents in sediments from the continental shelf and slope off Cape Cod: SEEP-I results. Cont. Shelf Res. 8, 925–946. Archer, D., Emerson, S., Reimers, C., 1989. Dissolution of calcite in deepsea sediments—pH and O2 microelectrode results. Geochim. Cosmochim. Acta 53, 2831–2845. Berg, P., Rysgaard, S., Funch, P., Sejr, M.K., 2001. Effects of bioturbation on solutes and solids in marine sediments. Aquat. Microb. Ecol. 26, 81–94. Boudreau, B.P., 1986. Mathematics of tracer mixing in sediments: I. Spatially dependent, diffusive mixing. Am. J. Sci. 286, 161–198. Boudreau, B.P., 1994. Is burial velocity a master parameter for bioturbation. Geochim. Cosmochim. Acta 58, 1243–1249. Boudreau, B.P., 1997. Diagenetic Models and Their Implementation. Springer. Boudreau, B.P., 1998. Mean mixed depth of sediments: the wherefore and the why. Limnol. Oceanogr. 43, 524–526. Boudreau, B.P., Bennett, R.H., 1999. New rheological and porosity equations for steady-state compaction. Am. J. Sci. 299, 517–528. Carslaw, H.S., Jaeger, J.C., 1959. Conduction of Heat in Solids. Clarendon Press. Cochran, J.K., 1985. Particle mixing rates in sediments of the Eastern Equatorial Pacific—evidence from Pb-210, Pu-239, Pu-240 and Cs-137 distributions at Manop sites. Geochim. Cosmochim. Acta 49, 1195– 1210. Crank, J., 1975. The Mathematics of Diffusion. Oxford University Press. DeMaster, D.J., Cochran, J.K., 1982. Particle mixing rates in deep-sea sediments determined from excess 210Pb and 32Si profiles. Earth Planet. Sci. Lett. 61, 257–271. Duursma, E.K., Hoede, C., 1967. Theoretical, experimental and field studies concerning molecular diffusion of radioisotopes in sediments and suspended solid particles of the sea. Part A: Theories and mathematical calculations. Neth. J. Sea Res. 3, 423–457.

973

Fornes, W.L., Demaster, D.J., Levin, L.A., Blair, N.E., 1999. Bioturbation and particle transport in carolina slope sediments: a radiochemical approach. J. Mar. Res. 57, 335–355. Goldberg, E.D., Koide, M., 1962. Geochronological studies of deep-sea sediments by the Io/Th method. Geochim. Cosmochim. Acta 26, 417– 450. Guinasso, N.L., Schink, D.R., 1975. Quantitative estimates of biological mixing rates in abyssal sediments. J. Geophys. Res.-Ocean. Atmos. 80, 3032–3043. Henderson, G.M., Lindsay, F.N., Slowey, N.C., 1999. Variation in bioturbation with water depth on marine slopes: a study on the Little Bahamas Bank. Mar. Geol. 160, 105–118. Iversen, N., Jorgensen, B.B., 1993. Diffusion-coefficients of sulfate and methane in marine-sediments—influence of porosity. Geochim. Cosmochim. Acta 57, 571–578. Legeleux, F., Reyss, J.L., Schmidt, S., 1994. Particle mixing rates in sediments of the northeast Tropical Atlantic—evidence from Pb210(Xs), Cs-137, Th-228(Xs) and Th-234(Xs) downcore distributions. Earth Planet. Sci. Lett. 128, 545–562. Li, W.Q., Guinasso, N.L., Cole, K.H., Richardson, M.D., Johnson, J.W., Schink, D.R., 1985. Radionuclides as indicators of sedimentary processes in abyssal Caribbean sediments. Mar. Geol. 68, 187–204. Meysman, F.J.R., Boudreau, B.P., Middelburg, J.J., 2003. Relations between local, nonlocal, discrete and continuous models of bioturbation. J. Mar. Res. 61, 391–410. Meysman, F.J.R., Boudreau, B.P., Middelburg, J.J., 2005. Modeling reactive transport in sediments subject to bioturbation and compaction. Geochim. Cosmochim. Acta 69, 3601–3617. Meysman, F.J.R., Middelburg, J.J., Heip, C.H.R., 2006. Bioturbation: a fresh look at Darwin’s last idea. Trends Ecol. Evol. 21, 688–695. Middelburg, J.J., Soetaert, K., Herman, P.M.J., 1997. Empirical relationships for use in global diagenetic models. Deep-Sea Res. Part IOceanogr. Res. Pap. 44, 327–344. Mulsow, S., Boudreau, B.P., Smith, J.N., 1998. Bioturbation and porosity gradients. Limnol. Oceanogr. 43, 1–9. Nie, Y.H., Suayah, I.B., Benninger, L.K., Alperin, M.J., 2001. Modeling detailed sedimentary Pb-210 and fallout Pu-239, Pu-240 profiles to allow episodic events: an application in Chesapeake Bay. Limnol. Oceanogr. 46, 1425–1437. Nozaki, Y., Cochran, J.K., Turekian, K.K., Keller, G., 1977. Radiocarbon and Pb-210 distribution in submersible-taken deep-sea cores from project Famous. Earth Planet. Sci. Lett. 34, 167–173. Officer, C.B., 1982. Mixing, sedimentation-rates and age dating for sediment cores. Mar. Geol. 46, 261–278. Officer, C.B., Lynch, D.R., 1983. Determination of mixing parameters from tracer distributions in deep-sea sediment cores. Mar. Geol. 52, 59–74. Pope, R.H., Demaster, D.J., Smith, C.R., Seltmann, H., 1996. Rapid bioturbation in equatorial pacific sediments: evidence from excess Th234 measurements. Deep-Sea Res. Part II-Top. Stud. Oceanogr. 43, 1339–1364. Robbins, J.A., Krezoski, J.R., Mozley, S.C., 1977. Radioactivity in sediments of Great Lakes—post-depositional redistribution by deposit-feeding organisms. Earth Planet. Sci. Lett. 36, 325–333. Santschi, P.H., Guo, L., Asbill, S., Allison, M., Britt Kepple, A., Wen, L.-S., 2001. Accumulation rates and sources of sediments and organic carbon on the Palos Verdes Shelf based on radioisotopic tracers (137Cs, 239,240Pu, 210Pb, 234Th, 238U and 14C). Mar. Chem. 73, 125–152. Thomson, J., Brown, L., Nixon, S., Cook, G.T., Mackenzie, A.B., 2000. Bioturbation and holocene sediment accumulation fluxes in the NorthEast Atlantic Ocean (benthic boundary layer experiment sites). Mar. Geol. 169, 21–39. Wheatcroft, R.A., Jumars, P.A., Smith, C.R., Nowell, A.R.M., 1990. A mechanistic view of the particulate biodiffusion coefficient—step lengths, rest periods and transport directions. J. Mar. Res. 48, 177– 207.