of 15
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
=
[SO
2–
4
]
o
[SO
2–
4
]
b

z
b
(6)
where [SO
2–
4
]
o
is the SO
2–
4
concentration at the sediment water interface, [SO
2–
4
]
b
is the
minimum SO
2–
4
concentration, and
z
b
is the depth of the minimum SO
2–
4
concentration.
The SO
2–
4
sedimentary diffusion coefficient was estimated assuming a constant molec-
ular diffusion coefficient of 4.88 x 10
6
cm
2
/s for SO
2–
4
as follows:
D
s
=
D
o
φ
1
n
(7)
where
D
o
is the molecular diffusion coefficient for SO
2–
4
and
n
= 1.8 (Berner, 1980).
106
Absolute nSRRs in mol/m
3
/yr at different sediment depths were estimated using
the following equation:
nSRR
(
z
) =
1
φ
"
D
s
d
(
φ
dC
dz
)
dz
+
d
(
φω
c
C
z
)
dz
#
(8)
Here,
ω
c
is the compaction-corrected solid phase burial rate at depth
z
and
C
z
is the cor-
responding SO
2–
4
concentration. In this case,
φ
was estimated based on exponential fits
–4–
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 shipboard measured porosity data with depth and
dC
dz
was estimated by differentiat-
ing a spline fit to the shipboard [SO
2–
4
] data.
ω
c
was calculated as follows:
ω
c
=
ω

1
φ
f
1
φ

