A Preliminary Comprehensive Model for the Chemistry of Sulfidic

35

A

Preliminary Comprehensive

of S u l f i d i c

Marine

Model

for the C h e m i s t r y

Sediments

LEONARD RORERT GARDNER

Downloaded by CORNELL UNIV on November 29, 2012 | http://pubs.acs.org Publication Date: March 19, 1979 | doi: 10.1021/bk-1979-0093.ch035

Department of Geology, University of South Carolina, Columbia, SC 29208

Diagenetic reactions i n subaqueous sediments commonly produce systematic v e r t i c a l changes i n the composition o f pore waters. In general these reactions involve the decomposition o f sediment organic matter, the regeneration o f nutrients and the p r e c i p i t a t i o n and d i s s o l u t i o n o f various mineral phases. As a r e s u l t o f the v e r t i c a l concentration gradients that develop i n pore waters, d i f f u s i v e exchange of dissolved substances between the sediment and overlying water can occur. In environments where the overlying water i s oxygenated, decomposition o f sediment organic matter i s c a r r i e d on by aerobic organisms that u t i l i z e dissolved oxygen d i f f u s i n g into the sediment. However, i f the oxygen demand i s greater than that supplied by d i f f u s i o n , the sediment w i l l become anoxic at some depth. At t h i s point the oxidation o f organic matter can be continued by anaerobic organisms that u t i l i z e s u l f a t e o r n i t r a t e as oxidants. In marine sediments the high concentrat i o n o f s u l f a t e i o n i n sea water enables s u l f a t e reduction to dominate the process o f anaerobic decomposition. In recent years various workers 11-7) have succ e s s f u l l y developed models based on the mathematics o f d i f f u s i o n ( 8 ) to describe v e r t i c a l p r o f i l e s o f selected chemical parameters i n marine sediments dominated by s u l f a t e reduction. Several papers ^9, 10) have also proposed models f o r nitrogen diagenesis i n the upper aerobic zone o f such sediments. Most o f these models, however, deal with only one or^two r e l a t i v e l y well behaved parameters, such as SO5 o r CO2$ which do not i n t e r a c t strongly with other components o f the sediment besides organic matter. A t r u l y comprehensive model f o r such sediment should deal simultaneously with a l l of the major chemical parameters of the system and i d e a l l y should be formulated as an i n i t i a l value prob0-8412-0479-9/79/47-093-795$05.00/0 © 1979 American Chemical Society In Chemical Modeling in Aqueous Systems; Jenne, E.; ACS Symposium Series; American Chemical Society: Washington, DC, 1979.


796

CHEMICAL MODELING IN AQUEOUS SYSTEMS

lem so that the depth at which the sediment becomes anoxic could be f r e e l y predicted rather than assumed as i n the nitrogen diagenesis models previously c i t e d . Formulation of such a model would require consideration of the k i n e t i c s of the many complex reactions that may occur near the oxidized-reduced boundary as a r e s u l t of i n t e r a c t i o n s between downward d i f f u s i n g oxygen with upward d i f f u s i n g reduced substances such as H2S and ΝΗ4· Furthermore the upper portion of most marine sediments are inhabited by burrowing organisms and thus the model would have to take into account the movement of substances r e s u l t i n g from bioturbation i n order to accurately predict the depth to the oxidized-reduced boundary. The author has attempted to formulate such a model but because of the many ad hoc assumptions that were required and because of the d i f f i c u l t i e s encoun tered i n obtaining reasonable numerical solutions t h i s work i s not ready f o r presentation. However i f one i s w i l l i n g to ignore the upper oxidized p o r t i o n of the sediment or to assume that the sediment i s anoxic a t the surface, then, as w i l l be shown below, i t i s pos s i b l e to develop a reasonably comprehensive and r e a l i s t i c model f o r sediment dominated by s u l f a t e reduction. Although previous models f o r s u l f i d i c marine sediments have s u c c e s s f u l l y simulated p r o f i l e s of s u l f a t e , carbon dioxide and ammonia, development of a model f o r H2S p r o f i l e s has been stymied by the f a c t that H2S may p r e c i p i t a t e as FeS and/or FeS2· This i n a b i l i t y to p r e d i c t H2S p r o f i l e s has i n turn prevented p r e d i c t i o n of pH p r o f i l e s . Only recently have experi mental studies of the k i n e t i c s of i r o n s u l f i d e forma t i o n ( 1 1 , 12) provided a basis f o r formulating models f o r the diagenesis o f H2S i n the presence o f i r o n oxides. This paper proposes a system of 10 non-linear, simultaneous d i f f e r e n t i a l equations (Table I) which; upon f u r t h e r development and v a l i d a t i o n , may serve as a comprehensive model f o r p r e d i c t i n g steady state, v e r t i c a l p r o f i l e s of chemical parameters i n the s u l f i d e dominated zones o f marine sediments. The major objec t i v e o f the model i s to p r e d i c t the v e r t i c a l concentra t i o n p r o f i l e s o f H 2 S . h y d r o t r i o l i t e (FeS) and p y r i t e (FeS2). As with any model there are a number of assumptions involved i n i t s construction that may l i m i t i t s a p p l i c a t i o n . In a d d i t i o n to steady state, the major l i m i t i n g assumptions of t h i s model are the assumptions that the sediment i s f r e e of CaC03, that the d i f f u s i o n c o e f f i c i e n t s of a l l dissolved s u l f u r species are equivalent and that dissolved oxygen does not penetrate into the zone o f s u l f a t e reduction.

