Research ArticleGEOCHEMISTRY

Strong local, not global, controls on marine pyrite sulfur isotopes

See allHide authors and affiliations

Science Advances  26 Feb 2021:
Vol. 7, no. 9, eabb7403
DOI: 10.1126/sciadv.abb7403


Understanding variation in the sulfur isotopic composition of sedimentary pyrite (δ34Spyr) is motivated by the key role of sulfur biogeochemistry in regulating Earth’s surface oxidation state. Until recently, the impact of local depositional conditions on δ34Spyr has remained underappreciated, and stratigraphic variations in δ34Spyr were interpreted mostly to reflect global changes in biogeochemical cycling. We present two coeval δ34Spyr records from shelf and basin settings in a single sedimentary system. Despite their proximity and contemporaneous deposition, these two records preserve radically different geochemical signals. Swings of ~65‰ in shelf δ34Spyr track short-term variations in local sedimentation and are completely absent from the abyssal record. In contrast, a long-term ~30‰ decrease in abyssal δ34Spyr reflects regional changes in ocean circulation and/or sustained pyrite formation. These results highlight strong local controls on δ34Spyr, calling for reevaluation of the current practice of using δ34Spyr stratigraphic variations to infer global changes in Earth’s surface environment.


In marine sediments, fractionation of sulfur isotopes during microbial sulfate reduction (MSR) leads to the production of 34S-depleted sulfide in the porewater environment, some of which may react with dissolved ferrous iron (FeII) or with ferric iron (FeIII)–bearing solids to form pyrite [FeS2; (1)]. The formation and subsequent burial of pyrite is the main sink of reduced sulfur from the ocean and, together with organic carbon burial, contributes to the maintenance of an oxidized Earth surface (2). As such, variation in sulfur isotope ratios of sedimentary sulfate and pyrite, as well as the offset between δ34S values of pyrite and coeval seawater sulfate (Δpyr ≡ δ34Ssw – δ34Spyr), is widely used to reconstruct the ancient biogeochemical cycles of sulfur, carbon, and oxygen (312). Within this framework, stratigraphic variations in δ34Spyr and/or Δpyr are often interpreted to reflect temporal changes in the global sulfur cycle, such as periods of enhanced pyrite burial [e.g., (6, 9)] and/or changes in the size of the marine sulfate reservoir [e.g., (4, 5)]. However, recent observations of large Δpyr variability on temporal or spatial scales too small to reflect changes in the global sulfur cycle raise concerns regarding the ability of both sulfate-bearing phases and pyrite to preserve a global sulfur signal (1319).

Sulfate is well mixed in the modern ocean, and its sulfur isotopic composition does not vary among or within oceanic basins (12). Consequently, variation in Δpyr reflects the microbial and physical processes of transport and mineralization preserved in δ34Spyr (1420). On average, δ34Spyr is higher in shallower depositional environments than in deep ocean sediments, and this difference is thought to reflect shallow-deep differences in the relative importance of diffusion and downward advection of sediment (hereafter referred to as burial) as sources of porewater sulfate (13). Proximal settings are typically characterized by more rapid sedimentation and high rates of organic matter and reactive iron delivery. In such settings, diffusive time scales are long relative to the rate of burial, and porewater sulfate is effectively isolated from the overlying water-column sulfate reservoir at a relatively shallow depth within the sediment. A rapid MSR-driven decrease in the concentration of porewater sulfate and the associated isotopic discrimination results in down-sediment enrichment in 34S of the residual sulfate and instantaneous sulfide. Furthermore, the relatively rapid delivery of reactive iron results in effective sequestration of MSR-produced sulfide in iron sulfide minerals. Such effects in steadily depositing shelf environments are convincingly demonstrated by large and systematic δ34Spyr variations at the temporal scale of glacial-interglacial cycles (i.e., 100 ka), which were forced by sea level–driven changes in sedimentation rate (16).

In addition to the effects described above, proximal settings are more prone to short-term, stochastic variations under depositional conditions, where episodic resuspension and/or reorganization of the benthic ecosystem promote reoxidation of the local sulfide pool (2, 21, 22). A prominent example is the Guyana mobile mud belt, where frequent remobilization of surface sediment in response to energetic currents is responsible for intense cycling of redox conditions, resulting in unusually 34S-enriched pyrite [δ34S ranging from 0‰ to >30‰ Vienna Canyon Diablo Troilite (VCDT); (23)]. In contrast with the above, fine-grained offshore sediments are usually deposited more slowly in a relatively quiescent environment, where diffusion time scales are shorter than the time scale for burial. In such depositional environments, effective communication with the overlying water column is expected to lead to a gradual decrease (with depth) in the concentration of porewater sulfate and a correspondingly gradual increase in the δ34S of porewater sulfate and the sulfide (and pyrite) produced from it by MSR.

Burrowing and mixing by sediment-dwelling fauna (bioirrigation and bioturbation) ventilate the sediment to depths that would, in the absence of this activity, communicate with overlying seawater on diffusive time scales that are long relative to those of bioventilation. In so doing, bioventilation affects the balance of diffusion and burial as agents of porewater sulfate transport and oxygenation [e.g., (2)]. Conventionally, in bioventilated sediment, more efficient buffering of the porewater sulfate concentration and sulfur isotopic composition to those of seawater sulfate is expected, leading to preservation of 34S-depleted pyrite (21). However, below the bioventilated zone, a dependence of the diffusion scale on the sedimentation rate is again expected, and so shelf-basin differences in sedimentation rate and continuity are expected to give rise to shelf-basin differences in Δpyr. As bioventilation depth and activity do not differ statistically between shelf and deep environments (24), such Δpyr differences, when observed, are unlikely to be the outcome of differing degrees of bioturbation and bioirrigation.