(9)
where
φ
f
is the porosity value at the bottom depth in the exponential fit to the poros-
107
ity data.
108
For estimating burial fluxes, initial porosity values for all DSDP, ODP, and IODP
109
sites with data at depths
<
5 mbsf were plotted and fitted with a linear regression against
110
water depth after removing outliers. Initial porosities across the global ocean were then
111
estimated by using this regression and the resulting values applied to ETOPO1 hypsom-
112
etry grid points (Amante & Eakins, 2009) binned into 1 m depth bins. Sedimentation
113
rates were similarly estimated by applying Equation 4 across the same depth increments.
114
Surface areas represented by the ETOPO1 grid cells were calculated based on the lat-
115
itude and longitude of grid cell boundaries (Kelly &
ˇ
Savriˇc, 2021) and used to calculate
116
average initial porosities and sedimentation rates within the shelf, slope, rise, and abyss
117
depth classes; these averages were weighted by the estimated surface areas included within
118
the associated 1 m depth bins.
119
2.2 Pyrite Compilation
120
To understand relationships among pyrite burial flux, pyrite
δ
34
S, and water depth,
121
we compiled literature measurements of pyrite abundance and
δ
34
S from marine sedi-
122
ments and imported the data into MATLAB
®
. Pyrite ages were estimated using pub-
123
lished age models or by combining sediment depth information with sedimentation rates
124
approximated from water depth (Middelburg et al., 1997). Abundance outliers and sam-
125
ples of
>
40 Ma estimated age were removed from the compilation to reduce the influ-
126
ence of unrepresentative samples (e.g. pyrite nodules) and samples unlikely to have been
127
deposited at water depths similar to modern. Pyrite S accumulation rates were estimated
128
using pyrite abundance data, sedimentation rate estimates (Middelburg et al., 1997), and
129
a linear fit of initial porosity data (from the JANUS and LIMS databases) with water
130
depth. When unreported, water depths were in rare cases estimated based on published
131
GPS coordinates or from locality maps within the corresponding studies.
132
3 Results
133
3.1
k
-means clustering
134
Profiles for the optimal number of clusters (four) at a bottom boundary depth of
135
100 meters below seafloor (mbsf) are displayed in Figure 1 and are notably different. Pro-
136
files in the most populated cluster (Cluster 1;
n
= 352) show little decrease in [SO
2–
4
]
137
with depth, consist with mostly oxic sedimentary conditions. In contrast, profiles in the
138
second most populated cluster (Cluster 2;
n
= 219) feature relatively quick SO
2–
4
deple-
139
tion with depth, consistent with anaerobic oxidation of methane (AOM) in the upper
140
part of the sediment. The centroid (i.e., within-cluster average) for this cluster shows
141
that the average depth of complete SO
2–
4
consumption is
30 mbsf. Cluster 3 (
n
= 153)
142
is intermediate in character between the two most populated clusters; profiles feature
143
a slight [SO
2–
4
] decrease with depth and 100 mbsf [SO
2–
4
] typically well above zero, con-
144
sistent with organoclastic MSR dominating the [SO
2–
4
] profiles. Profiles within the fi-
145
nal cluster (Cluster 4;
n
= 5) are anomalous and show an increase in [SO
2–
4
] suggesting
146
a source of SO
2–
4
at depth within the sediments (B ̈ottcher et al., 1998). As indicated by
147
the slopes of the [SO
2–
4
] gradients, the mean nSRR within the clusters decreases such
148
that Cluster 2
>
Cluster 3
>
Cluster 1.
149
–5–
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
0
10
20
30
40
[
S
O
4
2
-
]
(
m
M)
[
S
O
4
2
-
]
(
m
M)
0
10
20
30
40
50
60
7
0
80
90
100
Depth (mbsf)
C
l
u
s
t
e
r 1
(n = 352)
0
10
20
30
40
0
10
20
30
40
0
10
20
30
40
50
60
7
0
80
90
100
Depth (mbsf)
0
20
40
60
80
C
l
uste
r 2
(n = 219)
C
l
uste
r 3
(n = 153)
C
l
uste
r 4
(n = 5)
1
8
0
°
1
2
0
°
W
6
0
°
W
0
°
6
0
°
E
1
2
0
°
E
1
8
0
°
80
°
S
40
°
S
0
°
40
°
N
80
°
N
A
B
Figure 1.
(A) [SO
2–
4
] depth profiles for 0 to 100 mbsf for sites within Clusters 1 to 4. Dashed
lines denote the centroid for each cluster. (B) Global map of sites included within the cluster
analysis. Symbol shape and color denote cluster assignment as indicated in (A). Map created
using the M
Map mapping package (Pawlowicz, 2020).
Consideration of the geographic distribution of cluster sites (Figure 1B) yields fur-
150
ther insight into cluster characteristics. While the relatively high nSRR sites of Clus-
151
ters 2 and 3 are located predominantly in coastal regions, the low nSRR sites of Clus-
152
ter 1 are more commonly found in the open ocean (Figure 1B). This preference is reflected
153
in the median water depths for the sites within each cluster; while the median water depth
154
for Cluster 1 sites is
≥∼
3000 meters below sea level (mbsl), the median depths for all
155
other clusters are
≤∼
2500 mbsl. Cluster 4 sites are limited in their geographic distri-
156
bution and are not depicted in this figure due to the low number of sites. The geographic
157
disparity among the clusters is similarly mirrored by geochemical disparity, as demon-
158
strated by differences in mean TOC and CaCO
3
contents among the sites (Figure S2).
159
3.2 Pyrite compilation
160
The relationship between pyrite S burial flux and pyrite
δ
34
S resulting from our
161
literature compilation is shown in Figure 2. Pyrite
δ
34
S values vary by over 70
and
162
estimated pyrite S accumulation rates range over six orders of magnitude across differ-
163
ent marine environments. There is also some coherent structure to the variation; despite
164
substantial scatter, the data show a trend toward lower
δ
34
S as accumulation rate de-
165
creases from the maximum plotted accumulation rate toward a rate of
10
3
mol/m
2
/yr.
166
While data at accumulation rates below 10
3
mol/m
2
/yr are sparse, the apparent re-
167
versal of this trend at lower accumulation rates (Figure S3) could reflect rapid exhaus-
168
tion of labile organic substrate and preservation of pyrite exclusively formed under more
169
“closed-system” conditions (Sim et al., 2017).
170
Additional insight into the importance of different depositional environments can
171
be gleaned from the relationship between water depth and
∆δ
34
S, the difference between
172
the
δ
34
S of seawater SO
2–
4
and that of pyrite. Mean pyrite
δ
34
S weighted as a function
173
of pyrite S accumulation rate increases with water depth, while simultaneously, mean
174
sedimentary S concentrations for the same samples decrease with water depth (Figure
175
3). However, these data oversample productive coastal regions (Figure S3) with high S
176
concentrations, as evidenced by an expansion of geographic coverage and reductions in
177
mean S concentrations across nearly all water depth intervals when 1200+ additional mea-
178
surements without isotopic data from PANGAEA are considered (Figure 3 & S4). De-
179
spite this bias, we suspect that the increase in
∆δ
34
S with water depth would remain in-
180
–6–
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
10
-3
10
-2
10
-1
10
0
PyriteSAccumulationRate(mol/m
2
/yr)
-60
-40
-20
0
20
Pyriteδ
³⁴
S(‰)
1000
2000
3000
4000
5000
Water
Depth
(mbsl)
10
1
Figure 2.
Literature compilation of 343 measurements (colored circles) of pyrite
δ
34
S plotted
against estimated pyrite S accumulation rate. Color of symbols denotes modern water depth and
solid black line denotes a 30-point moving average. Accumulation rates were estimated based
on measured pyrite S content and sedimentation rate; rates ¡ 10
3
mol/m
2
/yr (n = 32) and
¿ 10
1
mol/m
2
/yr (n = 1) are not shown due to low data density (¡ 20 per order of magnitude
rate change). Sedimentation rates were taken from the literature or were estimated from water
depth based on the relationship of (Middelburg et al., 1997). Data from (Alt & Burdett, 1992;
Puchelt & Hubberten, 1980; Zak et al., 1980; Sweeney & Kaplan, 1980; Goldhaber & Kaplan,
1974; Scharff, 1980; Bonnell & Anderson, 1987; B ̈ottcher et al., 2004; A. Migdisov et al., 1980;
A. A. Migdisov et al., 1983; Lein et al., 1976; Krouse et al., 1977; Johnston et al., 2008; Vino-
gradov et al., 1962; Kaplan et al., 1963; Sweeney, 1972; Thode et al., 1960; Kaplan et al., 1960;
Lein, 1983; Hartmann & Nielsen, 2012; A. A. Migdisov et al., 1974).
tact in a geographically unbiased data set given our weighting of the mean pyrite
δ
34
S
181
by S accumulation rate (which accounts for the larger influence of sediments with high
182
S concentrations on the
δ
34
S of the overall pyrite burial flux within each depth interval).
183
4 Discussion & Conclusions
184
k
-means clustering identifies S cycle influences, but accurate prediction
185
of nSRR at unsampled sites remains difficult
186
Our application of
k
-means clustering to deep sea drilling data has demonstrated
187
its effectiveness at finding trends within large data sets in a less computationally com-
188
plex manner than in previous work (Bowles et al., 2014). The geographic locations (Fig-
189
ure 1B) and geochemical statistics associated with the clusters (e.g. Figure S2) also un-
190
derscore the strong influence of terrigenous input on MSR at a given site. Clusters with
191
higher nSRRs are associated with higher TOC content, lower CaCO
3
content, and closer
192
proximity to land than those with relatively low nSRRs. We note there is much less over-
193
lap among the clusters for TOC content than for CaCO
3
content and water depth, sug-
194
gesting that TOC input is the primary control on the [SO
2–
4
] profile at a given site.
195
The disparate mean TOC contents among the clusters (Figure S2) suggests some
196
potential for the nSRR at a given site to be estimated from TOC data. However, plots
197
of mean TOC content against water depth (Figure S5) and mean nSRR against mean
198
TOC burial flux (Figure S6) show enough scatter to render such estimates highly im-
199
precise; actual values routinely deviate by an order of magnitude or more from the best
200
fit relationships to these data sets. Estimates could likely be improved through substi-
201
tution of independently-measured sediment accumulation rates the addition of other re-
202
active transport variables affecting sedimentary sulfur cycling (e.g. temperature, sedi-
203
mentary iron content, and porosity) as independent variables.
204
–7–
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
Δδ
³⁴
S (‰)
Water Depth (mbsl)
60
50
40
30
20
10
0
<
50
> 50 &
<
300
> 300 &
<
1000
> 1000 &
<
2000
> 2000 &
<
3000
> 3000 &
<
4000
> 4000
0
0.5
1.0
1.5
S (wt %)
Figure 3.
Average
∆δ
34
S between seawater SO
2–
4
and pyrite (blue bars, left axis), average
solid phase S concentration for samples with S isotopic measurements (magenta bars, right axis),
and average solid phase S concentration for data including 1200+ samples without isotopic data
(light pink bars, right axis) binned as a function of water depth. Error bars denote the 2
σ
stan-
dard error of the mean.
Mass balance calculations confirm a biased or incomplete assessment of
205
the marine S cycle
206
Uncertainties in nSRR and other deep marine sedimentary parameters have im-
207
portant consequences for our understanding of isotope mass balance within the modern
208
marine S cycle. To demonstrate this, we constructed a modified steady state isotope mass
209
balance equation as follows:
210
δ
34
S
in
=
δ
34
S
evap
(1
f
py,tot
)
(
∆δ
34
S
py,shelf
f
py,shelf
+
∆δ
34
S
py,slope
f
py,slope
+
∆δ
34
S
py,rise
f
py,rise
+
∆δ
34
S
py,abyss
f
py,abyss
) (10)
where
δ
34
S
in
is the
δ
34
S of the weathering input,
δ
34
S
evap
is the
δ
34
S of buried sulfate
evaporites (assumed to equal
δ
34
S
sea
, the
δ
34
S of seawater here),
f
py,tot
is the total frac-
tion of the input flux buried as pyrite, and the remaining terms are the
∆δ
34
S and pyrite
burial fractions associated with shelf, slope, rise, and abyss sediments. We can also con-
struct a simple sum of pyrite burial fluxes from individual depositional environment classes
to obtain a total pyrite burial flux estimate:
F
py,tot
=
F
py,shelf
+
F
py,slope
+
F
py,rise
+
F
py,abyss
(11)
Using the averages from the S concentration compilation displayed in Figure 3, we as-
211
signed characteristic solid phase S concentrations (0
.
231
±
0
.
022, 0
.
250
±
0
.
018, 0
.
204
±
212
0
.
017, and 0
.
215
±
0
.
031 wt%; 2
σ
SE) and depth intervals (
<
300, 300-2000, 2000-4000,
213
and
>
4000 mbsl) to shelf, slope, rise, and abyss depositional environments, respectively.
214
We then used the ETOPO1 version of Earth’s modern hyposometry (Amante & Eakins,
215
2009) to estimate the modern pyrite S burial flux out of Earth’s oceans. After calculat-
216
ing initial porosities for each ETOPO1 depth bin based on water depth, assuming that
217
–8–
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
50
±
5% of shelf area constitutes permeable sands with no pyrite burial (Huettel et al.,
218
2014), and incorporating the root mean squared errors of the regressions of sedimenta-
219
tion rate and initial porosity against water depth, we arrived at an estimated total pyrite
220
S burial flux of 2
.
7
±
1
.
3 (1
σ
) Tmol S/yr. The large uncertainty in this estimate orig-
221
inates almost exclusively from the sedimentation rate regression (Middelburg et al., 1997)
222
and could be reduced through improved prediction of sedimentation rate from water depth.
223
Irregardless, this estimate is within error of the riverine SO
2–
4
input to seawater (2
.
8
±
224
0
.
4 Tmol S/yr (Burke et al., 2018)) and is compatible with a marine S cycle in steady
225
state. However, it is substantially higher than a prior estimate (1.2 Tmol S/yr; (Berner,
226
1982)). This disagreement and the implication that nearly all (
96%) of the riverine in-
227
put is being buried as reduced S suggests that biases remain. Comparison of the S burial
228
fluxes with aerial nSRRs derived from the
k
-means clusters further supports remaining
229
biases (Supporting Information).
230
One potential source of systematic bias in solid phase estimates of the reduced S
231
burial flux is sediment porosity. Although sediment porosity data are commonly collected
232
in deep sea drilling efforts, disturbance during the coring process can easily lead to loss
233
or compaction of core material (Morton & White, 1997). Such processes would cause the
234
in-situ porosity of recovered samples to be underestimated by post-coring measurements
235
and the corresponding solid phase-based burial estimates to be overestimated. Apply-
236
ing a linear regression between initial porosity and water depth that instead considers
237
the maximum initial porosity in a 50-point moving window with water depth (Figure S7)
238
lowers the pyrite burial flux estimate to 0
.
87
±
0
.
41 Tmol S/yr. This burial flux is now
239
within error of the prior estimate of 1.2 Tmol S/yr (Berner, 1982). While an inaccurate
240
estimate of the hydrothermal flux remains possible, this finding and preliminary com-
241
parison of deep sea drilling porosity data with similar multi-core data (Figure S7) strongly
242
suggests biased porosity measurements contribute to the high burial flux estimate. An
243
increase in dissolved-phase output flux estimates due to such a porosity bias suggests sam-
244
pling and analytical biases better explain the high global nSRR estimate of (Bowles et
245
al., 2014) (Supporting Information).
246
Assignment of characteristic
∆δ
34
S values to the aforementioned depth intervals
247
allows for further assessment of fluxes and isotope mass balance in the marine S cycle.
248
Assuming characteristic pyrite
δ
34
S values of
15
.
2
,
25
.
3
,
26
.
1
, and
23
.
1
and
249
f
py,i
values of 0.46, 0.35, 0.155, and 0.035 for the shelf, slope, rise, and abyss (respectively),
250
equation 10 yields a steady state seawater SO
2–
4
δ
34
S value of +17.3
assuming a pyrite
251
burial flux of 0.87 Tmol S/yr (
f
py,tot
= 0.31) and +28.7
if a pyrite burial flux of 1.3
252
Tmol S/yr (
f
py,tot
0
.
45, the maximum value considered likely in a prior study (Tostevin
253
et al., 2014)) is imposed; the sulfate evaporite burial flux is assumed equal to the pre-
254
anthropogenic riverine input (2.8 Tmol S/yr,
δ
34
S
in
= +4.6
(Burke et al., 2018)) mi-
255
nus the pyrite burial flux in each case. The bracketing of the measured modern seawa-
256
ter
δ
34
S of +21.24
(Tostevin et al., 2014) by these values is encouraging and suggests
257
general accuracy in both the pyrite burial flux magnitudes in each depositional environ-
258
ment and their isotopic composition. If
δ
34
S
in
remains unchanged, equation 10 can yield
259
a seawater
δ
34
S value in agreement with the observed value if (1)
f
py,tot
0
.
37 (pyrite
260
burial flux = 1.03 Tmol S/yr) and the pyrite
δ
34
S values are kept the same; (2) one or
261
more of the
∆δ
34
S values is increased while keeping
f
py,tot
0
.
31 (pyrite burial flux =
262
0.87 Tmol S/yr); (3) one or more of the
∆δ
34
S values is decreased while keeping
f
py,tot
263
0
.
45 (pyrite burial flux = 1.3 Tmol S/yr); or some intermediate combination of 1-3. We
264
favor scenario (3) given the absence of deltaic sediments with closed-system S diagen-
265
esis and high pyrite
δ
34
S in our
δ
34
S compilation and the relatively close match between
266
the corresponding pyrite burial flux and Berner’s (1982) independently-derived value (1.2
267
Tmol S/yr). Increasing the shelf pyrite
δ
34
S by about 17
to +2.0
(i.e., lowering the
268
total
∆δ
34
S to +34.0
) is sufficient to match the observed seawater
δ
34
S in this case.
269
–9–
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
15
20
25
A
δ
34
S
sea
(‰)
0
5
10
10
1
1
B
Shelf Pyrite Burial
Slope Pyrite Burial
Pyrite Burial
Flux (mol S/yr)
-
100
-
50
0
50
100
R
e
l
a
t
i
ve
S
e
a
L
eve
l
C
h
a
n
g
e
(
m
)
0
.
5
1
1
.
5
C
Relative
Shelf Area
Figure 4.
(A) Calculated steady-state seawater
δ
34
S for shelf
∆δ
34
S = 19.2 and slope
∆δ
34
S
= 46.5
. (B) Continental shelf (red), continental slope (blue), and total (black dashed) pyrite
burial fluxes. (C) Calculated change in the amount of submerged shelf area (0-300 mbsl) relative
to modern for a given change in sea level.
Seawater
δ
34
S is sensitive to pyrite burial locii changes
270
In the mass balance scenarios considered here, nearly half of pyrite burial in the
271
modern ocean occurs in shelf sediments. However, prior work (Berner, 1982) has discussed
272
a possible transfer of some shelf sedimentation to the deep ocean in conjunction with lower
273
sea levels during the last glacial period. Burial of pyrite in proportion to sedimentation
274
would concomitantly shift some shelf pyrite burial to deeper environments in this sce-
275
nario. How would such a transfer affect S isotope mass balance?
276
To estimate the sensitivity of seawater
δ
34
S to sedimentation shifts, we used the
277
ETOPO1 arc-minute global relief model (Amante & Eakins, 2009) to determine the ap-
278
proximate change in shelf area for changes in sea level (SL) ranging from -100 m to +100
279
–10–
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.