In Chemical Modeling in Aqueous Systems; Jenne, E.; ACS Symposium Series; American Chemical Society: Washington, DC, 1979.

35.

GARDNER

Sulfidic

Marine

Sediments

797

Also I t should be emphasized that the model developed here has not been v a l i d a t e d by q u a n t i t a t i v e comparisons with comprehensive f i e l d data because of my lack of access to such data. However, as w i l l be shown, the p r o f i l e s of H 2 S , FeS and FeS generated by the model compare favorably i n a q u a l i t a t i v e fashion with such p r o f i l e s as are a v a i l a b l e i n the l i t e r a t u r e . Thus I f e e l that the work presented here represents a s i g n i f i c a n t step towards the development of a t r u l y compre hensive model f o r the chemistry of s u l f i d i c marine sediments. Downloaded by CORNELL UNIV on November 29, 2012 | http://pubs.acs.org Publication Date: March 19, 1979 | doi: 10.1021/bk-1979-0093.ch035

2

Construction of the Model The formulation of t h i s model i s based on the mathematics of diagenesis developed by Berner (£)· Using the chain r u l e of p a r t i a l d i f f e r e n t i a t i o n (13, p. 3 3 5 ) Berner showed that f o r any property of a sedi ment, P, that i s a continuous d i f f e r e n t i a b l e function of depth, x, and time, t , *E = d£ ot dt

s

6£ οχ

(1)

where S(= dx/dt) i s the r a t e of sediment deposition minus the sum of the rates of erosion and compaction. Equation 1 i s known as the general diagenetic equation and i s the basis f o r the system o f equations presented i n Table I . I t should be noted here that horizontal v a r i a t i o n s i n Ρ are assumed n e g l i g i b l e . The t o t a l d e r i v a t i v e dP/dt represents the rates of a l l the dia genetic processes a f f e c t i n g P. Omitting the usually n e g l i g i b l e e f f e c t of the flow of water r e s u l t i n g from compaction, dP/dt can be represented i n general by (2)

where D i s the whole sediment d i f f u s i o n c o e f f i c i e n t (8) and CR represents the a d d i t i v e rates of a l l chem i c a l reactions a f f e c t i n g P. P o s i t i v e terms f o r CR imply production of Ρ whereas negative terms imply consumption. In the case where Ρ i s a s o l i d substance, D can be assumed to equal zero so that there i s no e f f e c t due to d i f f u s i o n . Also since steady-state con d i t i o n s are assumed, δp/δt = 0 and thus the general diagenetic equation reduces to



798

C H E M I C A L MODELING

IN AQUEOUS

SYSTEMS

In view of the f a c t that s u l f a t e reduction i s driven by the microbial decomposition o f organic matter [(OC) i n Table i j the diagenetic equation f o r t h i s substance w i l l be derived f i r s t . Although the mechanisms by which s u l f a t e reducing b a c t e r i a decom pose sediment organic matter are complex and not f u l l y understood, Berner (2, and others (2) have success f u l l y modeled steady-state p r o f i l e s of s u l f a t e on the assumption that the k i n e t i c s of decomposition are f i r s t order and that the organic substrate i s i n s o l i d form. Thus by s u b s t i t u t i n g D == 0 and CR = -K (OC) i n equation 3 the diagenetic equation f o r organic carbon i s found as

2 per mole and i s a f u n c t i o n o f sphere diameter and the density and molecular weight o f goethite. As f o r r e a c t i o n 8 , p y r i t e formation, Rickard reports that the k i n e t i c s o f t h i s process obey the following equation1

lf& = KFS2

di

CH+)(H S)A A 2

h

S

where (FeS ) 2

55

FeS in mol Γ"

1

2

KFS2

= r a t e constant (Table II)

A

= surface area of FeS i n cnr

h


U0)


Rate constant for FeS precipitation

Rate constant for FeS precipitation

Ratio of Ν to C in organic matter

Molar surface area of HFe0

Molar surface area of S°

KFS

KFS2

α

SHM

SOM

1

2

2

cm mol^l

cm mol""-

2

none

- -1 cm6sec-1 mol

mo 1 ~" ^ 1~% cm"*2 s e c

sec"*

1

1

cm sec"

2

cm sec"

Units

- 1

2.29 χ 10

4.53 χ 10

0.155

5

5

5

1.0 χ 10-13

1.5 χ 10-7

5.0 χ io-i°

3.0 χ ΙΟ"

6

4.S χ ίο-*

ValMe

c

c

assumed

b

assumed

a

source

4

C6

C

β

(SHM) (KFS)/S