Despite the anticipation of shelf-basin decoupling of pyrite sulfur isotopes, it has proven hard to observe this relationship over geologic time scales, mostly because of the absence of a long marine archive from which the dynamics in both onshore and offshore settings could be studied simultaneously. Here, we combine δ34Spyr records from shelf and basin settings in a single sedimentary system with a robust age model to clearly show the persistence of a regionally divergent sulfur isotope signal over the past 40 Ma. Coeval records of seawater sulfate δ34S values over the study period (11) allow calculation of Δpyr values, which also diverge between the shelf and basin records. Such divergent stratigraphic variations between two sites separated by only ~1000 km suggest local controls on δ34Spyr and Δpyr and call for a reevaluation of current practices of using variations in the sulfur isotope composition of sedimentary phases to infer substantial changes to Earth’s surface environment.


Depositional settings

The Eastern New Zealand Oceanic Sedimentary System (ENZOSS) is a large sediment recycling system, which extends over the eastern, 4500-km-long margin of the New Zealand microcontinent, a region that is currently sheared by the active collision between the Australian and Pacific plates [Fig. 1A ; (25, 26)]. Uplift and erosion produce copious fine-grained terrigenous material that is transferred to the coast [up to 109 metric tons per year of fluvial sediment discharge; (27)]. Some of this sediment is eventually captured within the shelf sedimentary system or guided through three submarine channels to the deep Southwest Pacific Basin (2830). There, the Deep Western Boundary Current (DWBC), initially generated by the Antarctic Circumpolar Current (ACC), redirects the material to a series of sediment drifts (2730).

Fig. 1 Location of studied boreholes.

(A) IODP 1352: 44.9374°S; 172.022692°W and ODP 1123: 41.786°S; 171.499°W. Also shown are the bathymetry (depth contours in meters) and the deep Pacific inflows; the ACC and the DWBC. Redrawn from (27). (B) Simplified stratigraphy of studied boreholes with corresponding lithological units. Age-calibrated stratigraphic correlations from this study are shown in Ma. MP, regional Marshall Paraconformity; E, Eocene.

We constrain the behavior of sedimentary pyrite sulfur isotopes preserved onshore on the shelf and offshore in the deep basin by using sediments from boreholes International Ocean Discovery Program (IODP) 1352 (44.9374°S, 172.022692°W, 343.8-m water depth) and Ocean Drilling Program (ODP) 1123 (41.786°S, 171.499°W, 3290-m water depth), respectively. The geochronology of the investigated sediment is strongly constrained by high-resolution δ18O measurements on benthic foraminifer tests (3133). The age model used here was obtained by tuning the δ18O data in the studied boreholes to the global δ18O stack over the past 1.9 Ma [see Materials and Methods and fig. S1; (34)]. The age of older sediments was constrained by biostratigraphic and/or magnetostratigraphic tie points [fig. S1; (35, 36)], which, while less robust, is suitable for the purposes of the current study.

Sediments drilled in the IODP 1352 borehole can be divided into three distinct sedimentological units (Fig. 1B). Unit I spans the Pleistocene and consists of a 520-m-thick mud-rich sediment with interbedded coarser beds of muddy fine sand layers reflecting cyclic glacial-interglacial sedimentation. Indications of bioirrigation are mostly preserved in the muddy glacial layers and often intensify toward the uncemented firm grounds (37). The transition between units I and II is gradual, reflecting a progressive change in water depth to a deeper slope setting. Unit II contains a gradual change from uncemented to lithified marlstone with less frequent sandstone layers, except at the base of the unit (Miocene), where glauconitic sandstone becomes more common. An erosional surface, known as the Marshall Paraconformity, occurs at the base of unit II [1750 mbsf (meters below seafloor)]. Below, unit III is composed of nannofossil micritic limestone (hemipelagic to pelagic foraminifera) of Eocene age (Fig. 1B).

The offshore ODP 1123 borehole recovered a 633-m-long succession of clay-rich nannofossil ooze, chalk, and limestone (Fig. 1B). These sediments can be divided into four distinct sedimentological units. Unit I preserves Pleistocene glacial-interglacial couplets. Unit II presents the characteristic sedimentary features of a drift deposit, such as cycles of coarsening-upward and fining-upward grain size, irregular lenses, and cross-laminated alignments (38). Unit III is composed of sediments similar to those in the overlying unit II, with a general increasing thickness of the silty-sandy layers and the terrigenous clay content with depth in the unit. This unit is classified as a nannofossil mudstone. Below the Marshall Paraconformity, the remainder of the core is composed of unit IV, an Eocene micritic limestone that appears correlative with the bottom of the shelf IODP 1352 core [Fig. 1B; (39)]. More detailed sedimentary information is given in Materials and Methods.

Onshore-offshore δ34Spyr

We performed 185 sulfur isotope analyses of chromium-reducible sulfur species (pyrite and other forms of reduced sulfur, notably iron monosulfides) along the two boreholes (see Materials and Methods and table S1). Extractions of acid-volatile sulfur (mostly iron monosulfide) from a subset of samples yielded no measurable sulfur (fig. S2 and table S2), suggesting that the chromium-reducible sulfur was almost entirely pyrite. The sulfur isotopic composition of this pyrite is reported as Δpyr, the offset (in per mil) between our δ34Spyr and extrapolated coeval δ34Ssw (fig. S3; Materials and Methods). Shelf and basin Δpyr values exhibit different long-term trends and short-term variability over the past 40 Ma (Fig. 2). The shelf site (IODP 1352) shows extreme variation in Δpyr, from 3.1 to 67.7‰, whereas its pyrite content varies between 0.01 and 0.97 weight % (wt %), with no appreciable long-term trend or correlation with Δpyr (figs. S4 and S5). The smallest Δpyr values of 3.1 and 3.8‰ occur at 0.44 and 5.75 Ma during intervals of global benthic δ18O positive excursions (34, 39, 40). Over the past 1.9 Ma, Δpyr values covary with climate and sea level, as deduced from the age model. Interglacial periods of high sea level are characterized by large Δpyr values and muted isotopic variability [average Δpyr = 54.2 ± 7.0‰ (1σ), n = 24], whereas glacial intervals of lower sea level display smaller and more variable Δpyr values [average Δpyr = 25.5 ± 10.9‰ (1σ), n = 18; fig. S6]. This correspondence of shelf Δpyr values with the amplitude and dynamics of global sea-level fluctuations (34, 39, 40) extends over the past 40 Ma (Fig. 2).

