Modelling early diagenesis of silica in non-mixed sediments

Deep-.Sru Rrsca~h. VoL 37. No. IlL pp. 15.,13-1567. lq~l. Prmtcd m Great Bt'tt~m.

01qM--Ol,ll~qO~.00 + 0.00 ~ LWOPerpmon Prc~ pk:

Modelling early diagenesis of silica in non-mixed sediments BERNARD P. BOUDREAU* (Received 15 November 1989: in revised form 30 April 1990; accepted 9 May 1990) Abstract--The diagenesis of dissolvable opal in non-mixed sediments of anoxie and dysaerobic basins is described by a set of coupled differential conservation equations. A semi-analytical solution is developed for the special case of steady-state constant-porosity diagenesis. This approach reduces the problem of predicting opal accumulation to one of solving a single nonlinear algebraic equation and provides an explicit equation for the critical opal flux below which no opal would be preserved. It also allows the derivation of a simple asymptotic model valid when a small fraction of incoming opal dissolves, i.e. <50%. Investigative calculations indicate that rapidly accumulating non-mixed sediments can easily preserve most of the opal that reaches them. In this case. the advective flux of ~lid sediment is great that it overwhelms the potential for dissolution. Bioturhation is not needed to promote accumulation, even when the opal content is modest. This contrasts with slowly accumulating (e.g. deep-sea) sediments wherein thc dissolution potential is great and most opal would dissolve were it not for mixing. Model fits to scdimcnt-porcwatcr data from the Carmen Basin and the San Pedro Basin indicate a value for the opal-dissolution ratc constant of between ItX)and 4(X)y- t. The San Pedro Basin data are particularly instructivc in this rcgard because the mode[ fit, ohtaincd by adjusting thc ratc constant, can be vcrificd by comparing the modcl-prcdictcd dissolvcd flux out of thcse sediments with an independently mc:|sured value. Application of the model to data from Peru Margin scdimcnts indicatcs that these cannot be in a stc:idy state, but arc adjusting to a changc in opal prt~luction that occurrcd about [(10 lip or a rcccnt suddcn deposition of sediments richcr in opal (c.g. by :t slide). This paper also discusses modilications and extensions that could bc made to the present model to account for multiple reactive types of opal. thermodynamic and kinetic inhibition of solubility equilibrium at depth and mineralogical transformations during late diagenesis

INTRODUCTION

T~E accumulation of biogenic opal in marine sediments constitutes the only major sink of the oceanic silica budget (DEMASTER,1979, 1981; LEDFORD-HOFFMANet al., 1986). The magnitude of opal accumulation in marine sediments is determined by a dynamic balance between supply, transport processes (e.g. burial and bioturbation) and dissolution. In sediments mixed by either biological or physical agents, these dynamics have been described by a variety of models that explicity account for both the distributions of dissolved silica and the solid opal with depth (e.g. BERNER,1974, 1980; SCmNKet al., 1975; SCHINK and GUtNASSO, 1977, 1980; WONG and GROSCH, 1978; BOUDREAU, 1990). In

"Department of Oceanography, Dalhousie University, Halifax, Nova Scotia B3H 4J l, Canada. 1543

1544

B . P . BOUDRF~U

particular, the studies by ScmNK et al. (1975) and SCmNK and G u I N ~ s o (1977. 1980) have convincingly illustrated the role of bioturbation in preserving opal in opal-poor sediments and have highlighted the crucial part played by the input rain of opal in forming siliceous oozes.

The studies listed above deal with diagenesis in sediments subject to biological or physical mixing. Yet, sediments wherein mixing is absent are also sites of opal accumulation, a process which sometimes can be rather intense (DEMaSTER, 1979, 1981). Such areas include regions of upwelling along continental margins (e.g, Gulf of California, Peru-Chile coast and Walvis Bay), and the deep relatively stagnant basins of the California Borderlands. Other non-mixed sediments, like those of the Black Sea and the Cariaco Trench, also store silica although in more modest amounts. Mathematical modelling of porewater silica and opal in non-mixed sediments has been far more limited than in their mixed counterparts. BOUCHEr (1984) and BeRELSONet al. (1987) use the solution of a simple diffusion-reaction model to describe the dissolved silica profiles and calculate fluxes from some California Borderland cores. There is certainly no equivalent for non-mixed sediments to the SCmNK et al. (1975) and ScmNg and GulNasso (1980) coupled opal-silica models. The existence of such a model would be quite useful, however. From a theoretical point of view. it would be possible to examine the conditions (i.e. parameter and input values) under which sediments could preserve opal to become an ooze or, conversely, be devoid of dissolvable opal. By coupling solid opal dissolution and the production of dissolved silica, the calculation of important geochemical variables such as flux become explicitly dependent on other observed geochemical properties and variables, such as the amount of opal preserved at depth. (This removes some of the arbitrariness associated with the normal curve-fitting technique.) It is further possible to judge if opal and silica distributions arc consistent with the asstmlcd conditions, e.g. steady state. Finally, sttch a model would idlow testing of laboratory-derived dissolution rate constants to see if they reproduce natural distributions. (The current lack of widely accepted laboratory rate constants precludes this particular usage, but this can change in the future.) In searching for a model for silica--opal diagcnesis in non-mixed sediments, one cannot simply let the biodiffusion coeffiicient in the diffusively mixed SCmNK and Gutr~asso (198(I) model go to zero, i.e. Du ---, 0, and hope to not encounter problems. The numerical method utilized by SCHINK and GtJINASSO (1980) tO solve their model is based on the centred finite-difference scheme, and this in known to be unstable if the ratio w h / D u becomes greater than two, where w is the burial velocity of solids and h the step size of the numcrical approximation (e.g. FIaDemO and VERONtS, 1977; GLASS and RODt, 1982; JaIN et al., 1984; SuYY, 1985). The ratio whlDB is known as the cell or grid Peclet number. As D~ ---, 0, it is inevitable that the condition for instability be encountered because a counter-balancing reduction of the step size is limited by numerical round-off. A solution requires that the Scmsg and GOtNASSO (1980) model be solved numerically with a different stable scheme or that an analytical solution be devised for the case with no mixing, i.e. Dt~ = 0. This paper adopts the second of these two approaches. It offers a general coupled model for the diagenesis of opal in non-mixed sediments (D u = 0), and then solves this model by a semi-analytical technique for the important special case of steady-state diagenesis. Calculations are performed to illustrate the manner in which the solution depends on the model parameters and inputs. Finally, the model is applied to data sets from three

Modellingearlydiagenesisof silica

1545

environments to help determine the processes that have contributed to the formation of the observed distributions and to illustrate the utility of the model. MODEL SPECIFICATION AND SOLUTION The present model deals with a three component system, i.e. dissolved silica, C (mM). solid opal, b (g cm -3 of total solids), and an inert phase, m (also g cm -3 of total solids).

which is most likely to be clay. In general, the increase (positive, negative or nil) of the dissolved silica concentration at a depth will result from the net interaction between production by dissolution and the transport by diffusion and advection of the porewater. This conservation is expressed mathematically as (BERNER. 1980)

Oq;C_ Ox O q;Ds~ ]

Oq;uC+q;ks(Cs-C) ~ ,

(1)

where x is the depth from the sediment-water interface (cm), ff is the porosity, D s is the tortuosity corrected diffusion coefficient for dissolved silica (cm 2 y-t), u is the advective velocity of the porewaters, k~ is the rate constant for opal dissolution (y-l), C~ is the "saturation" or asymptotic concentration of dissolved silica, and Po is the intrinsic density of opal (g cm -3 of opal). The opal dissolution kinetics represented by the term on the far right-hand side of equation (1) assumes a linear dependence on both the local saturation deficit, C~ - C, and the mass concentration of dissolvable opal. This assumed saturation functionality is consistent with experimental studies (fiorD, 1972; WOLLAST, 1974; LAWSONet al., 1978; KAMATANtand RILEY, 1979; KAMATANIet al., 1980; HURD and BtRDWmSTELL,1983). The linear b/pb-dcpcndencc originates from the idea that the specific surface area available for reaction is constant, i.e. the reactive sttrface area is proportional to the amount of opal. SCtIINKet al. (1975) and SClIINKand GUtNASSO(1980) have followed this treatment in their seminal papers. Ervz et al. (1982) fottnd this assumptkm quite adequate in describing their in situ dissolution data. liowcvcr, tlurD and Bmt}WmSTEt.L(1983) observed an inital increase in the spccilic surface area during their dissolution experiments, followed by a subsequent decrease its the experiments ran their course. KErn (1982), in modelling dissolution in the CaCO 3 system, created a dissolution model whereby the surface area was proportional to the number concentration of particles rather than their mass concentration. His model also predicted an early increase in the surface area per unit volume before a fall with progressive dissolution. BOUDrEAU (1985), again working in the carbonate system, investigated the differences between predictions made by models with arbitrary nonlinear equations linking b/pb to reactive surface area per unit volume of sediment, i.e. quadratic, cubic and piecewiselinear chapeau functions. This study showed that the choice of b/~-functionality had little effect on the predicted results when dissolution generated an initial increase in surhme area. For this reason and to achieve consistency with most previous studies, I adopt the linear b/pb-dependence with the caveat that it may engender some (probably small) error. All changes in opal concentration at a particular depth must be due to any imbalance between disappearance by dissolution and transport by burial, as mixing is absent, a(l - qo)b = dlt

0(1 - tp)wb 3x

rq~ks(C,- COb , Pb

(2)

1546

B.P. Bocmtr~u

where w is the burial velocity of the sediment, (I - ~) is the solid volume fraction, and y is a constant for unit conversion (0.0961 g mmol-t). The change in the concentration of the inert component is simply related to the change in the burial flux of this species with depth. a(l - ¢)m _

a(l - ~).'m

Ot

(3)

Ox

In addition, ScH1ro: and GU~NASSO(1980) correctly observe that total space must be conserved, which requires that b

m

Pb

Prn

-- + --

= I.

(4)

where p,, is the intrinsic density of the inert material. It is my intention to examine the solution and properties of the model, equations (1)-(4), for steady-state constant-porosity diagenesis in non-mixed sediments. It might be argued that this may not be a reasonable limiting case because sediments in dysaerobic and anoxic basins are commonly banded, which is indicative of fluctuation conditions of deposition. These bands have often shown to be varves (e.g. CALVERT. 1966; BRULAnD, 1974; DEMASTER. 1979). Adoption of a steady-state approach must be justified in this situation. LASAGA and HOLLAND(1976) have investigated the effects of a transient (sinusoidal) input of a reactive solid on porewater profiles of a product species. Specifically, they examined organic matter and its decomposition products, but there is little reason to doubt that their findings are equally applicable to the opal-silica problem. LASAC;Aand Hot.l.anl9 (1976) established that the transient profiles of product species will resemble those at steady state with a mean-value input if the frequency of modulation, a, satisfies the inequality a > 40~"

(5)

The right-hand sidc of equation (5) is at a maximum for the greatest burial velocity of interest. As employed below, w is no greater than I cm y- t. With this value, the right-hand side of equation (5) is of the order of 10 -2 y-I. Yearly laminated (varved) sediments alternate with a frequency of 1 y-l; consequently, equation (5) is satisfied and the steady-state description is useful. It would also be useful to ignore porewater advection; however, it is reasonable to ask if this assumption is valid in the more rapidly accumulating sediments dealt with below. Advective transport is a secondary effect in equation (I) if uL

D~

< 1,

(6)

where L is the depth-scale of interest. The left-hand side of equation (6) is the porewater Peclet number (BEAR, 1988; its inverse as given by BERNER, 1980). For all the cases considered in this paper, L -< 100 cm. With maximum u = w = 1 cm y - t and D, on the order of 100 cm z y - i, inequality (6) is true for all calculations in this paper, and porewater advection is negligible.

Modellingearly diagenesisof silica

1547

Assuming the validity of the steady state (at least as a time-averaged condition) and neglecting porosity variations with depth, equations (1)-(3) simplify to o,

dWbdx+ ~ -

b

= 0

(7)

(Cs - C) bpt,= 0

(8)

+ ks( Cs - C)

dwm d.r

- 0.

(9)

Many of the parameters in these conservation equations are not known with any degree of certainty. This applies particularly to the rate constant, ks and the asymptotic concentration, Cs. The value of ks is probably not known to better that an order of magnitude. This paper employs values between 3 and 500 y-t for comparative purposes (BoUDREAU, 1990). The asymptotic dissolved silica concentration. Cs. should attain theoretically the solubility of amorphous opal. i.e. about 1-1.5 mM. This value is seldom, if ever, obtained in surficial sediments, even in siliceous oozes. SCHINKand GUINASSO(1980) attribute this phenomenon to the formation of coatings, perhaps involving iron, manganese or aluminum. The work of LEWIN (1961), [LER(1973) and VAN BENNEKOMet al. (1989) all support this contention in regard to aluminum. As a result, SCmNK and GUtNASSSO(1980) adopt the ad hoc device of setting C~ equal to the observed asymptotic concentration in order to circunlvent this problem. An anonymous reviewer of the paper also suggests its an alternative nlcchanism that the release of dissolved silica is compensated at depth by the rcprecipitation of silica due either to diagenctic reactions with clays (e.g. MACKINand AI.I.I.:R, 1989) or quartz overgrowths. In addition, HESSE (1988), using Dccp Sea Drilling Project porcwatcr data, has shown that interstitial silica concentrations do reach levels approaching amorphous opal saturation during later diagencsis before falling duc to opaI-CT formation. This nlcans that the lower "asymptotic" concentrations seen in the upper metre of sediments arc not simple thermodynamic equilibrium values. Unforttmately, withottt ;t specific mathematical expression for the chemical mechanism that regulates C~, I am forced to continue to use the SCHtNK and GUtNASSO (1980) device. In any event, it will be shown that the model predictions are largely a function of the total saturation deficit, AC = C~ - C(0), and so depend less on the absolute value of C~. The model remains incomplete without the specification of boundary conditions. Formally, the boundary conditions for dissolved silica, opal and clay at the sedimentwater interface (x = 0) are, respectively, C(0) = Cw

(10a)

(1 - ~e).,(0)b(0) = F.

(10b)

(1 - q0w(0)m(0) = F m,

(10c)

where Fb and F m are the prescribed input fluxes of opal and clay. respectively. As equation (7) is a second-order differential equation, a second boundary condition must bc stated. The most natur,'d condition for C is

].548

B.P. BOUDIIEAU

dCJ

= 0.

(II)

All other conservation equations are first order, and no other boundary conditions can be imposed. Boundary conditions are simply statements of the environmental conditions of the system of interest. The range of values assigned to these input parameters should reflect the natural variability and uncertainties. For example, the dissolved silica concentration in the overlying waters, C , is reasonably constant at about 0.1 mM for all situations considered below, while the input of clay, Fm. varies widely, i.e. 0.005-0.5 gcm -z y- i. The rain of opal is similarly variable. The model calculations should account for this fact. The ordinary differential equations (7)-(9), subject to boundary conditions (10a,b,c,) and (11), can be integrated once to produce a single nonlinear ordinary differential equation with one initial condition. Further integration must be done numerically. The details are given in the Appendix. as well as my rationale for seeking an analytical solution. A few germane results are considered here. The critical flux of opal, F~r, is defined as the input rain of opal for which all incoming opal would dissolve and still succeed in saturating the porewaters for a given set of parameter values. Specifically. if Fb > F~r, then opal will accumulate and if Fb
Fml___~+ Pm

FmPb

+ (yq~)'k,D,ACZ

,

(12)

L~ P., /

where AC is the total saturation deficit. The critical flux also can bc interpreted as the maximum possible dissolutivc flux h~r given conditions. For F b > F¢.¢, the potential h~r dissolution is smaller than the actual flux, thus leading to preservation of opal. For Fh < F~,, the potential is greater and all the input dissolves. A second useful set of formulae to bc found in the Appendix is an explicit approximate solution to the nlodel valid when the fraction of the incoming opal flux that dissolves is less than about 50%, i.e. (1 -

q0)w(Qo)b(o0) ~ 0.5.

(13)

F. When equation (13) is true, then the dissolved silica concentration is approximately given by C ( x ) = C, - ACexp

k . b ~]t/2 \ - [k"b(~)/ x/ [ D~J J

(14)

which is a familar form. The concentration of opal that is preserved, b(~). is found as the correct root of the transcendental equation F,,,b( ~ )

The utility of this simple approximation is discussed below.

Modellingearly diagenesisof sii;ca

1549

PROPERTIES OF THE MODEL Without any knowledge of the properties of the model, it becomes difficult to understand the reason(s) why the model may succeed or fail to reproduce porewater or sediment profiles. This section of the paper establishes the sensitivity of the predictions to the model parameters. This examination begins by considering the relationship between the amount of opal preserved at depth in the sediment, b(~), as a function of the total saturation deficit, ~C = Cs - C,, and changes in various parameter values. Figure 1 displays the effect of changing the clay input, Fro, on the amount of opal preserved. The left-hand diagrams in this figure are for a ks value of 30 y-t and the right-hand plots are for ks = 300 y-t. The curves on the plots correspond to opal input values. Fb. needed to produce opal concentrations similar to those seen in nature, i.e. 0--40% opal. Two basic trends are evident in these diagrams. Dissolution becomes less noticeable as Fm increases and also as k~ decreases. Both trends are easily anticipated, but the magnitude of the effect of increasing Fm (at constant Fb/Fm) is nevertheless striking. The relatively horizontal attitude of the % opal curves for the higher Fro-values indicates that dissolution has become unimportant compared to burial in determining the fate of opal over the AC-range of 0-1 mM. This last assertion is verified by the results in Fig. 2. which plots the ratio of the critical opal flux. For as defined by equation (12), to the actual flux. Fb. for cases inclusive of those in Fig. 1. As stated previously, F,r is a measure of the potential for dissolution under a given set of conditions. At F,:JF, = 1, all opal will dissolve, while preservation will occur if this ratio is less than one. Figure 2 shows that the ratio F¢rlF, attains values of onc or greater for the smaller input rains of clay; thus, opal will totally disappear in such slowly accumulating scdimcnts. For larger Fro-values. this ratio becomes considerably less than onc, which is indicative of a great potential for preservation. For (yq~p,,/I~)-"k,,D,,AC'/F~, << 1. thc binomial expansion of equation (12) gives

b'~, ~ (y~p)"p,,,k,,D,, AC z th, F;,,

(16a)

Thus. in this limit, the potential h)r dissolution falls inversely with the rain of inert material. This behaviour can be seen in Fig. 2 in transition from F m = 50--500 mg cm-" y- t Notice also that the curvature of the plotted lines for these large Fm values is consistent with a squared dependence on AC. In the two diagrams of Fig. 2 for the smaller F~ values, the F,:,IFb curves are nearly linear in AC except near the origin. This is consistent with the asymptotic form of equation (12) for the case where (7¢:p,Jph)"k~D~AC'/F~, >> 1, i.e.

F,:r ~ yq~( k~O~) l~ A C

(16b)

which is also independent of Fm. The small divergence of advective flux implied in Figs 1 and 2 suggests that the conditions for the validity of equation (15), i.e. equation (13), may be fulfilled for large Fm values. Figure 3 illustrates that this is the ease. The dashed lines in the four diagrams of Fig. 3 are the solutions given by equation (15). These asymptotic solutions are essentially indistinguishable from the exact solution for the large F m value (Fig. 3b), but deviate appreciably at the smaller rain of clay. Based on clay inputs for non-mixed sediments in

1550

B . P . BOUDa~,U

(o)

(b)

50-

k,=30

50

k,=~O

Fo- 500

F. =500

40"_300

40' 3OO-

Opal

30. Fb ?

30" Fb

Preserved

Jv

%

20"

Jr 20"

~100

100

10"

10" --~ ~lO

0 0.0

. . 0.2

. . 0,4

.

. 0.6

0.8

50"

40"

1.0

-30 -IO 0. . . . . . . 0.0 0.2 0.4

k,:30

50

F.

401

: .SO

, 0.6

. 0.8

i

1.0

k~ = 300

F.=50

t-3o~

~ 3o---.-----~

t

% Opal Preserved 20"

F~ ~

201 ~

-,o-__._

~ , o ~ I

10"

101

0

"

0.0

50

•

-

,

0.2

0.4

'

J

-

•

0.6

-"

•

0

08

'

+ ....

1.0

k,

30 F,,, - 5 II

40

0.0

,

0.2

-

,

-

0.4

•

-

0.6

•

-

0.8

50

1.0

k, - 3(X) F,,, ~ 5

40

%, Opal

i

0

0.0

0.2

0.4 0.6 &C (raM)

0.8

1.0

-'"

0.0

0.2

-

•

0.4 0.6 /,C (raM)

-

,"

-

0.8

Fig. I. The steady-state opal concentration (% total solids) to bc found in non-mixed sediments as predicted by the model equations (4), and (7)-(9). The calculations arc done for two values of the dissolution rate constant, k, = 30 y- J for the diagram in (a) and k, -- 300 y- J for the diagrams in (b). The diagrams correspond to three values of the rain of inert materials, Fm = 5, 50 and 500 mg cm-2 y - I for the top, centre and bottom diagrams, respectively. The curves in thc diagrams arc gcncratcd with thc various inputs of opal, Fb, as indicated by the labels on the curves (in units of mg

cm-" y-').

1.0

Modelfing early diagenesis of silica (a)

ooo~ ,.=~o F. =500

/ /

0.002

1551 (b)

003 ,.=:o

/

/

F. ~0

/

0.02"

F,

Fb

0.001

Fb---~- IO"

/

0.01

o.~

- J "- : ~ - ~ . ~ ~ -

I ~ o.o

02

o,

o~

o.~

08

0.3 k,=30

,.o

o.o

/

o.2

o,

o.~

0.3 k,=~O

08

,o

/

= so

0.2

/

0.2"

F,,

N

F~

~' o

F,,---* 3 o

0.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.0

02

0.4

0.6

08

t,o

!1 _0/ ~

2

4" .

f 0

~'o""-::"~_" "

..........

0.0

0.4 0.6 AC (mMI

0.8

1.0

.

.

2O,

0.2

.

0.0

-

•

0.2

-

.

-

,

0.4 0.6 AC (raM)

, 0.8

Fig. 2. Ratio of the critical opal flux, F~, to the actual opal flux, Fb, for the cases displayed in Fig. I. The critical opal flux is not only the opal input that would cause the porewaters to saturate while preserving no opal. it is also the maximum possible dissolutive flux for the chosen conditions. The dashed line in the lower two diagrams indicates a ratio of one. No opal will be preserved under conditions that have ratios greater than one.

. 1.0

1552

S.P. aot, o=e_,~u

(a) so I/ .......

E,=ct

k, = 3O

Asymptotic

Fm = 5

40 ~

k, = Fm = 5

40

.0 ~

~ ] . O ~ 4 r - - F.

30

~ ~ '

30

~ "~

% Opal

'~'~L0 =

Preserved

20

~ ~

.

0.0

0.2

0.4 ~C

_

0.6 (raM)

~ ......

Exact Asymptoli¢

o ~ __\ x ~.,~ x " ~ .

~..

- - . _

0.5

~,

~X~' 2.0 \\ x

20

-~ -~

"~.

Fb

0.8

-..

9 1.0

0.0

02

0.4 0.6 AC (raM)

0.8

1.0

(b) ,50

6"

5"

"300---,.

3 &30 30

% Opml Preserved

3oO

3 20 2" 10 t ........ 0 " "" 0.0 0.2

Exacl Asymptotic • • 0.4 06 AC (raM)

-

Fb = 30 Fm = ~ 0 • 0.8 t .0

~ Exact ......... Asymptotic 0 . . . . . . • 0.0 0.2 0.4 0.6 &C (raM,)

-

Fb = 300 Fm = 500 • 0.8 1.0

Fig. 3. A comparison o f the exact semi-analytical solution to the silica di,tgcncsis model in non-mixed sediments, equations (A20) o f the A p p e n d i x . with an appropriate solution, i.e. equation ( 15)o which is valid when the silica burial flux is greater than twice the dissolufiv¢ flux. i.e. equation (13). The set o f diagrams marked (a) represent a slowly accumulating (deep-sea type) sediment where F,, = 5 mg ¢m -2 y-J. T w o dissolution rate constants./q, are compared and the curves correspond to various inputs o f opal (in mg crn -2 y - i ). The set o f diagrams marked as (b) is for a rapidly accumulating sediment (F,n = 500 mg cm -2 y - =) and for two inputs o f opal. The curves in (b) correspond to three different k, values.

such basins as those of the California Borderlands (KOIDE et al., 1973; SCHWALBACHand GORSLINE, 1985), the Gulf of California (v^s ANDEL, 1964; DEMASTER, 1979), and the Peru-Chile Margin (DEMASTER, 1979; REIME~S, 1982), it is anticipated that the much simpler asymptotic solution should provide an accurate description of opal diagenesis in these environments. The role of mixing in promoting opal preservation can be appreciated by comparing the results in Fig. 1 to those calculated with a model in which mixing is present. The difference

Modellingearly diagenesisof silica

1553

is most pronounced if mixing is infinitely fast, which is a reasonably accurate approximation to the SCmNK and GUINASSO(1980) model (BouDtEAU, 1990). These results are displayed in Fig. 4. In comparing Figs 1 and 4 it can be seen that the effects of mixing are observable only for the smallest clay input rains, Fro- Under these conditions, FcrlFb ~ 1 and, with no mixing, all opal dissolves with increasing AC. As explained by SCHINK et al., (1975), mixing removes opal from the zone of dissolution and causes its preservation, i.e. the long tailing of the curves in Fig. 4. Mixing has little infuence at the higher F m values because transport by advection is already sufficiently large to overcome dissolution. One fundamental difference between the mixed and non-mixed situations is the amount of opal at the sediment-water interface. Without mixing, b(0) is independent of the extent of dissolution. This can be seen by combining equations (4) and (10b.c) to obtain b(0) =

Fb

(17)

This relationship states that the opal concentration at the interface is always equal to that in the absence of dissolution. The amount of opal preserved, b ( ~ ) . will fall with 5C, and this can lead to sharp depth gradients in the opal concentration. With diffusion-like mixing, however, similar manipulations of the appropriate boundary conditions (ScmNK and GUINASSO, 1980; BOU~REAU, 1990) lead to Fb + (1 - q0D,, ~----bL~,, b(O) =

(18) P,, /

Because the gradient in opal, dbldx, at x = 0 depends on the intensity of dissolution, b(0) will drop as the integrated amount of reaction increases. The extreme case occurs as Du--* ~, and b(0)---* b(~). The invariant b(0) property of non-mixed systems can lead not only to sharp gradients in opal concentration, but also in the burial velocity. This is illustrated in Figs 5 and 6 for clay fluxes of 5 and 500 mg cm-2 y-I, respectively. For the low clay flux, the surficial opal concentration and burial velocity are no less than twice that at great depth, while for the larger clay flux, the gradient is no greater than a 10% variation.

APPLICATION OF THE MODEL AND DISCUSSION Model simulations of opal diagenesis in three non-mixed sediments are discussed beh)w in order to illustrate both the strengths and shortcomings of the model and to display the insights the model can provide regarding diagenesis at these locations.

C a r m e n Basin

The Gulf of California contains many small basins that accumulate large amounts of opal (CALVERT, 1964, 1966). In particular, the sediments of the Carmen Basin have been

1554

B.P. BOUD~.AU

(a)

(b)

GO " M i x e d

F~

•o

% Opal

50

k s ~ 30

k, = 300

Mixed

F,= =500

~

~

,o-

300

30'

F,,,, = ~ 0

-300

30"

Preserved 20'

20" |00

100--

10

-

10" 3O

30

I0 0

•

O0

-

0.2

•

I0 . . . . . .

0.4

0.6

0

08

50" Mixed

1.0

k s = 30

, 0.0

50

-

•

0.2

0,4

0.6

0,8

Mixed

F , = .50

40 ~3O

40 ~ 3 0 ~ _ ~

30"

30

? Fb

20 ~

Jr 1

20 0

~

~

I0

I0-.~

I0

~1-0 0.0

• 0.2

.

. . 0.4

.

. 0.6

0.8

50" Mixed

1.0

k , - 30

- - I ~ 0 . . . 0.0 0.2

50

.

. 04

.

. . 0.6

. 0.8

Mixed

F.-~

40"

40

pr%esO/;:ed

3.0 ~

Fb

,Ot.o':,° 0.0

1.0

k, m 300

F.-5

o

1.0

k I = 300

Fm '= 50

% Opal Preserved

-

,o

i

o 0.2

0.4 &C

0.6 (raM)

0.8

1.0

0,0

.

.

0.2

.

.

.

0.4

.

.

0.6

.

.

0.8

&C (raM)

Fig. 4. P e r c e n t a g e of o p a l to bc found in b i o l u r b a t e d s e d i m e n t s as p r e d i c t e d by a s s u m i n g infinitely fast m i x i n g of the solids with all o t h e r c o n d i t i o n s the s a m e as those in Fig. I.

1.0

Modellingearlydiagenesisof silica

0

0

10'

Depth (cm)

ql Opal 20 40 60 . ~ _ _. ~ _, ~ , , I

80 - "

~

l < i - - b ( ,--)

30'

I

40

Depth

(cm)

% Opal 10 20 30 40 50 60 . . . . . . . . . . .

k=:300

70 80 '

20"

4oll

F,,, = 5

Burial Velecily (cm/yr) 0.000 0.002 0.004 0.006 0.008 0.010 0 . • - • , = , !_

20 I0

k,= :~0

40

to

J

003

(---all b(--)

30

lO"

30

Burial Veleci¢7 (cndyr) 0.01 0.02 * ...~ '

201

F,,, =5

0

0.00 Ot

10t

20"

0

100

1555

30" k,-30

k,,, ]0

F.,,5

40- F.:-5

Fig. 5. Depth profiles of opal concentration (as percent of total solids) and burial velocity (cm y-I) in rclativcly slowly accumulating non-mixed sediment. Fm - 5 mg cm -z y-I and for two/q values, i.e. 30 y-= in the Icft-hand diagram and 300 y-i in the diagram on the right. Each curve corrcsponds to a different final opal content indicated by the label (in units of mg cm -z y-I).

found to be varved (DEMASTER,1979) and to accumulate at a velocity of about 0.2 cm y-t (CALVERT, 1966; DEM^STER, 1979). A porewater dissolved silica profile from the sediments in this basin has been obtained by J. GIESKES (personal communication), and the data are shown in Fig. 7. No corresponding opal profile is available, but the final opal content should be between 10 and 30%, based on the opal content maps of CALVERT(1966) and the data of DF-MASTER(1979). Figure 7 also displays the model curves for two possible final opal concentrations, i.e. 12% in (a) and 30% in (b). The various curves in these diagrams correspond to different k, values. Overall, the data is well represented by the model by choosing a ka value between 300 and 400 y- i. Such close coincidence is visually impressive, but hardly diagnostic of the worth of the model. The rate constant and the final opal content act as free parameters, which generate the close fit. Nevertheless, the model is successful in this application and can be used to calculate other important information such as the interfaeial flux.

1556

B.P. BOUOIEAU

San P e d r o Basin

As part of their bottom lander experiments, BErEt,SON et al. (1987) have reported porewater and solid silica data from the San Pedro Basin of the California Borderlands. These sediments are unbioturbated and accumulate at a speed of about 0.02 cm y-t (BRULAND et al., 1981; SCHWALBACHand GORSLXNE, 1985). The sediments contain 3.5 +_2% opal by weight. BERELSONet al. (1987) have also determined the dissolved silica flux out of the sediments by direct measurement with the lander, by diagenetic modelling and by modelling the water column data. They state values for this flux of 0.026 ___0.004 mmol cm-Z y-t (0.7 mmol m -2 day-t). The data set provided by BEREt,SON et al. (1987) is particularly valuable in that the measured flux constitutes an independent check on the model predictions that also eliminates some of the arbitrariness in setting the value of k,. The visual best-fits of the model to the data are shown in Fig. 8a. Figure 8b illustrates the relationship between the model-predicted flux of dissolved silica at the sediment-water

10

% Opal 20

Ilurial V¢lecil7 (cm/yr) 30

° i i"

40

S

10

~ ~

0

30"

% Opal 20 '

,

k= = 300 F.. -=~oo

40 J

Burial Vclo¢ily (¢ndyr) 30

40

10

40

t.

b(--)

k= = 300 I~ = 5o0

10

30'

0.110

o.--, . . . .

I

40

20

0.105

20" I0

Depth (cm)

0.100

10"

Depth (cm) 20"

30"

0.095

0.0975 0

10

I

O,I (300

O.1025

O.1050

I b(--)

20 ,

I0

30 ~

1

I

(--)

ks=30 F.=~O

30" 40"

k~=30 F.=~O

Fig. 6. Depth profiles of opal concentration (as percent of total solids) and burial velocity (cm y - t ) for a non-mixed sediment accumulating clay 100 times more rapidly than that in Fig. 5, i.¢. F m = 500 mg cm -z y - I , and for the same two k, values, i.¢. 30 y-t in the Icft-hand diagram and 300 y-n in the diagram on the right. Each curve corresponds to a different final opal content indicated by the label (in units of mg cm-Z y - t).

1557

Modelling early diagencs/s of silica

Silica 0.0 0

-

'

0.2 '

.....

(raM) 0.4 '---'m

Sili¢ll 0.6 ~

0.8 '

0.0 0

_- -

10

10"

:,'0 ~

20"

Depth (cm)

' 30"

C a r m e n Basin

(b(~)- IZ%)

30"

ks - 20(l/yr 40-

- ..... •

50

ks-40~/yt ks - IOO(l/yr dam

i

(raM) 0.4 n

~

0.6 i

0.8

Carmen

Basin (b(--) = 30%)

......... 40"

0.2 I

~ ..... ' •

ks - ~O/y¢ ks=3(Xl/la, ks - ~ q l / y r alum

50

Fig. 7. Dissolved silica data (from J. GiESK~, personal communication) obtained from the varved sediments of the Carmen Basin of the Gulf of California and some model-predicted curves based on the assumption of 12% final opal content in (a) and 30% in (b). The model curves are generated with varying values of the dissolution rate constant./q as indicated in the legend of each diagram.

interface as a function ofk,, for three AC values. For the San Pedro Basin, AC ~ 0.4 mM if C, is taken to be about 0.5 raM. The model fits to the silica data suggest a ks value of 150-200 y-i. This results in a predicted dissolved silica flux at the sediment-water interface of between 0.02 and 0.03 mmol cm -2 y - i or a figure within the range stated by BERELSONet al. (1987). The agrccmcnt is highly encouraging. Peru Margin

REtMERS (1982) conducted an extensive geochemical and gcotechnical analysis of the sediments off central Peru. She obtained profiles of both dissolved silica and opal at two stations, one of which is shown in Fig. 9. The opal distribution is distinguished by a sharp drop at a depth of about 12 cm, i.e. ~3.3 to ~1%. D~M^sTEa (1979) has determined that this sediment is accumulating at about 0.13 cm y-l. The dissolved silica increases with depth, but with a gentler gradient than in the Carmen Basin. (There is indication of a subsurface maximum in the dissolved silica concentration profiles. This may be due to mineral precipitation reactions; alternatively, some minor irrigation may be present. For the purposes of this paper, the inflection will be ignored.) Figure 9 also contains model predictions for both the dissolved silica and solid opal generated with values of ks between 10 and 90 y-I and b(oo) set to 1%. The best correspondence between model and the dissolved silica data demands a k~ value somewhere in the range of 30-60 y - t , but such a fit is only marginally acceptable. More to the point, the model fails to produce passable simulations of the corresponding opal profile. While an adequate fit to the solid data could be achieved by increasing the value of k,, this would lead to a poor prediction of the dissolved data. This disagreement is best resolved by admitting to the existence of an essential nonsteady state condition due either to an increase in the opal flux about 100 eP or the recent rapid deposition of sediment richer in opal, e.g. by a slide. Although the model has not proven to be formally successful in this case, it has provided valuable confirmation on the

1558

B . P . Bouotr~u

state of the system by confirming the change in opal concentration with depth is not likely diagenetic. Any diagenetic model is a highly idealized abstraction of nature. Any fit to data will be at best imperfect, and a failure to accurately reproduce a particular data set should come as no surprise. Judging the value of a model solely on its successful conformity with a data set is dubious. Specifically, the emergence of an unfavourable comparison between model and data, i.e. a negative result, is sometimes taken as an excuse to dismiss the model. Thus, while the good fits to the Carmen and San Pedro Basin data may be pleasing, the lack of agreement for the Peru Margin data is equally edifying. EXTENSION AND MODIFICATIONS TO THE MODEL The inherently incomplete character of a diagenetic model encourages its modification and extension. The present model certainly contains some approximations and assump(raM)

Silica 0.0

0.1

0.2

0.3

(O) 0 . . . . . . . 10

Depth

0.4

'

-"

0.5

0.6

200

t

"

"

00

20

•

•

(cm)

•

30

I

•

40

San Pedro 50

San Pedro Basin

(b) o.10

J

0.O8

0.45

0.02

0.00

-

0

,

200

-

,

-

400

,

600

-

,

800

-

1000

k. (yr "l)

Fig. 8. Dissolved silica (a) from the non-mixed sediments of the San Pedro Basin of the California Borderlands, as reported by BERELSONel oi. (1987). The curves are model-generated predictions, equation (A16), with b ( ~ ) : 3.5%, w ( ~ ) = 0.02 cm y-S, and the ks values indicated by the labels in the diagram. (b) illustrates the dependence of the model-predicted flux of dissolved silica out of the sediment as a function of the rate constant for dissolution,/%, and three values of the overall saturation deficit, A C = C, - Cw. The in situ AC in the San Pedro Basin is near 0.4 raM. The dashed line is the value of the mean flux as measured by BEaELSOt'~etal. (1987).

1559

Modelling early diagenesis of silica (a)

stiles

(raM)

0.0 0.1 0.2 0.3 0 . . . . . . . . . .

0.4

(b) 0.5 0.6 ' . . . .

0.7

%

0

1

2

0

3

4 •

I0 10

•

10

•

9O 20'

20

Depth tom)

• 30'

40

30

Peru

•

40 •

50

Peru

50

Fig. 9. Dissolved silica data (a) and solid opal data (b) collected by REIMER$(1982) from the non-mixed sediments of the central Peru Margin. The curves are model predictions for k., values indicated by the labels (k, in units of y-~).

tions that should be removed. An obvious change is the incorporation of a porosity variation with depth, as porosity changes can lead to modifications in the expected shape of profiles (e.g. JORGENSEN, 1978). This will necessitate the use of a numerical method of solution. The manner in which adsorption could be incorporated can be based on the dcvelopment found in BERNER(1980) and need not be repeated here. I now discuss more extensively three other possible changes to the present modch the first deals with the input of multiple opal types, the second with the inability to reach the theoretical solubility at shallow depths, and the third rehttes to recrystallization during later diagenesis.

hzput of muhiple opal types SCmNK and GUINASSO(1981)) have advanced the idea that opal entcring the sediment is not of one reactive type, but of many. ANDREWS and FIARGRAVE(1984) present strong evidence for the existence of at least two types: one of high reactivity that dissolves near the sediment-water interface and another which dissolves more slowly. [In addition, the reactivity of the opal may change with time (depth) due to the "aging" of the material, as discussed in the next section.] The model of this paper can be amended to account for n-types of opal. This requires stating a conservation equation for each of the "n" reactive types (as modified from SCHINK and GUINASSO, 1980), d(1 - q))wbi + yqgki(Ci _ C) bi = 0 dx Pi

(19)

along with new boundary conditions F; = (l - (p)w(0)bi(0).

(20)

where Ci and Pi are the solubility and density of the ith opal type. respectively. There is no a priori reason to believe that the opal types entering the sediment have the same solubility

1560

B.P. BouDRr~u

or that they have the same density (although the latter assumption may introduce little error even if it is not true). The dissolved silica equation must be modified to reflect the multiple sources

d [~osd-~-~] + q: ~ k,(C, - c) b--2= O. i~l

(21)

Pi

Because some o f the more reactive opal types may have quite large rate constants of dissolution, the characteristic depth for their dissolution may approach the thickness of the diffusive sublayer of the benthic boundary layer (BouDrEAU and GUINASSO, 1982; ANDREWS and HARGRAVE.1984). This would require alteration of the dissolved boundary condition at x = 0 to D~ dC =

-

c.),

(~.~) . .

where C(0) is the dissolved silica concentration at the sediment-water interface, Cw is now the concentration in the bottom water immediately above the boundary layer, and fl is the mass transfer coefficient, also called the piston velocity (BouDREAU and GUINASSO, 1982). Overall. it is unlikely that this set of equations will be assessible through analytical methods.

The asymptotic concentration The fact that the theoretical solubility of amorphous opal is rarely attained at shallow depths during early diagcnesis suggests that the apparent solubility is altered by some other chemical process (e.g. precipitation of another mineral) or that the kinetics of dissolution become inhibited, it is beyond the scope of this paper to identify which process is responsible, but the conscqt, enccs of these processes to modelling can nevertheless be discussed. If the solubility is reduced by an inhibition process with depth. C~ or Ci must be replaced by an appropriate depth-dependent function that preferably is calculated from some mechanistic model. Likewise, if the rate of reaction is inhibited, k~ or k i must be written as a function of depth. In either case, prospects for an analytical solution are dim and numerical methods must be considered. Note that the asymptotic behaviour of the system with a depth-dependent solubility is different from that with kinetic inhibition. In the solubility case, the concentration at depth is equal to a lower prescribed value, while in the kinetically inhibited system, the asymptotic concentration is still the solubility of amorphous silica. In the later case, the slowed pace of dissolution would allow attainment of the solubility of amorphous silica only at great depth. Distinguishing between these limiting cases may prove to be difficult, but the data in HESSe (1988) appears to favour the second possibility. If the early diagenetic precipitation of another silicon-bearing mineral of lower solubility (e.g. quartz overgrowths or a clay) is the culprit in preventing porewaters from attaining opal saturation, then this can be treated in a manner similar to that discussed in the next section.

Modellingearly diagenesis of silica

1561

Later diagenesis Whereas our interest in this paper was directed at early diagenesis, mineralogical and chemical changes in the silica-opal system are documented to continue well below the top metre of sediments (CALVERT. 197 i; HESSE, 1988). Specifically, a sequence of mineralogical transformation can be observed with increasing depth (CALVERT, 1971; KASrNE~ et al., 1977; WILUAMSet al., 1985; WtLUAMS and CRERAR, 1985; HESSE, 1988) opaI-A ---, opal-CT---, cryptocrystalline quartz,

(23)

where opal-A is the biogenic amorphous opal found in tests raining to the sea floor and opaI-CT is a less soluble disordered cristobalite-tridymite phase. These transformations are said to occur via dissolution-reprecipitation reactions (KASTr~ERet al., 1977), rather than solid phase changes. The porewater dissolved silica will adjust itself to the changing saturation concentrations associated with the solid phase transformations. Although the present model was devised for non-mixed sediments, it also applies to diagenesis below the zone of mixing in bioturbated sediments; therefore, it can be used to describe later diagenesis in all types of sediments. The incorporation of late diagenetic reactions into the model is formally facile. If there is only one type of opaI-A and the precipitation of the opaI-CT is a first order process (abundant nuclei), then the conservation equation for dissolved silica is (at constant porosity) D~ ~ + kA(C A t..1.1f-

C)

b.AA_ kc-r(Ct--r - C) = 0 Pl,

(24)

and for the opaI-CT dwbc-r

yrFkc'r ( C c l -

C) = O,

(25)

where subscripts A and C'7I"indicate the mineral. Meanwhile, tile opai-A concentration is governed by equation (8). if there is more than one reactive type of opal, then equation (24) must be modilicd in the image of equation (21). Method of solution While no attempt is made in this paper to solve the models given by equations (19)-(25), I would like to recommend a numerical technique likely to be successful and useable by anyone with a computer with a FORTRAN compiler. ASCHERet al. (1981) and BADERand ASCtlER (1987) have developed a software code for the solution of two-point boundary value problems by the collocation method. The code is powerful and relatively easy to adapt to the problems at hand. (Partial derivatives needed for the so-called Jacobian can be calculated analytically, which requires some basic calculus skills, or numerically following the suggestion in ASCHERet al., 1981.) Copies of the code are readily available from the creators (see the appendix in ASCHERetal., 1988) or failing that, from this author. I have successfully employed this code to solve the SCHINK and GUINASSO (1980) model. Acknowledgements--This research was supported by a Natural Sciencesand Engineering ResearchCouncil of Canada (NSERC) Grant (URF0043963). Dr Joris Gieskcs and Dr William Bcrclson kindly made their data available to me. Comments from the anonymousreviewerssubstantially improved the manuscript.

1562

B.P. Bovaa~u REFERENCES

ANDJEWS D. and B. T. l-haGR^vt (1984) Close interval sampling of interstitial silicate and porosity in marine sediments. Geochimica et Cosmochimica Acta. 4 , 711-7,--r2. A s c n t a U., J. Cn~s'~NstN and R. D. RUSSELL (1981) Collocation software for boundary-value ODEs. Association for Computing Machinery (ACM) Transactions on Mathematical Software, 7.2(D- ~222. AscnEa U. M., R. M. M. MArmEU and R. D, RUSSELL(1988) Numericalsolution of boundary value problems for ordinary differential equations. Prentice Hall. Englewood Cliffs. New Jersey. 595 pp. BADER G. and U. AscxEx (1987) A new basis implementation for mixed order boundary value ODE solver. Society for Industrial and Applied Mathematics Journal of Scientific and Statistical Computation, 8.483500. B~,~R J. (1988) Dynamics offluids in porous media. Dover Publications. New York. 764 pp. BEgE~ON W. M.. D. E. HAMMONDand K. S. JOHNSON(1987) Benthic fluxes and the cycling of biogenic silica and carbon in two southern California borderland basins. Geochimica et Cosmochimica Acta. 51, 1345-1363. BERNEa R. A. (1974) Kinetic models for the early diagenesis of nitrogen, sulfur, phosphorus, and silicon in anoxic marine sediments. In: Thesea. Vol. 5. E. D. GOLDaERG. editor. John Wiley. New York. pp. 427--450. BERNER R. A. (1980) Early diagenesis: a theoretical approach, Princeton University Press. Princeton, New Jersey. 241 pp. BoucaEa J. M. (1984) Silica dissolution and reaction kinetics in southern California Borderland sediments, M.S. thesis. University of Southern California, Los Angeles. 149 pp. BOUDREAU B. P. (1985) Diagenetic models of biological processes in aquatic sediments. Doctoral dissertation, Yale University. New Haven. Connecticut. 526 pp. BOUDREAU B. P. (1990) Asymptotic forms and solutions of the model for silica-opal diagenesis in bioturbated sediments. Journal of Geophysical Research. 95. 7367-7379. BOUDREAU B. P. and N, L. GtrIN^sso. Jr (1982) The influence of a diffusive sublayer on accretion, dissolution, and diagcnesis at the sea floor. In: The dynamic environment of the ocean floor. K. A. FA,~NL~C,and F. T. MAN,ElM. editors. Lexington Books. Lexington, Massachusetts. pp. 115-t45. BauLANt~ K. W. (1974) Lcad-21(J disequilibria and geochronologies in the marine environment. Ph.D. thesis. Scripps Institution of Oceanography. University of California. San Diego. 106 pp. BRUI.AND K. W.. R, P. FRANKS. W. M. LANDING :lHd A. SOUTAR (198l) Southern California intcrbasm trap calibratkm. Earth and Planetary Science Letters. 53, 40(b-4(J8. BraNt" G. D. and A. C. | hnDmaaSH (1987) Stiff ODE solvers: :t review of current and coming attractions. Journal of Computational l'hysics, 70, 1-62. C^l, VEar S. E. (1964) Factors affecting distribution of hminatcd diatomaceous sediments in the Gulf of C;tlifnrnia. In: Marine geology of the Gulf of Califi,rnia, T. t l. vAN AIqDEL and G, G. SmJa. Jr, editors, American Association of Petroleum Geologists, Memoir 3, pp. 311-330. CALVEar S. E. (1966) Accumulatkm of diatomaceous silica in the sediments of the Gulf of California. Geological Society of America Bulletin. 77. 569-596. CAt.YEar S. E. ( 1971) Nature of silica phases in deep sea cherts of the North Atlantic. Nature, 234. 133-134. DEMAsrER D. J. (1979) The marine budgets of silica and 3:Si. Ph. D. thesis. Yale University. New Haven, 308 pp. DEM^sTEa D. J. (1981) The supply and accumulation of silica in the marine environment. Geochimica et Cosmochimica Acta, 45. i 71 5-1732. EREz J., K. TAKAI-IASHIand S. Homo (1982) In-situ dissolution experiment of Radiolaria in the central North Pacific Ocean. Earth and Planetary Science Letters, 59,245--254. FJADEtRO M. E. and G. VEaONtS (1977) On weighted-mean schemes for finite-difference approximation to the advection--diffusion equation. Tellus. 29,512-522. GEAR C. W. (1969) The automatic integration of stiff ordinary differential equations. In: Information processing 68. A. J. H. MORaELL. editor, North-Holland, Amsterdam, pp. 187-193. GLASS J. and W. Root (1982) A higher order numerical scheme for scalar transport. Computer Methods in Applied Mechanics and Engineering, 3 I, 337-358. GRADSHTEYN 1. S. and I. M. RYzlIIg (1980) Table of integrals, series, and products. Academic Press. Orlando, Florida. 1160 pp. H.v.sst: R. (1988) Diagenesis # 13. Origin of chert: diagenesis of biogenic siliceous sediments. Geoscience Canada, i$. 171-192. HtNDMARSH A. C. (1972) G E A R : ordinary differential equation system solver. Technical report UCID-30001, revision 2, Lawrence Livermore Laboratory, University of California. Livermore. California.

Modelling earlydiagenesisof silica

1563

Huav D. C. (1972) Factors affeaing the solution rate of biogcnic opal in seawater. Earth and P/an~ary Science Letters. 15.411-417. HuaD D. C. and BmDWHLSTr~.L (1983) On producing a more general model for biogenic silica dissolution.

American Journal of Science. 283. 1-28. [I.ER R. K. (1973) Effect of adsorbed alumina on the solubility of amorphous silica in water. Journal of Colloidal

and Interracial Science, 43. 399-.408. J~aN M. K., S. R. K. lWNG^R and G. S. SuBam~tANVaM(1984) Variable mesh methods for the numerical solution of two-point singular perturbation problems. Computer Methods in Applied Mechanics and Engineering. 42.273-286. JORGEr~SEr4B. B. (1978) A comparison of methods for the quantification of bacterial sulfate reduction in coastal marine sediments. Geomicrobiology Journal. I. 29--47. KAMAtANI A. and J. P. RxLeY (1979) Rate of dissolution o[ diatom silica walls in seawater. Marine Biology. 55. 29--35. KAMATA~I A.. J. P. RtLEV and G. SKIR~OW(1980) The dissolution of opaline silica of diatom tests in sea water.

Journal of the Oceanographical Socie~ of Japan, 36.201-208. KASTNEr M.. J. B. KEErqE and J. M. Gtr.SKES (1977) Diagenesis of siliceous oozes----I. Chemical controls on the rate of opaI-A to opal-CT transformation--an experimental study. Geochimica et Cosmochimica Acta. 41. 1041-1059. KErn R. S. (1982) Dissolution of calcite in the deep-sea: theoretical prediction for the case of uniform size particles settling into a well-mixed sediment. American Journal of Science, 282. 193-236. KOIDE M.. K. W. BRULANDand E. D. GOCOBERO(1973) Th-228/Th-232 and Pb-210 geochrono[ogies in marine and lake sediments. Geochimica et Cosmochimica Act^. 37, [ 171-1187. LAS^G^ A. and H. D. HOLLAND (1976) Mathematical aspects of non-steady-state diagenesis. Geochimica et

Cosmochimica Act^. 40.257-266. Lawsos D. S.. D. C. HURD and H. S. PANK~AtZ (1978) Silica dissolution rates of decomposing phytoplankton assemblages at various temperatures. American JournalofScience, 278, 1373-[393. LeDFoRo-i-lorrm^s P. A.. D. J.. DEM^stE~ and C. A. N,rrRou~a (1986) Biogcnic-silica accumulation in the Ross Sca and the importancc of Antarctic continental-shelf dcposits in Ihc marine silica budgct. Geochimica et Co.smochimica Act^. 50.2(~--2110. LEWIS J. C. ( 1961 ) The dissolution of silica from diatom walls. Geochimica et Cosmochimica Act^. 21. 182-198. MacKtn J. E. and R. C. AI.LeR (1989) The ncarshorc marinc and cstuarinc chcnfistry of dissolved aluminum and rapid suthigcnic mincral prccipilathm. CRC Critical Reviews in Aquatic Sciences, I. 537-554. MuKPIIv G. M. (19(~}) Ordinary differential equations and their soh~tions. D. Van Nostran. Princeton, 449 pp. Rl-,M,~rs C. E. (1c}82) Organic matter in anoxic sediments off central Peru: rclations of porosity, microbizd dccomp~silion and dcform.'tth~n propcrtics. Marine Geology, 46. 175--[97. ScmnJ< D. R. and N. L. (;uis^sso, Jr (1977) Effects of bioturbation on sedhncnt-scawatcr interaction. Marine Geology, 23, 133-154. SCmnK D. R. and N. L. GUtnASSO.Jr (1980) Proccsses affecting silica at thc abyssal sediment-water interface. [n: Biogt~ochimie de la matidre organique d l'interfiwe eau-s~diment marin, Colloqucs [ ntcrnationaux du Centre National de Recherche Scicntifiquc no. 293, pp. 81-92. SCmnK D. R., N. L. GU,NASSO, Jr and K. A. FANNXNG([975) Processes affecting the concentration of silica at the sediment-water interface of thc Atlantic Ocean. Journal of Geophy.s'ical Research. 80, 3013-3031. SCHWAt.aaCH J. R. and D. S. GO~SL,SE (1985) Holoccne sediment budgets for the basins of the California contincnlal borderland. Journal of Sedimentary Petrology, 55,829--842. S~l^Metr41~L. F. and C. W. G r . ~ (1978) A user's view of solving stiff ordinary differential equations. Society for Industrial and Applied Mathematics Review, 21, 1-17. SHvv W. (1985) A study of finite difference approximations to steady-state, convection-dominated flow problems. Journal of Computational Physics. 57.415-438. VAN Asoel. T. H. (1964) Rcccnt marine sediments of the Gulf of California. In: Marine geology of the Gulf of California, T. H. VAN ASDEL and G. G. S,~o~, Jr, editors, American Association of Petroleum Geologists, Memoirs 3. pp. 216--310. VaN BeNne~oM A. J.. J. H. F. J^,'qS~N. S. J. VaN DEr GAAsr, J. M. VAN IP£~EN and J. PtErE~s (1989) Aluminiumrich opal: an intcrmcdialc in the preservation of biogcnic silica in the Zairc (Congo) deep-sea fan. Deep-Sea Re.search. 36, 173-1¢~). W,LLIAMS L. A. and D. A. C ~ r ^ ~ (1985) Silica diagcncsis, ft. General mechanisms. Journal of Sedimentary Petrology, 55, 312-321

B . P . BOUO~EAU

1564

WiLt.tAMS L. A., G., A. P.~Ic.sand D. A. Clttl,,Lit(1985) Silica diagenesis, [. Solubility controls. Journal of Sedirnenta~ Petrology, 55. 301-311. WoLL.ssr R. (1974) The silica problem, In: The sea, Vol. 5, E. D. GotosEm¢, editor. John Wiley. New York. pp, 359-392_. WoNt G. T. F. and C. E. GROSCH (1978) A mathematical model for the distribution of dissolved silicon in interstitial waters--an analytical approach. Journal of Marine Research. 36. 79..3-750.

APPENDIX This appendix contains the solution, in the form of a first integral, of equations (4). (7), (8) and (9) when subject to boundary conditions (10a)-(I 1). In searching for a solution, it is disturbing to note that equations (7)-(9) contain the four dependent variables in nonlinear combinations. This must be resolved and the number reduced if there is to be any hope of an analytical solution. First note that the product of variables w and b can be considered as a single new variable, wb(x). Then b(.r) is written as a function of this new variable by combining equations (4). (9) and ([0c). wb

b=(i

+"hi I - q:)p., ~, F~

(All

Substitution of equation (A I ) into equations (7) and (8) produces d:C D, ~

~ - ~(1'

+ -

~)Pm

~

=0

(A2)

= 0.

(A3)

:

and

dwh

,:.k,, (C', - C)

+

Although these cqualhms may look more complicated titan thosc they replace, thcy contain only two dependent wkriablcs, i.c.C and wb. [t now rcmains to express one of these in terms of thc other and dcrivc a sing[c equation in a single dcpcndcnt variable. To do this.multiply cquation (A2) by yq,/(I - 9") :rod subtract cqu:~thm (A3) from the result,

IY-..~=D ~ C _ dwb =0.

( - ~j

"

"

(A4)

dx

Integrate cquation (A4) from x to ~ ,

wb(x) =

1/¢ D, dC+ wb(oo) (l - ~,)

(A5)

d.r

where wb(~) is the value of this variable as x ~ ~ and is an unknown value (cigenvalue) that must be calculated as part of the model. Substitution of equation (A5) into equation (A2) produces a single equation containing only the unknown C,

"

~b

\(|--(10) d.x

(|--~)Pm

p t , ~ ( l - g Q d x ewb(oo)

=0

(A6)

which is formidable in appearance, but manageable in reality. To facilitate the intcgration of equation (A6). a new variable is introduced, O = C~- C

(A7a)

Modelling early diagenesis of silica

1565

dO dC . . . . . dz dz

(ATb)

! also introduce the Lagrangian notation for derivatives where one prime indicates a first derivative with respect to x and two primes a second derivative. Equation (A6) then reads after minor algebra

a~O"

+

wb(~)O"

w b ( ~ ) - a:O'

a t O - - 0.

a:O'O" wb( ~ ) - a,O'

w b ( = ) - azO'

(AS)

where a t • kJD~. a 2 - 7 ¢ D , / ( 1 - q'). and a3 - Fa~oJ[(l - qr)p~J. It is necessary to introduce yet another new variable: let • ' = O'.

(A9)

This may look trivial, but it is the key to the solution (MURPHY. 1960). From equation (Ag) O " = W' = ~, ~V.V. dO

(Ai0)

So that equation (AS) reads. a~'d~F w b ( ~ ) - a2O/

+

wb( ~ )Wd~

a:~-'dqt

w b ( ~ ) - ald/

w b ( ~ ) - a2O/

- atOdO.

(All)

The three terms on the left-hand side of equation (AI 1) can be expanded through the use of partial fractions. I

q'

wb( ~ )

W b ( ~ ) - a:q'

a:(wb(~o) - a:q')

(AI2)

--

a:

and

q': =

-

(wh(~)): ,,:,

=

) -

wh( = ) ,,,

-

,,"

'~ -

(Al3)

,7..

Substituting equations (A 12) and (A 13) into equation ( A l l ) and cancelling terms gencr:ttcs ~F -

. dqt + a,(wb(o=).

- a:~') = atlld()

(A 14a)

or with the dciinition of the derivative of a logarithm

¢t2I

Integrating equation (AI4h) 2

- U3~wh(=)Intwb(o0)

- u,'l')

= a,TO: _ ",

(AIS)

=~-(C,-C)'-c,.

(Alfi)

a 2]

or, reverting to the original variables, 2~dx]

Z'~-~

wb(°~)lnwb(=)+"2

where ct is an integration constant that can take on various values depending on the circumstances. Specilically. it there is a sufficient opal flux to saturate the porewaters and accumulate opal. then Cl = ~ w b ( ~ )

In(wb(==)).

If all the incoming opal dissolves without saturating the porewaters, then

(A 17a)

1566

B.P.

Bouoaeau

el = 2 ( £ , - C.)-"

(AlTb)

where C,, is the non-saturated asymptotic value of C. Finally. if the opal flux is equal to the critical flux. For. then CI =0. To establish which of these three possibilities will prevail for a given set of input and parameter values, the critical flux must be calculated. At the critical value of the flux. the diffusive flux out of the sediments exactly

balances the input flux. F~, =

rC,D,~--~r,.

(A18)

Substitution of equation (A 18) into equation (A 16) generates a quadratic for F~,. al

1 2~r:.,+

~ a:y~D,

F~:,_ ~-~L(C - C~) 2 = 0

(AI9)

"

which has as its solution equation (12) of the text. Equation (A 16) that governs the dissolved silica is a highly nonlinear implicit first order ordinary differential equation. A general analytical solution for equation (A 16) appears impossible. This does not mean that all the above work has been for nothing. Equation (A 16) and related results are far more useful than the original set of equations, and i detail the reasons below. First. equation (A t6) with equation (10a) constitutes an intial value problem, which is far easier to integrate numerically than the original two-point lyoundary-value problem. This is particularly true with the present equations as they are numerically "stiff" for the larger values of the rate constant (SUaMetNE and GEAa. 1979; BYaNEand H ~ o u ^ a s m 1987). Actual integration of these equations to obtain profiles was accomplished with the thNoM,XaSH (1972) implementation of the GEAa (1969) method for stiff equations. Secondly. by deriving equation (A 16). an equation for the critical flux. F~r. was established, i.c. equation (i2). which would otherwise not have been possible. Thirdly. it is not necessary to integrate cquation (A 16) if the only question of intcrest is the amount of op;,I preserved for a given sct of conditions. The amount of opal preserved as x ~ ~ is obtained simply by applying equation (A 16) at x -- 0 and substituting equation (AS) fi,r the dcrivative of file silica c o n c e n t r a t i o n ,

2(~.,fu-~l~;.,}. l

{{ '~

- (l - ,/')..'/,(~)1:

+

•

a,},rl, l ) ,

l/:.

- (I - ,:,)..b(~)

I

+ a.~,:b(~o)In

(I

b(~o)

((, - Cw)" = 0

(A20)

a2

which is a nonlinc:tr algebraic equation for wb( o0). After solving eq uatkm (A20) hw wb( ~ ). equation (A l ) can be used to calculate b('x), the opal preserved. Fourthly. the development of a senti-analytical solutkm will prove valuable as a means of verifying the solution given by fully numerical methods, which will be needed for the study of more general cases, i.e. depth-dependent porosity, multiple opal types, precipitation of a silicon mineral at depth, and non-steady-state calculations. Finally. analytical solutions for concentration profiles arc possible for two limiting cases. The first is when the dissolved flux out of the sediment is less than about .50% of the final burial flux. i.e. equation (13). Rewrite equation (A 16) as

! ( ~ ] : + "~ dC + 2~dxl

a, dx

e~ wt,(oo)j. a~.

+ wb(=)

= ~(C,

C)-'.

(A21)

Now employ the expansion (GaaDSHT~YN and Rvzmg. 1980, p. 44)

In

+

wb(a:"-"~

"

wb(~) ¢"~

(A22)

Modelling early diagenesis of silica

1567

to simplify equation (A21) to wb(=)/

z

which, with equation (A l ). has as its solution equation (14) of the text. The amount of opal preserved at depth. b(~=), is then found by substituting equations (4) and (10b.c) and the first derivative of equation (A~J) into equation (AS), i.¢. equation (15). The other limiting case for which analytical solutions are available is when all the opal dissolves, b(~=) = 0. Anyone interested in this case can contact the author for these solutions.