« -1.5 CSHM)CKFS)/D

K/S 2

3

C » -K/D

7

2

C = UOM) U T M ) U F s 2 ) / s

5

2

C = -0.5 (SOM)VSTM) UFs2)/ϋ

C = S/D

Molar surface area of FeS cm mol~"l 2.67 χ 1U c a. Typical of coastal sediments, southeastern U.S.A. b. Consistent v/ith Berner* s {,22) empirical relation between S and K. c. Based on spherical particles 2 \im i n diameter. DERIVED CONSTANTS t

STM

Rate constant for sulfate reduction

Κ

2

Diffusion coef·

0

2

Sedimentation rate

S

PRIMARY CONSTANTS ι Svmbol DescriDtion

Table II· Primary and derived constants for model equations.


802

CHEMICAL

A

s s

MODELING IN

AQUEOUS

SYSTEMS

surface area of S ° i n cnr.

Again assuming spherical p a r t i c l e s A^and A expressed as

can be

g

A

h

= STM

(FeS)

Ul)

A

s

=» SOM

(S°)

(12)


where STM and SOM are the surface areas per mole of FeS and S°, r e s p e c t i v e l y . Thus i n summary the CR expression f o r H2S i s given by KFS(H ) (H S) / SHM(HFe0 ) +

CR - 0.5(OC) - 1.5 -0.5

2

3

2

2

2

KFS2(H+)(H S)STM (FeS) SOM(S°) 2

2

(13)

2

where the c o e f f i c i e n t s , -1.5 and -0.5, are derived from the stoichiometry of reactions 6 and 7. Again i t i s l e f t to the reader to v e r i f y that s u b s t i t u t i o n of the above r e s u l t f o r CR i n equation 3 leads to the d i f f e r e n t i a l equation f o r t o t a l dissolved H2S shown i n Table I . From t h i s d i s c u s s i o n the o r i g i n of the d i f f e r e n t i a l equations f o r H F e 0 , FeS, S° and FeS2 should also be obvious· With regard to the stoichiometry o f reactions 6 and 7 i t should be mentioned that Rickard was unable to a s c e r t a i n the exact mechanisms and reactions underlying h i s k i n e t i c r e s u l t s . I have written r e a c t i o n 6 as an o v e r a l l d e s c r i p t i o n of the process of FeS formation because the c o e f f i c i e n t s of each reactant i n equation 6 matches i t s corresponding exponent i n equation 85. Furthermore r e a c t i o n 6 i s compatible with R i c k a r d s i n t e r p r e t a t i o n of the k i n e t i c mechanisms involved i n t h i s process and with h i s report of the production o f native s u l f u r i n h i s experiments. In the case of reaction 7, however, the achievement of harmony between k i n e t i c s and stoichiometry i s more problematic. Again Rickard was unable to a s c e r t a i n the exact mechanisms and reactions involved i n p y r i t e formation. He d i d not attempt to write an o v e r a l l r e a c t i o n but suggested that the process involves d i s s o l u t i o n o f both FeS and S°, the formation of p o l y s u l f i d e ions and t h e i r subsequent r e a c t i o n with F e * ions to form p y r i t e plus native s u l f u r . As with r e a c t i o n 6 I have chosen to write react i o n 7 i n such a fashion that the stoichiometric c o e f f i c i e n t s of the reactants correspond with t h e i r respective exponents i n equation 10· I f one assumes that R i c k a r d s r e a c t i o n v e s s e l s were devoid, as seems l i k e l y , of strong o x i d i z i n g agents such as f e r r i c i r o n 2

1

2

1


35.

GARDNER

Sulfidic Marine Sediments

803

or 0 and that the N? c a r r i e r gas used i n h i s experi ments was not reduced to ammonia, then the constraints of reaction 7 force one to postulate the reduction of H to Η2· This i s a somewhat unattractive requirement i n view of the lack of supporting experimental evidence f o r Ho production but i t i s at l e a s t thermodynamically f e a s i b l e i n that the equilibrium concentration of aqueous H2 would have to be at l e a s t 10^ times greater than the concentration of H2S. Given t y p i c a l concen t r a t i o n s of H2S i n anoxic sediments i t i s u n l i k e l y that an equilibrium pressure of H2 could develop at one atmosphere t o t a l pressure. I f one i s w i l l i n g to admit the presence of O2 i n the system, then r e a c t i o n 7 could be modified as follows 1 2


+

2FeS(s) + S°(s) 2FeS2(s) +

+ H

+

+ HS~(aq) + kOz

-> (14)

H 0. 2

T h i s , however, i s undesirable because O2 does not occur i n the rate equation and because, as assumed, i t i s u n l i k e l y that there i s s u f f i c i e n t O2 present i n sedi ments at depths where s u l f a t e reduction and p y r i t e formation are a c t i v e . Furthermore the often c i t e d reaction between FeS and S°, FeS(s) + S°(s) -» FeS2(s)

(15)