Fig. 2 Divergent stratigraphic Δpyr records.

Shelf and deep basin Δpyr values (red and blue, respectively, this study) over the past 40 Ma. Filled markers have strong stratigraphic age constraints (foraminiferal δ18O measurements), whereas empty markers have less robust age constraints (biostratigraphic and/or magnetostratigraphic tie points). For global climatic context, deep-sea records of δ18O [light gray circles and dark line; (34, 40)] are shown. Note the change of scale of δ18O associated with the change of time scale between the middle and rightmost parts of the figure, as represented with the inner gray y axis on the right.

In contrast, basinal Δpyr values display a mean of 72.1 ± 4.8‰ (1σ), with the standard deviation being controlled by a secular increase from ~53 to ~73‰ over the past 40 Ma rather than by shorter-term variation (Fig. 2). Over the entire duration of the record, no correlation is observed between sea-level and basinal Δpyr values (Fig. 3A). The interval 0.9 to 0.02 Ma marks the onset of a long-term plateau of larger Δpyr values (73.5 ± 1.6‰ on average). The long-term monotonic increase in Δpyr values is not accompanied by any observable change in pyrite content (4 ± 1 × 10−3 wt % from chromium-reduction extractions, 7 ± 4 × 10−2 wt % after accounting for post-sampling pyrite oxidation, see Materials and Methods) until ~5 Ma, when an increase in pyrite abundance up to a maximum of 0.6 uncorrected wt % occurs at approximately the same time that Δpyr reaches ~73‰ (figs. S2, S3, and S7).

Fig. 3 Strong sea-level control.

(A) Relationship between change in sea level (m relative to present day) and Δpyr values. The green band represents the theoretical equilibrium sulfate-sulfide sulfur isotope fractionation (ΔSO4-H2S) at the bottom-water temperature of the study sites (43). Filled markers have strong stratigraphic age constraints (foraminiferal δ18O measurements), whereas empty markers have less robust age constraints (biostratigraphic and/or magnetostratigraphic tie points). (B) Box and whisker plots of Δpyr values grouped by systems tract associated with fourth-order sea-level variation. Red horizontal lines denote the median; the box upper and lower boundaries denote the 75th and 25th percentiles of the data, respectively; and the black horizontal lines contain 99.3% of the data. Data include only samples with strong stratigraphic age constraints (i.e., 0 to 1.9 Ma). HST, highstand systems tract; LST, lowstand systems tract; TST, transgressive systems tract.

From these observations alone, reconstructing a coherent stratigraphic onshore-offshore Δpyr framework appears difficult. Changes in the depositional environment associated with glaciations may explain the observed swings in shelf Δpyr values but fail to explain the absence of such swings alongside a longer-term secular Δpyr increase recorded at the deep basin site. This heterogeneity in the sulfur isotope composition within a single sedimentary system is remarkable, considering the common use of stratigraphic variations in δ34Spyr or Δpyr values, sometimes from only a few localities, to reconstruct global changes in pyrite burial, marine sulfate levels, or the type of microbial metabolisms present (312).


Shelf pyrite δ34S values controlled by sea level–driven changes in sedimentation

The highly variable onshore Δpyr record provides clear insight into the effect of global eustatic sea-level variations on the local control of the New Zealand shelf diagenetic system and pyrite sulfur isotopes. The mid-Miocene marks a period of high sea level known as the Miocene Climate Optimum, which was followed by a general decrease in sea level, starting ~13.9 Ma and culminating with the onset of northern hemisphere glaciation (39, 40). Over this period, uplift of the Southern Alps led to high rates of sediment input to the Canterbury shelf, initiating drift deposition of silt near the slope, until ~3.25 Ma (28). This long-term decrease in sea level was punctuated by several major increases in ice volume, accompanied by drops in global sea level, as attested by positive benthic δ18O steps or excursions (34, 39, 40). The termination of drift development may have been related to initiation of Late Pliocene to Pleistocene high amplitude sea-level variations, which could have enhanced downslope processes relative to shelf sediment redistribution by bringing fluvial sources and/or shorelines closer to the shelf edge during glacial times (2730).

We propose a direct link between sea-level drop, an increase in sedimentation rate, and/or burial of reactive organic matter and iron, leading to more complete conversion of porewater sulfate to pyrite under “burial-dominated” conditions and small Δpyr intervals. Initiation of late Pliocene to Pleistocene high-amplitude sea-level variations resulted in high-frequency migration of the sedimentary system along the shelf. By controlling the depositional properties (organic carbon and iron loading, sedimentation rate, and porosity), sea-level changes are directly responsible for the bimodal Δpyr distribution observed over the past 2 Ma in the shelf sediments. The amplitude of Δpyr oscillation increases from ~15‰ to >60‰ over the mid-Pleistocene transition (i.e., between 1.25 and 0.65 Ma), when the dominant periodicity of climate cycles changed from ~40 to ~100 ka, and the ice ages intensified and acquired their asymmetrical character [i.e., gradual glacial growth and abrupt deglaciation; (33)]. This Δpyr variation amplitude increase is caused by both a downward shift of the smallest Δpyr values and an upward shift of the largest Δpyr values, and we suggest that it is driven by higher-amplitude glacial-interglacial cycles of melting and buildup of continental ice sheets (Figs. 2 and 3A). The frequency of these Δpyr variations (~40 to 100 ka) implies a dominant local or regional control rather than changes in the global sulfur cycle, and we suggest that sea-level fluctuations and attendant sedimentary changes serve as this control on shelf Δpyr values. Similar variations in Δpyr with sea level have been observed in an onshore site in the Gulf of Lion in the Mediterranean Sea (16) and in an inner shelf site in the East China Sea (15). The occurrence of sea level–Δpyr covariation in three basins so distant from each other and so different in their sedimentological characteristics suggests that such covariation may be general.

