manuscript submitted to
Geophysical Research Letters
Balancing the marine sulfur cycle
1
D.L. Johnson
1
∗
, J.F. Adkins
1
2
1
Division of Geological & Planetary Sciences, California Institute of Technology, 1200 E California Blvd,
3
Pasadena, CA 91125
4
Key Points:
5
•
Cluster analysis identifies geographically distinct pore-water [SO
2–
4
] profiles
6
•
Pyrite literature compilation confirms shelf and slope dominate global burial
7
•
Sampling and measurement biases can account for prior imbalance in S flux es-
8
timates
9
∗
Current address: Department of Earth, Environmental and Planetary Sciences, Rice University, 6100
Main St., MS-126, Houston, TX 77005
Corresponding author: D.L. Johnson,
Dan.Johnson@rice.edu
–1–
ESSOAr | https://doi.org/10.1002/essoar.10511114.1 | CC_BY_NC_4.0 | First posted online: Mon, 18 Apr 2022 11:14:37 | This content has not been peer reviewed.
manuscript submitted to
Geophysical Research Letters
Abstract
10
Sulfate (SO
2–
4
) reduction is a major SO
2–
4
output flux from Earth’s oceans, but an im-
11
balance between recent estimates of this flux and corresponding input fluxes suggests in-
12
accuracy in our understanding. Here, we combine global geographic trends in aqueous
13
and dissolved-phase sedimentary data to resolve this inaccuracy. [SO
2–
4
] profiles from 700+
14
sites partition into geographically-distinct
k
-means clusters based on net sulfate reduc-
15
tion rate (nSRR). Pairing nSRRs with literature-derived pyrite accumulation rates con-
16
firms that shelf and slope pyrite burial dominate burial globally. Our results also sug-
17
gest that sampling and measurement biases have led to erroneously high prior global out-
18
put estimates and can account for the flux imbalance. Disparate mean
δ
34
S values for
19
shelf versus deeper ocean pyrite indicate that sea level change may be an overlooked mech-
20
anism for forcing past changes in seawater
δ
34
S.
21
Plain Language Summary
22
The marine sulfur (S) cycle helps regulate Earth’s atmospheric O
2
concentrations
23
and is commonly assumed to operate in a steady-state condition over millions of years.
24
However, recent input and output flux estimates differ by a factor of three, suggesting
25
a drawdown of marine sulfate (SO
2–
4
) concentrations. We reexamine one of the main out-
26
put fluxes (pyrite, FeS
2
) by combining pore water and solid phase data from seafloor sed-
27
iments. Distinct geographic trends in these data allow us to make a new output flux es-
28
timate compatible with a concentration steady-state in the marine S cycle. Our anal-
29
ysis suggests that sea level changes may have affected our records of the ancient S cy-
30
cle and resulting estimates of past O
2
concentrations.
31
1 Introduction
32
The marine sulfur (S) cycle plays critical roles in regulating ocean-atmosphere O
2
concentrations (Berner, 2005) and balancing CO
2
releases associated with terrestrial sul-
fide oxidation (Torres et al., 2014). The burial of iron sulfide minerals (e.g. pyrite, FeS
2
)
formed through the reaction of sedimentary iron with aqueous sulfide derived from sul-
fate (SO
2–
4
) reduction acts as a key process within both of these roles. Here, electrons
are transferred from organic matter to sulfate via the following simplified stoichiometry:
2 CH
2
O + SO
2–
4
+ 2 H
+
H
2
S + 2 CO
2
+ 2 H
2
O
(1)
H
2
S + Fe
2+
+ S
0
FeS
2
+ 2 H
+
(2)
Although previous studies (Alt & Burdett, 1992; Puchelt & Hubberten, 1980; Zak et al.,
33
1980; Sweeney & Kaplan, 1980; Goldhaber & Kaplan, 1974; Scharff, 1980; Bonnell & An-
34
derson, 1987; B ̈ottcher et al., 2004; A. Migdisov et al., 1980; A. A. Migdisov et al., 1983;
35
Lein et al., 1976; Krouse et al., 1977; Johnston et al., 2008; Vinogradov et al., 1962; Ka-
36
plan et al., 1963; Sweeney, 1972; Thode et al., 1960; Kaplan et al., 1960; Lein, 1983; Hart-
37
mann & Nielsen, 2012; A. A. Migdisov et al., 1974) have revealed significant diversity
38
in the character of S cycling across global deep ocean sediments, most have focused on
39
samples collected from a narrow geographic range. A need for better integration of data
40
from multiple regions into assessments of the modern marine S cycle is underscored by
41
a factor of
∼
3 imbalance between recent global net SO
2–
4
reduction rate (nSRR) esti-
42
mates (11.3 Tmol/yr (Bowles et al., 2014)) and the pre-anthropogenic riverine input flux
43
of SO
2–
4
(2
.
8
±
0
.
4 Tmol/yr; (Burke et al., 2018)). If the pre-anthropogenic S cycle was
44
approximately in steady state, such imbalance suggests that this nSRR estimate is too
45
high and that prior S output estimates (Berner, 1982) are more accurate. Alternatively,
46
the marine S cycle could be out of steady state or input fluxes from rivers and assumed
47
minor sources (e.g. volcanism/hydrothermalism) could be grossly underestimated.
48
Better understanding of the modern marine S cycle is also important for reconstruct-
49
ing past S cycle changes. An isotopic fractionation associated with microbial SO
2–
4
re-
50
–2–
ESSOAr | https://doi.org/10.1002/essoar.10511114.1 | CC_BY_NC_4.0 | First posted online: Mon, 18 Apr 2022 11:14:37 | This content has not been peer reviewed.
manuscript submitted to
Geophysical Research Letters
duction (MSR) typically depletes product aqueous sulfide in the rare heavy isotope
34
S
51
(Kaplan & Rittenberg, 1964; Wing & Halevy, 2014). This
34
S-depletion is passed on to
52
sulfide minerals that form an important marine output flux. The observed SO
2–
4
-sulfide
53
δ
34
S offset has commonly been utilized create one-box steady-state model reconstruc-
54
tions of the relative amount of reduced S buried throughout Earth’s history (Canfield,
55
2004; Kampschulte & Strauss, 2004; Canfield & Farquhar, 2009). Although recent work
56
(Pasquier et al., 2017, 2021) has demonstrated local variability in this offset due to phys-
57
ical sedimentary parameters, few studies have interrogated the quantitative impacts of
58
this variability on a global scale. Furthermore, studies of the offset across time (Leavitt
59
et al., 2013; Wu et al., 2014) have relied almost exclusively on shelf sediments due to preser-
60
vation biases in the geologic record and may underestimate the true magnitude of the
61
offset if deep marine pyrite burial significantly affects isotope mass balance.
62
Here, we re-examine global SO
2–
4
reduction and the diversity of sedimentary S cy-
63
cling across deep ocean environments by applying cluster analysis to pore water [SO
2–
4
]
64
profiles from 700+ DSDP, ODP, and IODP sites. We subsequently construct a litera-
65
ture compilation of 300+ paired pyrite abundance +
δ
34
S measurements and 1200+ to-
66
tal sulfur abundance measurements from recent marine sediments to explore the coher-
67
ence of sedimentary data with existing global nSRR estimates and the potential effects
68
of depositional environment shifts on isotope mass balance in the marine S cycle.
69
2 Materials and Methods
70
2.1 Cluster analysis
71
To examine S cycling in deep ocean sediments across the globe, we downloaded data
72
collected during previous DSDP, ODP, and IODP cruises from the JANUS (
http://iodp
73
.tamu.edu/janusweb/links/links\
all.shtml
) and LIMS (
http://web.iodp.tamu
74
.edu/LORE/
) databases. These data included all available geochemical data on inter-
75
stitial waters, gases, and solid phases as well as physical properties data on sedimentary
76
porosity. Once downloaded, data were imported into MATLAB
®
and saved into MATLAB
®
77
structures.
78
Applying
k
-means clustering to pore water profiles requires interpolation of all data
79
to the same depth resolution over an identical depth interval. To enable cluster analy-
80
sis, pore water [SO
2–
4
] data were extracted from the structures and interpolated to an
81
identical depth resolution of 1 m across all sites. Due to variation in the length of core
82
over which pore water data were collected, we experimented with several bottom bound-
83
ary depths to find an appropriate balance between the total number of profiles included
84
in the analysis (which is reduced as the bottom boundary extends beyond the maximum
85
depth of pore water data) and the length of the depth interval for clustering; our choice
86
of a bottom boundary of 100 mbsf reduces the number of sites included in the analysis
87
from 780 to 729, but makes the differences between the clusters more visually clear. Fol-
88
lowing interpolation, [SO
2–
4
] data were inserted into a master array in which each row
89
represented data from a different site and each column data from a specific depth (i.e.,
90
column 1 houses all [SO
2–
4
] for 0 mbsf at different sites, column 2 houses 1 mbsf data,
91
etc.). Sites with
<
five pre-interpolation [SO
2–
4
] data points or with no [SO
2–
4
] data be-
92
low the chosen bottom boundary were excluded from our analysis. Outlying [SO
2–
4
] mea-
93
surements differing more than 10 mM from measurements above and below the data point
94
were also removed before clustering.
95
After creation of the master [SO
2–
4
] array, we used MATLAB
®
’s evalclusters func-
tion to evaluate the optimal number of clusters based on the entries in the [SO
2–
4
] ar-
ray. We chose the optimal number of clusters by maximizing the Calinski-Harabasz In-
dex (also known as the variance ratio index). This index ratios the overall variance (sum
–3–
ESSOAr | https://doi.org/10.1002/essoar.10511114.1 | CC_BY_NC_4.0 | First posted online: Mon, 18 Apr 2022 11:14:37 | This content has not been peer reviewed.
manuscript submitted to
Geophysical Research Letters
of the squares) of observations between clusters to the variance within the clusters:
Calinski-Harabasz Index =
SS
B
SS
W
x
N
−
k
k
−
1
(3)
where
SS
B
is the variance between the clusters,
SS
W
is the variance within the clusters,
96
N
is the number of observations (sites in this case), and
k
is the number of clusters. In-
97
dex values for totals of two to eight
k
-means clusters are depicted in Figure S1 and in-
98
dicate an optimal number of four. Clustering was accomplished using MATLAB
®
’s kmeans
99
function with a squared Euclidian distance metric and 1000 maximum iterations. Clus-
100
ters were numbered such that the cluster with the most sites was designated as Cluster
101
1 and progressively less populated clusters were assigned chronologically increasing num-
102
bers. Use of other bottom boundary depths does not dramatically affect the results, as
103
the clusters remain similar in terms of optimal number (3-6) and general character across
104
the explored bottom boundary range (10 to 400 mbsf in 10 m increments).
105
Following clustering, water depth, porosity, total organic carbon (TOC) content,
and calcium carbonate (CaCO
3
) content data from the deep ocean drilling databases were
extracted and examined to understand the controls on the [SO
2–
4
] profiles across all sites.
Sedimentation rates at each site were estimated based on water depth using an empir-
ical relationship (Middelburg et al., 1997) as follows:
ω
= 10
a
+
bZ
∗
CF
(4)
where
ω
is the sedimentation rate in cm/yr,
a
= -0.87478367,
b
= -0.00043512,
Z
is the
water depth, and
CF
= 3.3. Mean aerial nSRRs at all sites were also estimated using
the following equation (Canfield, 1991):
Mean aerial nSRR =
φD
s
dC
dz
+
φ
o
ωC
o
−
φ
b
ωC
b
(5)
Here,
φ
is the mean porosity within the zone of SO
2–
4
reduction,
D
s
is the sedimentary
diffusion coefficient for SO
2–
4
,
dC
dz
is the SO
2–
4
concentration gradient with depth,
φ
o
is
the initial phi at the sediment-water interface,
ω
is the sedimentation rate,
C
o
is the ini-
tial SO
2–
4
concentration at the sediment-water interface,
φ
b
is the porosity at the base
of the zone of SO
2–
4
reduction, and
C
b
is the minimum SO
2–
4
concentration (i.e., the SO
2–
4
con-
centration at the base of the zone of SO
2–
4
reduction). Porosities at various depths were
calculated using a one- or two-term exponential fit to site porosity data with depth. SO
2–
4
con-
centration gradients were estimated as follows:
dC
dz
=