i s unappealing because i t bears no obvious r e l a t i o n to the rate equation and because i n conjunction with reac t i o n 6 i t would upon completion ( i . e . at a s u f f i c i e n t l y great depth) produce FeS and FeS2 i n at l e a s t a one to one mole r a t i o which contradicts f i e l d studies that commonly show almost t o t a l conversion o f FeS to FeS2 at depth. On the other hand the combination of reactions 6 and 7 allow but, as w i l l be seen, do not demand com p l e t e conversion to FeS2» Although other reactions involving s i d e r i t e formation, methane formation, N£ reduction and/or reactions with organic complexes have been proposed f o r p y r i t e formation (14b) none of them have been studied k i n e t i c a l l y and thus siinply cannot be considered at t h i s time. F i n a l l y although Rickard reports that framboidal p y r i t e was not produced during his experiments t h i s does not n e c e s s a r i l y mean that framboidal p y r i t e formation does not obey the s t o i c h i ometry and k i n e t i c s of equations 7 and 10· At t h i s point construction of the model i s com p l e t e except f o r incorporation of H* into the system of equations. As discussed above,H i s an important v a r i a b l e i n the k i n e t i c s of i r o n s u l f i d e formation. +



804

CHEMICAL

MODELING IN

AQUEOUS

SYSTEMS

However because the rate law expressions f o r i r o n s u l f i d e formation (equations 8 and 10) are non-linear the d i f f e r e n t i a l equations f o r H2S and the i r o n s u l f i d e s are not amenable to e x p l i c i t s o l u t i o n . Thus i t i s important to develop an equation f o r ¥r that can be incorporated i n a numerical s o l u t i o n technique such as that of Runge-Kutta (U>)· Fortunately an appropriate d i f f e r e n t i a l equation f o r H* can be developed from charge balance considerations. Here i t i s assumed that dissolved substances other than those l i s t e d i n Table I are not a f f e c t e d by diagenesis. I f t h i s i s true, then a charge balance d i f f e r e n c e equation can be written (16): (NH4)