We explored the possibility that changes in the oxygen penetration depth were involved in driving the observed Δpyr variations. When the oxygen penetration depth is shallow, the onset of MSR is close to the sediment-water interface, and we suggest that this translates into two opposing effects on Δpyr. On one hand, a short diffusive length scale from the water column to the zone of MSR should result in more effective buffering of porewater sulfate concentrations and isotopic compositions to those of seawater sulfate. On the other hand, a shallower onset of MSR, at depths in the sediment where organic matter is more abundant and more reactive, should lead to higher overall rates of sulfate drawdown and less effective buffering to the seawater sulfate concentration and sulfur isotopic composition. With a shallow oxygen penetration depth, less pronounced isotopic distillation of porewater sulfate due to a smaller diffusive distance should lead to larger ∆pyr, whereas more pronounced isotopic distillation due to higher overall rates of MSR should lead to smaller ∆pyr. The net effect on ∆pyr is expected to depend on the system being studied. Conversely, with a deep oxygen penetration depth, the distance for diffusion of seawater sulfate to the zone of MSR is larger, but the organic matter is older, less labile, and less abundant. Again, the expectation is of opposite effects on ∆pyr, with a system-dependent net effect.

Porewater oxygen diffuses from the overlying water column and is consumed mostly by reaction with organic matter through aerobic respiration (41, 42). Therefore, it is expected that with more rapid sedimentation or with higher total organic carbon (TOC) concentrations, oxygen is drawn down over a shallower depth in the sediment. In other words, the oxygen penetration depth is expected to decrease with increasing deposition velocity and TOC concentration or reactivity. As no systematic variation is observed in TOC concentration or isotopic composition (a proxy for its source) at the shelf site (figs. S4 and S8), it appears that variation in sedimentation rate is the main driver of variation in the oxygen penetration depth. One might then expect that glacial intervals, which are characterized by higher sedimentation rate, should be times of a shallower oxygen penetration depth at the study site. As discussed above, if the oxygen penetration depth were a major driver of variation in ∆pyr through its effect on the diffusive length scale, then the expectation would be that glacial intervals display larger ∆pyr. However, the opposite behavior is observed (smaller ∆pyr during glacial intervals; Fig. 2), suggesting that at our study site the effect of the oxygen penetration depth on ∆pyr via the length for diffusive exchange is secondary to the effects of variations in the sedimentation rate on the isotopic distillation of porewater sulfate (more rapid sedimentation, leading to more rapid disconnection of porewater from seawater sulfate).

Basin pyrite δ34S values controlled by changes in deep oceanic circulation

Over the same time period, the deep basin borehole preserves very different Δpyr variability from the onshore borehole (Fig. 2). Between ~20 and ~0.5 Ma, Δpyr in the basin sediments increases monotonically from around 53‰ to around 73‰. Over the past ~400 ka, when the shelf site displays wild swings in Δpyr, the basin site displays large and almost invariant Δpyr (73.3 ± 3.2‰). If the preserved TOC correlates with the initial organic matter content, then the relatively similar organic matter concentrations and carbon isotope compositions at both the shelf and basin sites suggest that substrate limitation and source are of minor importance in driving the observed divergent stratigraphic variations (figs. S4 and S8). We suggest that because of its distance from land (~1000 km) and water depth (~3300 m), the offshore depositional setting remained insensitive to the ~100-m variations in Pleistocene sea level and the associated variations in sediment supply (Fig. 3A). Consequently, sedimentary pyrite formation at the basin site remained diffusion-dominated, and Δpyr values remained high, relatively invariant, and close to the theoretical calculation of isotopic equilibrium between sulfate and sulfide [Fig. 3A; (43)] over glacial-interglacial sea-level cycles.

The long-term Δpyr increase (from ~50 to 70‰ between ~20 and 0.5 Ma) spans the timing of regional, tectonically driven changes in deep oceanic currents, and we suggest that these changes are a possible cause for the Δpyr increase. From the Late Cretaceous to the late Eocene, rifting created the South Pacific Ocean and marked the formation of the New Zealand Pacific continental margin (25, 26). A post-rift transgressive phase resulted in the deposition of widespread siliceous-calcareous fine-grained sediments, as reflected by our low-TOC, high-CaCO3 wt % samples (figs. S4 and S5). Deposition was interrupted by the onset of abyssal flow ~33.5 Ma (27), associated with the opening of the Tasmanian gateway and inception of a circum-Antarctic flow (precursor to the ACC). Widespread erosion, recorded in the region as the so-called Marshall Paraconformity, was followed in the late Oligocene–early Miocene (~25 Ma) by sedimentary drift accumulation under the influence of the DWBC (27).

Associated with the change in deep ocean currents is a decrease in bottom-water temperature, which is expected to have increased the temperature-dependent sulfate-sulfide equilibrium sulfur isotopic fractionation (43). If the sulfur isotopic fractionation associated with sulfate reduction is close to the equilibrium fractionation, as suggested by the low δ34Spyr values observed in the deep core, then an increase in the sulfate-sulfide equilibrium isotopic fractionation is expected to translate into a similar increase in the microbial fractionation of sulfur isotopes. The onset of delivery of cold circum-Antarctic water to the deeper parts of the ENZOSS is thought to have led to a decrease in water temperature from ~15° to ~7.5°C across the Marshall Paraconformity, followed by a gradual decrease to Pleistocene temperatures between 2.0° and −1.5°C (33, 44). The temperature-dependent sulfate-sulfide equilibrium sulfur isotope fractionation would have thus increased by about 3‰ across the Marshall Paraconformity, consistent with the ~4‰ change observed in Δpyr values over this time period (Fig. 2 and table S1). The remaining temperature decrease from ~20 Ma to the Pleistocene would cause an additional increase in the equilibrium fractionation of 3 to 4‰, perhaps explaining an increase in Δpyr values of similar magnitude but leaving ~13 to 14‰ of the observed Δpyr increase unexplained.

We suggest that the sedimentary drift accumulation, associated with the change in deep ocean circulation and supported by sedimentary structures in the studied basin borehole (Materials and Methods), may explain the remaining ~13 to 14‰ increase in Δpyr. Drift accumulation is highly dynamic, leading to nonsteady diagenetic conditions that could result in accumulation of 34S-enriched sulfide in a mobile sediment reservoir and preservation of that sulfide in pyrite. Examples of such nonsteady diagenetic conditions include, but are not limited to, chemocline migration in response to varying sedimentation rate and remobilization of sediment to expose relatively 34S-enriched sulfide to new reactive iron. Similar dynamic mechanisms have been invoked to explain the existence of 34S-enriched pyrite in shallow shelf and deltaic environments (23, 45). In the basin environment studied here, surface remobilization is less intense and more infrequent. Expected δ34Spyr values are, therefore, lower (i.e., larger Δpyr) than in the shallow environments described in previous studies, consistent with the values we observe.

The decreasing influence of reworking from the onset of redeposition after the Marshall Paraconformity (~20.5 Ma) to the Pleistocene, as suggested by the long-term monotonic increase in Δpyr, may reflect the decreasing vigor of the deep oceanic currents. Vigorous deep inflow of ACC and DWBC suggested in the Miocene (and expected also in the Oligocene) are responsible for the high mobility of the sedimentary drifts (46) and the persistence of reworking. The redirection of the ACC during the Early Pliocene and the reduction in the intensity of DWBC formation during the global Plio-Pleistocene cooling interval are jointly responsible for a substantial decline in drift activity. We suggest that this long-term decline in reworking may have caused part of the Δpyr increase between ~20 and ~0.5 Ma. The large and invariant Δpyr over the past 400 ka (δ34Spyr = 73.3 ± 3.2‰, n = 77) suggests that glacial-interglacial environmental changes did not cause substantial variation in the deep oceanic currents at our basinal site, which were too weak to cause frequent sediment remobilization and the associated 34S enrichment of porewater sulfide and pyrite.

We explored an alternative explanation that the long-term increase in Δpyr values may record the preservation in pyrite of increasingly 34S-enriched sulfide as the sediments get progressively buried. Exposure of FeS to sulfide and polysulfide species in aqueous solution results in the formation of pyrite (4750), and the isotopic composition of the pyrite formed this way is expected to reflect the δ34S values of both the FeS precursor and the porewater sulfide and/or polysulfide (51). To constrain the δ34S values of porewater sulfide, we note that the δ34S values of porewater sulfate in the deep core increase by ~35‰ over the upper 220 m of the core due to Rayleigh distillation during MSR [fig. S7; (52)]. If the microbial fractionation of sulfur isotopes remained approximately constant with depth in the sediment, then the δ34S values of the sulfide instantaneously produced by reduction of the porewater sulfate are likewise expected to increase by ~35‰ over this depth interval. Polysulfide species equilibrate rapidly with sulfide (53, 54), and their sulfur isotopic composition is also expected to increase by ~35‰ [and be offset from that of porewater sulfide by a few per mil; (51)]. The increase in δ34Spyr values observed over the upper 220 m of the core is ~15‰, approximately half the increase in porewater sulfate δ34S values (and those of instantaneously produced sulfide and/or polysulfide). Considering sulfur isotopes alone, the down-core gradients of pyrite and porewater sulfate δ34S values would suggest approximately equal proportions in the ultimate pyrite of sulfur from initial FeS (with approximately constant δ34S values) and porewater sulfide/polysulfide (with increasing down-core δ34S values). However, we observe neither an increase in pyrite abundance nor evidence for a change in the proportions of pyrite and FeS in the total chromium-reducible sulfur (i.e., no evidence for down-core transformation of FeS into pyrite; Materials and Methods). Furthermore, no systematic changes in the availability of reactive iron species are observed, possibly suggesting that pyrite formation pathways did not vary substantially over time. Conclusively quantifying the proportions of pyrite originating from different formation pathways requires information on the abundance and sulfur isotopic composition of porewater sulfide, porewater polysulfide species, solid FeS, elemental sulfur, and pyrite, not all of which are available here. However, on the basis of the available data, we suggest that progressive transformation of FeS to pyrite was a minor contributor to the observed down-core increase in δ34Spyr values.

On the basis of the available information, we suggest that the long-term, monotonic increase in Δpyr at the basinal site reflects a combination of a decrease in bottom-water temperature and a decrease in reworking intensity and frequency, both driven by changes in deep oceanic currents. In the basin, as on the shelf, changes in the sulfur isotopic composition of pyrite reflect regional oceanic, sedimentological processes and local diagenetic processes, complicating the use of δ34Spyr and/or Δpyr to infer changes in the global biogeochemical cycle of sulfur.

Rethinking spatial and temporal patterns in pyrite δ34S values