- (Ν*φ

0

+ (H+)

+

- (H )

0

= (OH")

- (0H")

+ 2(S0|")

- 2(SO^)

Q

+ (HC05) - ( H C 0 5 )

+ 2(C0|")

- 2(C0^")

Q

+ (HS~) - (HS~)

o

o

0

(16)

where the subscripted q u a n t i t i e s represent the concen t r a t i o n s o f these substances i n the overlying sea water. Here i t should be pointed out that i n develop ing the equations i n Table I i t has been t a c i t l y assumed that pH dependent species of C02t H2S and NH3 (such as HC0§, HS", etc.) have equivalent d i f f u s i o n c o e f f i c i e n t s . In p r i n c i p l e a separate d i f f e r e n t i a l equation could have been w r i t t e n f o r each species. However each of these would incorporate a pH dependent d i s t r i b u t i o n function. The programing complexity that t h i s would generate i s not j u s t i f i e d given the prelim inary nature of the model and the other rough assump tions involved. In equation 16 each of the carbonate species can be expressed i n terms of t o t a l CO2 and H* by means o f the appropriate carbonate equilibrium con stants. The other pH dependent species such as HS* can be treated i n a s i m i l a r fashion. In t h i s manner equa t i o n 16 can be expressed i n terms o f the v a r i a b l e s (H S), ( S 0 ) , (CO2), (NH ) and (H+), the i n i t i a l con centrations and the various equilibrium constants. Using the chain-rule of p a r t i a l d i f f e r e n t i a t i o n (12, p. 335) the t o t a l d i f f e r e n t i a l o f (H ) can be obtained and i s o l a t e d on the l e f t hand of the equation. D i v i sion o f both sides of the r e s u l t i n g equation by dx y i e l d s a f i r s t order d i f f e r e n t i a l equation o f the form 2

4

4

+

= f[(H S), (S0 ), (CO2), (NH4), 2

4

(H S)\ 2

( S 0 ) , ( C 0 ) \ (NH4) J e

4

1

2


(17)


35.

GARDNER


805

which can be incorporated i n a Runge-Kutta algorithm. Space l i m i t a t i o n s unfortunately do not permit a more e x p l i c i t d e r i v a t i o n of t h i s equation than that sketched here. Development of t h i s model has neglected several substances that i n p r i n c i p l e p o s s i b l y could have been included. The most obvious omissions are dissolved phosphate released by the decomposition of organic matter and the p o s s i b l e d i s s o l u t i o n or p r e c i p i t a t i o n of CaC03. Both of these substances would a f f e c t the pH p r o f i l e . The e f f e c t of phosphate on pH would probably be small because i t s concentration i n organic matter i s low compared to Ν and C. The e f f e c t of CaC03 on pH however could be important and thus a p p l i cation of the model should probably be r e s t r i c t e d to CaC03 free sediments. The Simulation Algorithm The equations of the system presented i n Table I are both non-linear and simultaneous. Thus there i s no obvious way of obtaining e x p l i c i t solutions and one must therefore resort to numerical techniques. Because of i t s e f f i c i e n c y and wide use, I chose a Runge-Kutta method with fourth-order accuracy (1£)· The d e t a i l s of the algorithm developed f o r t h i s problem w i l l eventu a l l y be published i n an appropriate j o u r n a l . The algorithm can i n p r i n c i p l e handle the system of equa tions presented i n Table I . In p r a c t i c e however prob lems a r i s e because some of the equations are second order and thus the algorithm must be supplied with i n i t i a l values ( i . e . at χ = o) f o r both concentration and the f i r s t d e r i v a t i v e . As w i l l be discussed below the algorithm may generate u n r e a l i s t i c concentrations i f inappropriate values are chosen f o r the f i r s t deriv a t i v e . I t should be emphasized that t h i s does not r e s u l t from d e f i c i e n c i e s i n the algorithm, such as error buildup, because such e f f e c t s can be minimized by choosing s u f f i c i e n t l y small depth increments. The performance of the algorithm with respect to error buildup was evaluated by having the algorithm compute p r o f i l e s f o r a hypothetical sediment devoid of r e a c t i v e i r o n oxide. For such a case the equations i n Table I e i t h e r vanish or uncouple and become l i n e a r so that e x p l i c i t solutions can be obtained f o r each ^except H" ) · Using e x p l i c i t l y determined i n i t i a l values f o r the required f i r s t derivatives and a step s i z e of 0.05 cm, a l l of the numerically computed values agreed with those computed from t h e i r corresponding e x p l i c i t func tions to at l e a s t four d i g i t s to a depth of 100 cm. h



806

CHEMICAL

MODELING

IN

AQUEOUS

SYSTEMS

Furthermore at each step the r a t i o of the change i n c a t i o n i c charge to the change i n anionic charge (see equation 16) was 1.0000. Thus i t appears that the c h i e f Limitations of the algorithm are those due to p o s s i b l e d e f i c i e n c i e s i n the conceptual framework of the model and to the appropriateness of the i n i t i a l conditions supplied to the algorithm. The problems associated with c o r r e c t l y s p e c i f y i n g the i n i t i a l f i r s t d e r i v a t i v e values became evident luring e a r l y attempts to simulate p r o f i l e s of H2S and Lron s u l f i d e s . As pointed out e a r l i e r the d i f f e r e n t i a l equation f o r organic carbon has an e x p l i c i t s o l u t i o n . This s o l u t i o n can be substituted f o r (OC) i n the equa tions f o r S04, NH4 and CO2 which are l i n e a r and can thus be solved a n a l y t i c a l l y . Thus the e x p l i c i t s o l u tions obtained can be used to c a l c u l a t e the c o r r e c t I n i t i a l f i r s t d e r i v a t i v e s f o r these parameters. The equation f o r H2S, however, i s non-linear and coupled to :hose f o r HFe02, FeS, S° and FeS2« Thus there i s no )bvious a n a l y t i c a l way to c o r r e c t l y s p e c i f y the i n i t i a l i e r i v a t i v e f o r H2S. However since FeS p r e c i p i t a t i o n :an not occur at the zero concentration of H2S assumed it χ o, i t seemed reasonable to assume that at χ ο :he d e r i v a t i v e of H2S should be equal i n magnitude but >pposite i n sign to that of S04» U t i l i z a t i o n of t h i s issumption i n the algorithm r e s u l t e d i n H2S coneentra vons that g r e a t l y exceed the concentrations of H2S :hat would be produced i n a system devoid of r e a c t i v e .ron oxide but otherwise s i m i l a r . Furthermore simula:ion runs showed that although the algorithm obeyed lass balance constraints on Fe i t i n v a r i a b l y converts it depth e s s e n t i a l l y a l l of the i n i t i a l l y supplied ^ a c t i v e HFe02 to F e S . In cases where the i n i t i a l concentration of organic carbon i s l e s s than four times :he i n i t i a l concentration o f HFe02, t h i s behavior of :he algorithm i s also u n r e a l i s t i c because according to •eaction 5 the maximum number o f moles o f reduced s u l "ur i n a l l forms (H2S, FeS, FeS2t S°) that can be proluced i s equal to one h a l f of the i n i t i a l number moles >f decomposable carbon. Thus one i s faced with the »roblem of e i t h e r f i n d i n g an i n i t i a l d e r i v a t i v e f o r hS that gives acceptable r e s u l t s or incorporating into he algorithm the appropriate organic carbon substrate ionstraint on i r o n s u l f i d e formation. The f i r s t a l t e r lative was attempted on a t r i a l and e r r o r basis but oon proved f u t i l e . Fortunately however I was able to iodify the algorithm by incorporation o f an appropriate rganic carbon c o n s t r a i n t . The carbon constraint on i r o n s u l f i d e formation i s ased on r e a c t i o n 5 and Berner s (12) demonstration s

3

2

1


35.

807


GARDNER

that the t o t a l s u l f a t e reduced during b u r i a l to depth x, SR, i s given by L S where ( 0 C ) = i n i t i a l concentration of organic carbon. Since the s o l u t i o n of the d i f f e r e n t i a l equation f o r (OC) i s o

(OC) = ( O C ) ^ e x p ( ^ * ) ]

U9)

i t follows that SR = U . 5 [ ( 0 C ) - (OC)]* In the case of a sediment devoid of i r o n oxide, such that none of the reduced s u l f a t e i s incorporated i n s o l i d phases such as FeS, the concentration of t o t a l dissolved HoS equals (SO^)^ - (SO4) providing the d i f fusion c o e f f i c i e n t s f o r H S and SO4 are equivalent. Thus i t i s reasonable to assume that i n a system where h a l f of the reduced s u l f a t e has been incorporated i n s o l i d phases as a r e s u l t of reactions with i r o n oxide, the concentration of H S should be equal to one h a l f of H S concentration that would p r e v a i l i n a system devoid 0 1 i r o n oxide but otherwise s i m i l a r . This idea can be expressed mathematically as


o

2

2

2

SR - 2(FeS ) - (FeS) - (S°)

(H S 2

(S0)

o

-

2

(SO4)

=

SR

< 2 0 )

*

Thus i f no i r o n s u l f i d e genesis occurs because of the absence of i r o n oxide, then (H S) = (SO^Jp - (SO4)· Otherwise (H S) w i l l be some f r a c t i o n of ( S 0 4 ) - (SO4) depending on the amount of SR that has not been incorporated into s o l i d phases as reduced s u l f u r . This constraint on H S and i r o n s u l f i d e genesis was incorporated into the o r i g i n a l algorithm i n the following manner. At the end of the Runge-Kutta c a l c u l a t i o n s f o r each depth increment the following tolerance r a t i o was computed* 2

o

2

2

SRC

- 2(FeS2)

» [(SR

- (FeS)

-

(S°)

SR(H2S) (S0 ) 4

o

-

(SO4)

J/SR

.

(21)

I f the carbon constraint i s obeyed exactly SRC equals zero. I f SRC l i e s between -0.05 and +0.05 the algorithm accepts the r e s u l t s and proceeds to the next increment. I f SRC l i e s outside t h i s range, the algorithm repeats the c a l c u l a t i o n s f o r the increment with


CHEMICAL

808

MODELING

IN

AQUEOUS

SYSTEMS

appropriately adjusted values of (H S)' u n t i l tolerance is satisfied. 2

Simulation of Iron S u l f i d e Genesis In order to i l l u s t r a t e the q u a l i t a t i v e agreement of the carbon constrained model developed here with f i e l d data c i t e d i n the l i t e r a t u r e , several simulations of H2S, FeS, FeS and pH p r o f i l e s are presented i n Figures 1 and 2· Because of t h e i r appearance i n the r a t e equations f o r i r o n s u l f i d e k i n e t i c s these simulations have a l s o been designed to i l l u s t r a t e the poss i b l e e f f e c t s of pH and reactant molar surface areas on i r o n s u l f i d e genesis. In both figures s u l f a t e reduction and i r o n s u l f i d e formation begin about one cm below the sediment surface. The p r o f i l e s of pH above t h i s depth have been simulated on the basis of a model developed (but not discussed here) f o r the o x i dized zone of the sediment. This model incorporates reactions pertaining to the aerobic decay of organic matter, bioturbation and nitrogen diagenesis. Since reactions i n the aerobic zone are assumed to be uncoupled from those i n the zone o f s u l f a t e reduction the model f o r the aerobic zone has been used here simply as a means o f varying the i n i t i a l pH at the top of the zone of s u l f a t e reduction. As can be seen,the H2S p r o f i l e s f o r a l l four simul a t i o n s shown i n Figures 1 and 2 show maximum concent r a t i o n s i n the submillimolar range and subsequent decreases with depth. In sediment hosting a c t i v e s u l f a t e reduction and p y r i t e formation, H2S concentrations a t t a i n maximum values of about 10"3 moles l i t e r * 1 and commonly decrease t h e r e a f t e r with depth (IS, 12). In general incorporation of the carbon constraint i n the algorithm produces simulations that agree with q u a l i t a t i v e predictions of H2S p r o f i l e s by Goldhaber and Kaplan (18, p. 638). The simulations o f FeS shown by jobs 433, 1502 and 1673 are s i m i l a r to published p r o f i l e s (12, 20) i n that they show abrupt peaks j u s t below the oxidized-reduced boundary. Also, as i s commonly the case i n natural sediments, the simulated peak concentrations of FeS are about an order o f magnitude lower than the ultimate concentration of FeS2» Thus i n both magnitude and shape the simulations shown here bear strong resemblance to those r e c e n t l y described i n the l i t e r a t u r e , suggesting that the carbon constrained model has p o t e n t i a l a p p l i c a t i o n to the study of s u l f u r diagenesis i n marine sediments. The simulations shown i n Figure 1 are designed to i l l u s t r a t e the possible e f f e c t of the i n i t i a l pH at the


2



Figure 1. Simulations of sediment chemistry showing effects of differences in pH at the top of the zone of sulfate reduction. Values for model constants are listed in Table II unless otherwise on figure. Ρ and ANP are parameters of the aerobic model that control the depth and pH at the oxidized-r educed boundary.


oo

810

CHEMICAL MODELING IN AQUEOUS

SYSTEMS


:1 .si CO

KH

co

Ο Co Qù ^

be β •Si

III S, s * ν. CJ

s-

I §

•Si

-s 8

is

ο .

Ci)

co

Ο

tuD §

fi δ



35.

GARDNER


811

s t a r t of s u l f a t e reduction on the course of i r o n s u l f i d e genesis. The pH at the oxidized-reduced boundary must be c o n t r o l l e d at l e a s t i n part by reactions i n the oxidized zone, such perhaps as ammonia oxidation. Thus i t i s conceivable that a comprehensive understanding of i r o n s u l f i d e genesis might require an understanding of the factors that control the pH at t h i s boundary. In both simulations the surface area parameters (SHM, STM, SOM) and the i n i t i a l organic carbon concentration are i d e n t i c a l . Thus the only d i f f e r e n c e between the simul a t i o n s as f a r as i r o n s u l f i d e genesis i s concerned i s that i n job 1502 the s t a r t i n g pH i s 6.47 whereas i n job 1673 i t i s 7.75. As can be seen t h i s d i f f e r e n c e i n s t a r t i n g pH has no e f f e c t on the p r o f i l e of FeS and only s l i g h t l y a l t e r s the maximum FeS concentration. Although the r e s u l t s of t h i s simulation experiment are preliminary they suggest that the pH at the oxidizedreduced boundary i s not l i k e l y to be an important factor i n i r o n s u l f i d e genesis. This r e s u l t s from the fact that the rapid depletion of SOç ion j u s t below the oxidized-reduced boundary i s not compensated by an adequate increase i n dissolved H^S. Thus charge balance consideration require a rapid r i s e i n pH which soon d i s s i p a t e s the i n i t i a l k i n e t i c advantage of low pH. The simulations shown i n Figure 2 are designed to i l l u s t r a t e the e f f e c t of surface area. In both, the depth to the zone of s u l f a t e reduction as well as the i n i t i a l conditions at the beginning of s u l f a t e reduct i o n are i d e n t i c a l except f o r the assumed molar surface areas. In job 433 the molar surface areas of S°, SOM, and FeS, STM, are about an order of magnitude l a r g e r than the molar surface area of HFe02» SHM. Although the magnitudes of SOM, STM and SHM chosen f o r job 433 are l a r g e r than f o r jobs 1502 and 1673 (Figure 1) the r a t i o s of SOM and STM to SHM are s i m i l a r i n a l l three jobs. As can be seen by comparing Figures 1 and 2 the p r o f i l e s of FeS and FeS i n job 433 are very s i m i l a r to those f o r jobs 1502 and 1673. On the other hand i n job 230 (Figure 2) the molar area of HFe02 i s about an order of magnitude l a r g e r than the molar areas of S° and FeS so that the r a t i o s of SOM and STM to SHM are about one hundred times smaller than they are i n jobs 1502, 1673 and 433. This dramatic d i f f e r e n c e r e s u l t s i n the p r e c i p i t a t i o n of approximately equal amounts of FeS and FeS2 at depth i n job 230. Thus the s i z e (and r e a c t i v i t y ) of i r o n oxide p a r t i c l e s supplied to the sediment could be an important f a c t o r i n producing FeS enriched sediments such as those discussed by Berner ( 8 ) from the Black Sea. 2

2


812

CHEMICAL

MODELING IN AQUEOUS

SYSTEMS


Concluding Remarks Although the above discussion does not constitute a rigorous v e r i f i c a t i o n o f the model ( £ L ) · the simil a r i t y of the simulations presented here to observed p r o f i l e s o f H 2 S , FeS and FeS 2 i n marine sediments i n d i cates that the model may be o f value i n the study and i n t e r p r e t a t i o n of v e r t i c a l patterns i n s u l f u r diagenes i s · Comprehensive multiparameter analyses o f sediment p r o f i l e s from a v a r i e t y o f s i t e s w i l l be required to v a l i d a t e the model* In t h i s endeavor techniques w i l l have to be devised to a s c e r t a i n the molar surface areas of the various s o l i d phase reactants. Eventually i t may be p o s s i b l e to expand the model presented here to include processes i n the aerobic zone so that the depth to the oxidized-reduced boundary can be predicted as well as the pH p r o f i l e through t h i s boundary. This achievement would c o n s t i t u t e a t r u l y comprehensive model. Acknowledgements The author g r a t e f u l l y acknowledges the assistance of R. Underwood and G. Johnson who guided the author i n developing the Runge-Kutta algorithm. W. £· Sharp, W. S . Moore and M. Goldhaber reviewed the manuscript and provided h e l p f u l suggestions and discussions.

Abstract This paper proposes a system of ten non-linear, simultaneous differential equations which, upon further development and validation, may serve as a comprehensive chemical model for sulfide dominated marine sediments. These equations were developed in accordance with R. A. Berner's mathematics of diagenesis and incorporate D. T. Rickard's work on the kinetics of iron sulfide formation. A Runge-Kutta algorithm has been developed to provide numerical solutions for the equations. When utilized with an organic carbon substrate limitation on the production of reduced sulfur, the algorithm generates profiles of H S, FeS and FeS that agree qualitatively with measured profiles reported in the recent literature. Experiments with the algorithm suggest that the ratio of FeS to FeS at depth depends strongly on the i n i t i a l molar surface area of goethite but that the profiles of FeS, FeS , and H S are not greatly affected by the i n i t i a l pH. 2

2

2

2

2


35.

GARDNER


813

Literature Cited 1. 2.


3.

4. 5.

6. 7. 8. 9.

10.

11. 12. 13. 14.

Berner, R.A. Stoichiometric models for nutrient regeneration in anoxic sediments. Limnol. Oceanogr. 22, 781-786 (1977). Berner, R.A. Diagenetic models of dissolved species in the i n t e r s t i t i a l waters of compacting sediments. Amer. Jour. Sci. 275, 88-96 (1975). Berner, R.A. Kinetic models for the early diagenesis of nitrogen, sulfur, phosphorus, and silicon in anoxic marine sediments, p. 427-450, in Goldberg, E.D. (ed.), "The Sea," v o l . 5, John Wiley & Sons, New York, 1974. Berner, R.A. An idealized model of dissolved sulfate distribution in recent sediments. Geochim. Cosmochim. Acta 28, 1497-1503 (1964). Martens, C.S. and Berner, R.A. Interstitial water chemistry of anoxic Long Island sound sediments. l. Dissolved gases. Limnol. Oceanogr. 22, 10-25 (1977). Lasaga, A.C. and Holland, H.D. Mathematic aspects of non-steady-state diagenesis. Geochim. Cosmochim. Acta 40, 257-266 (1976). Toth, D . J . and Lerman, A. Organic matter reactivity and sedimentation rates in the ocean. Amer. Jour. Sci. 277, 465-485 (1977). Berner, R.A. "Principles of Chemical Sedimentology" 204 p. McGraw-Hill Book Co., New York, 1971. Vanderborght, J., Wollast, R., and Billen, G. Kinetic models of diagenesis in disturbed sediments. Part 2. Nitrogen diagenesis. Limnol. Oceanogr. 22, 794-803 (1973). Vanderborght, J. and Billen, G. Vertical d i s t r i bution of nitrate concentration in i n t e r s t i t i a l water of marine sediments with n i t r i f i c a t i o n and denitrification. Limnol. Oceanogr. 20, 953-961 (1975). Rickard, D.T. Kinetics and mechanism of the sulfidation of goethite. Amer. Jour. Sci. 274, 941-952 (1974). Rickard, D.T. Kinetics and mechanism of pyrite formation at low temperature . Amer. Jour. Sci. 275, 636-652 (1975). Hildebrand, F.B. "Advanced Calculus for Applications, " 646 p. Prentice-Hall, Inc., Englewood C l i f f s , New Jersey, 1965. Berner, R.A. Sedimentary pyrite formation. Amer. Jour. Sci. 268, 1-23 (1970).


814

15. 16. 17.


18. 19. 20.

21.

22.

CHEMICAL

MODELING IN AQUEOUS

SYSTEMS

Froberg, C.D. "Introduction to Numerical Analysis," 433 p. Addison-Wesley Publishing Co., Reading, Massachusetts, 1969. Gardner, L.R. Chemical models for sulfate re duction in closed anaerobic marine environments. Geochim. Cosmochim. Acta 37, 53-68 (1973). Berner, R.A. Sulfate reduction, pyrite formation, and the oceanic sulfur budget, p. 347-361, in Dryssen, D. and Jagner, D., ed., "The Changing Chemistry of the Oceans," Nobel Symposium 20, Almquist and Wiksell, Stockholm, 1972. Goldhaber, M.B. and Kaplan, I.R. The sulfur cycle, p. 569-655, in Goldberg, E . D . , ed. "The Sea," vol. 5, John Wiley & Sons, New York, 1974. Jorgensen, B.B. The sulfur cycle of a coastal marine sediment (Limfjorden, Denmark). Limnol. Oceanogr. 22, 814-831 (1977). FOAM (FriencTs of Anoxic Mud): Goldhaber, M.B., Allen, R . C . , Cochran, J . K . , Rosenfeld, J . K . , Martens, C . S . , and Berner, R.A. Sulfate reduction, diffusion, and bioturbation in Long Island sound sediments: Report of the Foam Group. Amer. Jour. Sci. 277, 193-237 (1977). Behrens, J.C., Beyer, J.E., Madsen. O.Β.G. and Thomsen, P.G. Some aspects of modelling the long-term behaviour of aquatic ecosystems. Ecological Modelling 1, 163-198 (1975). Berner, R.A. Sulfate reduction and the rate of deposition of marine sediments. Earth Planet. Sci. Lett. 37, 492-498 (1978).

RECEIVED November16,1978.


A Preliminary Comprehensive Model for the Chemistry of Sulfidic

Recommend Documents