Our findings highlight the roles of the local to regional depositional environment in early diagenetic processes, and specifically in controlling Δpyr. The evidence presented here reinforces previous observations in shelf environments. Moreover, our findings indicate that deep-sea environments may be likewise prone to producing stratigraphic variation in Δpyr values that is unrelated to the global sulfur cycle, albeit because of different processes. In the shallow depositional environment, rather than recording widespread variation in the global sulfur cycle, the variations in Δpyr reflect sedimentological drivers that vary on time scales much shorter than the residence time of sulfate in seawater. In the deep depositional environment, although the time scale of variation is comparable to or longer than the residence time of sulfate in seawater, variations in Δpyr reflect regional changes in bottom currents or perhaps local diagenetic dynamics.

As these data demonstrate, a detailed understanding of the depositional context and how it affects Δpyr is essential for robust interpretations of geochemical data in both marine sediments and the geologic record of sedimentary rocks. A way forward may be a more systematic use of sedimentology and sequence stratigraphy in facies-dependent studies of Δpyr values (and other geochemical properties), which are more likely to compare similar depositional environments in a stratigraphic section. Given the ubiquity of shallow marine strata in the geologic record (55), we tested this approach by defining the transgressive-regressive sequences in the IODP borehole 1352 (Fig. 3B and Supplementary Text). Within this framework, large Δpyr values and low isotopic variability occur during highstand systems tracts, whereas transgressive and lowstand systems tracts display more variable and smaller Δpyr values. Similarly, we suggest that Δpyr excursions preserved in the geologic record could easily be compared to evidence of changing depositional environment or known variations in sea level (56, 57), especially when the duration of such Δpyr excursions is shorter than or comparable to the residence time of sulfate in seawater. In addition, the use of sequence stratigraphy allows possible correlation of stratigraphic Δpyr variations among or within sedimentary basins. Such an integrated approach may help to disentangle the multiple local to regional components that contribute to archives of Δpyr values, and their relative importance in obscuring global sulfur-cycle signals preserved in the geologic record.

Notably, throughout the time interval recorded by the studied cores, the ocean has remained well mixed, oxygenated, and sulfate rich, yet Δpyr displays large variations. There are many periods in Earth history, including major global anoxic events and mass extinctions, in which high-amplitude variations in Δpyr (and/or δ34Spyr) values have been observed and interpreted to reflect global-scale events or drivers. Close examination of these periods reveals that most of them are associated with variations in sea level (table S3), as attested by sedimentary facies changes and stratigraphic evidence. In many of these intervals, there is an apparent connection between variations in Δpyr and periods of presumed anoxia or euxinia. However, anoxia and/or euxinia are often inferred on the basis of iron speciation in sedimentary rocks or other proxies related to iron-sulfur chemistry, which may suffer from similar controls by local depositional and diagenetic processes rather than water-column chemistry (58). We therefore suggest that it is possible that many (if not most) of the geochemical perturbations in table S3, as well as others, may have been more directly related to changes in the depositional and diagenetic environment than to biogeochemical changes in the global ocean. Such local changes may be triggered globally (e.g., by a synchronous change in sea level), leading to spatially correlated shifts in Δpyr that, nevertheless, do not represent global changes in sulfur cycling or ocean chemistry. This idea is consistent with the growing number of broadly coeval carbon and sulfur isotopic records that diverge in the timing, shape, or magnitude of variability. With this in mind, the detailed consequences of variation in the properties of the depositional environment on early diagenetic processes should be examined not only for sulfur isotopes in pyrite but also for other redox-sensitive proxies closely related to iron-sulfur chemistry (e.g., the concentrations, speciation, and/or isotopic compositions of Mo, Cr, U, and Fe).



Chronological framework. Over the past 1.9 Ma, for both drill cores, we constrain the geochronology of the investigated sediment by tuning previously reported benthic δ18O measurements to the Lisiecki and Raymo (34) δ18O stack (fig. S1). These age models are, for the most part, similar to previously constructed age models at the shelf location (31, 32) and at the offshore location (33). For older sediment, where no foraminifera δ18O measurements were available, we used biostratigraphic and/or magnetostratigraphic information reported by the IODP/ODP sampling parties [fig. S1; (35, 59)]. Note that in the deep location, we removed the debris flow interval between ~543 and 550 mbsf to obtain an undisturbed record. For each sedimentary record, sedimentation rates were inferred at the depth of each δ34Spyr measurement by calculating the linear sedimentation rates implied by the age models.

Sedimentological description.
IODP 1352

The IODP borehole 1352 can be divided into three distinct sedimentological units. Unit I spans the Pleistocene and consists of a 680-m-thick mud-rich sediment. Detailed inspection of sedimentary facies and foraminiferal stratigraphy reveals that dominant clay sedimentation was interrupted by the deposition of several sandy mud layers, rich in planktic foraminifera and other biogenic components, during each of the major interglacial isotope stages. These coarser-grained, more biogenic sediments reflect highstand sedimentation, and contrast with the fine-grained, terrigenous muds deposited during glacial periods. The basal contacts of the coarser-grained units with the fine-grained muds are sharp, and their contacts with overlying muds are characterized by sandy mud with abundant shells and shell fragments, interpreted as a firm ground that marks the initiation of transgressive facies. Evidence for bioturbation and bioirrigation is mostly described in the muddy units, corresponding to periods of low sea level, and often intensified below the reworked or erosional uncemented firm grounds (37).

Although the transition between units I and II is located around 680 mbsf, the transition is gradual, reflecting a gradual deepening to a slope setting. Unit II contains a gradual change from uncemented calcareous sandy mud to lithified marlstone with less frequently intercalated sandstone, except at the base of the unit (Miocene) where glauconitic sandstone becomes more common. An erosional surface, known as the Marshall Paraconformity, occurs at the base of unit II (1750 mbsf), where there is an abrupt change to lithologic unit III, which is composed of nannofossil micritic limestone (hemipelagic to pelagic foraminifera) of Eocene age.

ODP 1123

The offshore ODP borehole 1123 penetrated 633 mbsf, recovering a succession of clay-rich nannofossil ooze, chalk, and limestone. At the top, unit I corresponds to a Pleistocene drift sediment sequence, showing glacial-interglacial couplets, composed of 188 m of clayey nannofossil oozes with numerous tephra layers.

Unit II displays lithological features that are characteristic of sediment drift deposits, i.e., composite sequences of a few centimeters to decimeters in thickness showing overall negative grading from muddy through silty-sandy lenticular beds and then positive grading back to muddy drift facies. The thickness of the silty-sandy unit generally increases with depth. Silt mottles, lenses, irregular to subhorizontal layers, and cross-laminated alignments are also preserved. Unit II is sufficiently indurated to be classified as chalk, and it contains fewer tephra beds.

Unit III consists mainly of similar sediments as the above units, except that terrigenous clay is sufficiently abundant to classify the unit as a nannofossil mudstone. The identification of sedimentary features in those units is made difficult by abundant fractures, often along bedding planes or biscuiting in the drilled cores. Nevertheless, we do observe abundant sandy lenticular beds alternating with mud-rich layers, silt mottles, lenses, irregular to subhorizontal layers, and cross-laminated alignments in Miocene sediment, supporting the inference of drift sedimentation (35). The identification of sediment drifts typically relies on the interpretation of reflection seismic datasets. Spectacular, regional-scale, dune crossbedding, with sets up to 5 m thick and tens of meters long, has been identified in the seismic profile through the North Chatham Drift. These features have been interpreted to result from a semipermanent flowing current (27, 60, 61). In addition to the presence of these thick drift deposits, which commonly contain mud waves, furrows, and other active bedforms, there is abundant evidence in the seismic profiles of erosion or nondeposition in the form of scoured bedrock and marginal mots.

A singularity of unit III is the presence of a 7-m-thick chaotic assemblage of plastically deformed clasts of nannofossil chalk at ~543 mbsf, resulting from a debris flow. Below the Marshall Paraconformity, the remainder of the core is composed of unit IV, an Eocene micritic limestone that seems to be correlative with the bottom of the shelf IODP 1352.


Sulfur extraction. Sulfur was extracted using the chromium-reduction method (62). This method allows recovery of all reduced sulfur present in sedimentary samples (denoted CRS and including pyrite, elemental sulfur, and iron monosulfide phases). During extraction, approximately 2 g of bulk sample was reacted with ~25 ml of 1 M reduced chromium chloride (CrCl2) solution and 25 ml of 6 N HCl for 4 hours just below the boiling point, in a specialized extraction line under a N2 atmosphere. Note that for samples with very low CRS content located at the bottom of ODP 1123 borehole, up to 4 g of bulk rock was used during CRS extraction. The liberated hydrogen sulfide was reacted in a silver nitrate (0.1 M) trap, recovering the sulfide as Ag2S with a reproducibility better than 5% for repeated analyses. Precipitated Ag2S was rinsed three times using Milli-Q water, centrifuged, completely dried, and homogenized before analysis.

At least five bulk samples at each location and at a depth range covering the entire core were treated with 25 ml of 6 N HCl for 2 hours to release acid volatile sulfur (AVS; mainly iron monosulfide). The H2S liberated was then converted to Ag2S via reaction with silver nitrate (0.1 M). None of the acid-treated samples yielded enough Ag2S to be weighed. Following these exploratory extractions, we assume that the AVS component is negligible in our samples.

Iron speciation analyses. Iron speciation measurements were made at the Weizmann Institute of Science using the calibrated extraction of (63), as modified for application to recent sediments [e.g., (64, 65)]. The extracted fractions are defined operationally, after (66, 67). The first step consists of a 0.5 N HCl extraction for 1 hour, targeting ferrous iron (FeII) phases, including AVS, surface-bound FeII, and FeII carbonate/phosphate phases. As the sulfur extractions indicated that the AVS component is negligible in our samples, we assume that no FeII is associated with iron monosulfides and assign the FeIIHCl pool to nonsulfidized FeII minerals (i.e., mostly carbonates). The 0.5 N HCl also extracts the most reactive FeIII minerals in the sediment (termed FeIIIHCl), in particular, ferrihydrite. Because of its rapid kinetics of sulfidization (68) and transformation to more crystalline FeIII oxides, it is highly unlikely that the ferrihydrite extracted by 0.5 N HCl was a part of the original sediment. Instead, we consider the FeIIIHCl fraction in this core, which was not stored under anoxic conditions, to most likely represent the oxidation product of pyrite. Accordingly, the FeIIIHCl fraction was added to the chromium-reducible fraction to correct the pyrite abundance for postsampling oxidation. A second sequential extraction step uses a buffered Na-dithionite solution to extract more crystalline FeIII (oxyhydr)oxide minerals, such as goethite and hematite (termed Fedi-ct). A third extraction step targets magnetite, using Na-oxalate (termed Feoxa). A fourth extraction step targets pyrite by reduction with chromium (termed FeCRS). The extraction is similar to that described above for CRS (62). As mentioned above, we added FeIIIHCl to FeCRS to estimate the original amount of pyrite present in the cores (termed Corr.Fepyr).

The FeIIHCl, Fedi-ct and Feoxa and FeCRS fractions were determined by spectrophotometry using a ferrozine assay (69), immediately after completion of the extraction. Determination of the FeIIIHCl abundance requires reduction of all FeIII to FeII, which is achieved by reaction with ascorbic acid. The total FeHCl is then measured by spectrophotometry, and the FeIIIHCl abundance is determined by subtraction of FeIIHCl from the total FeHCl (69).

Sulfur isotopic analysis. For analysis, 450 μg of Ag2S was loaded into tin capsules with excess V2O5. The 34S/32S ratio was measured on a Thermo Fisher Scientific Delta V Plus, following inline combustion in a Costech ECS 4010 elemental analyzer at Washington University in St. Louis. CRS sulfur isotope compositions are expressed in standard delta notation as per mil (‰) deviations from the VCDT standard, with an analytical error of <0.5‰.

Organic carbon concentration and isotopic analysis. Before organic carbon extraction, the carbonate mineral fraction was removed from bulk samples using excess 1.5 N HCl digestion for 48 hours following (70). During digestion, centrifuge tubes were placed in an ultrasonic bath to increase the mechanical separation of clay and calcium carbonates. Total dissolution residues were washed three times with distilled water, centrifuged, and then completely dried at 50°C. The residual powders were homogenized, and before analysis, 30 mg was loaded into tin capsules. Analyses were performed using an elemental analyzer (Thermo Fisher Scientific Flash EA 2000) coupled to an isotope ratio mass spectrometer (Thermo Fisher Scientific Delta V+ EA-IRMS) at Washington University in St. Louis. TOC concentrations were measured using the thermal conductivity detector of Flash EA 2000. Carbon isotope ratios are reported in standard delta notation as per mil deviations from the Pee Dee Belemnite standard, with an analytical error of <0.2‰ (1σ) for organic carbon isotopes.

Comparison with existing data. The sulfur isotopic composition of a few total reducible inorganic sulfur (TRIS; n = 5) samples at the ODP 1123 site has been measured on board the research vessel and reported (71). As these data have not been released in tabulated form, to allow comparison with our measurements, we digitized figure 5 of (71) using DigitizeIt ( To account for the uncertainties associated with such an exercise, we performed three independent digitizations and report the mean and SD of the three replications (table S3).

As shown in fig. S3, there is overall good agreement between our data (in blue) and the sparser TRIS data (in purple), although we do observe a few per mil differences between the two datasets at depths where they overlap. Differences in the reference frame between the datasets and inaccuracies because of the data digitization are possible but minor compared to the differences in the reported δ34Spyr values. We suggest that the observed mismatch may reflect the freshness of the on-board samples relative to our “aged” samples. The fresh samples probably contained porewater H2S and/or AVS (mostly as iron monosulfides), both of which are sensitive to oxidation and would not survive long-term exposure to oxygen during storage. We performed a few exploratory acidifications of samples from both cores, which yielded immeasurable amounts of AVS. Consequently, in the rest of the samples, we did not attempt to extract AVS, and all reduced sulfur was extracted as CRS. We suggest that our CRS extractions were composed mostly of pyrite, whereas the on-board TRIS results also contained H2S and/or AVS. Isotopic differences between these sulfur species are the likely cause of differences in reported δ34Spyr values. Despite the effect of years of exposure to oxygen on the pyrite abundance, the negligible sulfur isotopic fractionation associated with aerobic oxidation of pyrite (72) allows interpretation of the sulfur isotopic variability preserved in the remaining pyrite as negligibly altered [e.g., (73)].

One main point of disagreement seems to be preserved in the data point immediately overlying the Marshall Paraconformity, which is the only one that has a δ34Spyr value that is more negative than our nearest sample (by ~10‰). Without access to complementary analyses from the on-board study (e.g., sediment description–TOC–δ13Corg–Fe speciation), we cannot confidently discuss the possible causes of this mismatch. Adding to the uncertainty, this single sample appears to be deposited at the onset of drift deposition (i.e., immediately after the Marshall Paraconformity). This is a time when one may expect more energetic conditions, more frequent and more intense reworking, and also greater uncertainty in the age model.

Δpyr calculation. The offset between sulfur isotopes in sulfate and pyrite preserved in the sediments (Δpyr) was calculated as followsΔpyr=δ34Sswδ34Spyrwhere δ34Ssw is the coeval sulfur isotopic composition of seawater sulfate, which was obtained by linear interpolation of binned and locally averaged δ34SCAS values from (11). The δ34Spyr values and the interpolated δ34SCAS values are shown in fig. S3.

Sequence stratigraphic framework. We use the stratigraphic concepts and terms in which a seismic depositional sequence is defined by related sediment packages bounded above and below by unconformities or sequences boundaries, presumably related to variations in sea level (74). These genetically related packages are further classified by the relative sea level during their deposition into lowstand, transgressive, and highstand systems tracts (fig. S9).

Following the work of (37), stratigraphic trends in δ18Oflem, Ca/Sr, and Ca/Ti elemental ratios permit correlation of the lithofacies succession to global sea-level variability. Briefly, δ18Oflem was used to estimate past variation in sea level. We used the Ca/Ti ratio to estimate down-core variation in the CaCO3 content. Carbonate-poor intervals are commonly related to low Ca/Ti ratio, whereas high ratios are interpreted to represent carbonate-rich sections (75). The Ca/Sr ratio was used to distinguish between biogenic and detrital carbonate. As Sr is fixed by calcifying organisms at the same time as Ca, decoupling between Ca/Ti and Ca/Sr is indicative of a high detrital carbonate input (75).

In the boreholes, sequence boundaries occurred systematically at high δ18Oflem, when the Ca/Ti and Ca/Sr ratios return to background values. The maximum flooding surfaces were defined by a coeval increase in δ18Oflem, Ca/Ti, and Ca/Sr, whereas transgressive surfaces are marked by minimum δ18Oflem, and low Ca/Ti and Ca/Sr ratios (fig. S9).


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We acknowledge J. Houghton for help with sample preparation and analysis. Funding: V.P. was supported by a Postdoctoral Fellowship from the Dean of the Chemistry Faculty at the Weizmann Institute of Science. Author contributions: V.P., R.N.B., D.A.F., and I.H. conceived the work. V.P. and R.N.B. carried out geochemical analyses. V.P., R.N.B., D.A.F., and I.H. contributed to interpretation of the results and wrote the paper. Competing interests: The authors declare that that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. All geochemical and isotopic data on IODP 1352 and ODP 1123 samples included in this study are available in the main text and in the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article