Quantifying the Seawater Sulfate Concentration in the Cambrian Ocean

Although the earliest animals might have evolved in certain “sweet spots” in the last 10 million years of Ediacaran (550–541 Ma), the Cambrian explosion requires sufficiently high levels of oxygen (O2) in the atmosphere and diverse habitable niches in the substantively oxygenated seafloor. However, previous studies indicate that the marine redox landscape was temporally oscillatory and spatially heterogeneous, suggesting the decoupling of atmospheric oxygenation and oceanic oxidation. The seawater sulfate concentration is controlled by both the atmospheric O2 level and the marine redox condition, with sulfide oxidation in continents as the major source, and sulfate reduction and pyrite burial as the major sink of seawater sulfate. It is thus important to quantify the sulfate concentration on the eve of the Cambrian explosion. In this study, we measured the pyrite contents and pyrite sulfur isotopes of black shale samples from the Yurtus Formation (Cambrian Series 2) in the Tarim Block, northwestern China. A numerical model is developed to calculate the seawater sulfate concentration using the pyrite content and pyrite sulfur isotope data. We first calibrate some key parameters based on observations from modern marine sediments. Then, the Monte Carlo simulation is applied to reduce the uncertainty raised by loosely confined parameters. Based on the geochemical data from both Tarim and Yangtze blocks, the modeling results indicate the seawater sulfate concentration of 8.9–14 mM, suggesting the seawater sulfate concentration was already 30–50% of the present level (28 mM). High seawater sulfate concentration might be attributed to the enhanced terrestrial sulfate input and widespread ocean oxygenation on the eve of the Cambrian explosion.


INTRODUCTION
The seawater sulfate concentration is a critical indicator of the redox condition in the atmosphereocean system. On the one hand, the seawater sulfate concentration is controlled by both marine redox condition and atmospheric O 2 level, because oxidative weathering of sulfide in continents is one of the major sources, and sulfate reduction and pyrite burial represent one of the major sinks of seawater sulfate (Canfield, 2004;Canfield and Farquhar, 2009). On the other hand, the seawater sulfate concentration should be globally homogeneous, reflecting the overall global redox condition of the atmosphere-ocean system. Thus, reconstruction of seawater sulfate concentration would provide a direct constraint on the global ocean redox condition and the atmospheric O 2 level.
The second rise of atmospheric O 2 level occurred in the late Neoproterozoic, coined the Neoproterozoic oxygenation event (NOE) (Shields-Zhou and Och, 2011). NOE is supported by several lines of geochemical evidence. The enrichment of redoxsensitive elements (e.g., V, U, Mo) in the early Ediacaran black shales implies the oxidation of the deep ocean immediately after the Marinoan Snowball Earth glaciation (Sahoo et al., 2012;Sahoo et al., 2016). The global occurrence of Shuram Excursion, the largest negative carbon isotope excursion in Earth's history, has been interpreted as massive oxidation of dissolved organic carbon (DOC) in the Ediacaran deep ocean (Fike et al., 2006;Kaufman et al., 2007;McFadden et al., 2008;Grotzinger et al., 2011). Furthermore, the decrease in the reactive Fe content in deepsea deposits also indicates oxidation of the deep ocean after the Ediacaran Gaskiers glaciation (580 Ma) (Canfield et al., 2007). NOE was associated with the dramatic change in the biosphere. For example, biomarker data indicate the increase in eukaryotic primary productivity in the nonglacial interlude between the two Cryogenian (720-635 Ma) Snowball Earth glaciations, while paleontological data indicate the diversification of eukaryotes in the earliest Ediacaran and the subsequent evolution of multicellular organisms, e.g., macroscopic algae and Ediacara biota (Glaessner, 1984;Zhang, 1989;Narbonne, 2005;Yin et al., 2007;Yuan et al., 2011;Liu et al., 2014;Brocks et al., 2017).
Although the biological evolution and the atmosphereocean oxygenation were broadly coincident in the geochemical and paleontological records, more and more studies indicate inconsistent or even contradictory results drawn from different proxies. It is proposed that the redox landscape in the Ediacaran and early Cambrian ocean might be temporally dynamic and spatially heterogeneous Sahoo et al., 2016;Li et al., 2018), and the ocean was predominantly anoxic and was frequently punctuated by episodic or sporadic oxidation or euxinia (Sahoo et al., 2016;Zhang et al., 2018;Ding et al., 2019).
The seawater sulfate concentration would provide the key evidence to resolve the inconsistency between different proxies. However, the seawater sulfate concentration cannot be directly measured from sedimentary rocks. Based on the stratigraphic variation of sulfur isotopes of carbonate-associated sulfate (CAS, δ 34 S CAS ) (Kah et al., 2004), low seawater sulfate concentration of ∼2 mM in the early Cambrian ocean was proposed (Loyd et al., 2012;Thompson and Kah, 2012). In contrast, the marine sulfur mass balance model indicates a higher seawater sulfate concentration of ∼10 mM. In the latter scenario, it is suggested that the high seawater sulfate concentration might be attributed to the invention of bioturbation during the Cambrian explosion (Canfield and Farquhar, 2009). Such contradictory results prevent further discussion of marine-atmosphere redox coupling on the eve of the Cambrian explosion.
In this study, we develop a new method to quantify the seawater sulfate concentration by using Fe speciation and pyrite sulfur isotope data. We analyzed the black shale of the lower Cambrian Yurtus Formation in the western Tarim Block, northwestern China. Combining with geochemical data from the Yangtze Block, the seawater sulfate concentration in the early Cambrian ocean was quantified.

Pyrite Sulfur Isotope Analysis
Pyrite sulfur isotope ratios were determined at the State Key Laboratory of Biogeology and Environmental Geology, China University of Geosciences (Wuhan). The purified Ag 2 S precipitate (after chromium reduction) was mixed with an excessive amount of V 2 O 5 and was wrapped in a tin cup. S isotope ratios were determined by a Thermo Instruments Delta V Plus isotope ratio mass spectrometer coupled with a Costech elemental analyzer. S isotope values are reported by delta notation as per mil (‰) deviation relative to the V-CDT (Vienna-Cañon Diablo Tribolite) international standard. Samples were calibrated by international standards: IAEA S1 (−0.3‰), IAEA S2 (22.65‰), and IAEA S3 (−32.5‰). The analytical precision is ∼0.1‰ (1σ), which was determined by repeated analyses of IAEA international standards.

RESULTS
The Fe speciation and pyrite sulfur isotope data are tabulated in Table 1, and the Fe speciation data have been reported in Zhu Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 767857 2 et al. (2021). The stratigraphic profiles of geochemical data are illustrated in Figure 2. Pyrite sulfur isotope (δ 34 S py ) ranges from -8.2‰ to +16.1‰ (mean +4.2‰, n 17), with most being higher than 5.0‰ in the lower black shale interval. The δ 34 S CAS values vary between +26.2‰ and +36.1‰ (mean +30.4‰, n 12). Pyrite contents range from 0.05 to 8.62% (mean 2.26%, n 17). These samples have wide ranges of variation of Fe T contents, ranging from 0.49 to 4.40% (mean 1.51%, n 17). Fe HR contents range between 0.25 and 4.48% (mean 1.27, n 17), with the maximum value at the uppermost black shale. Fe py contents range between 0.02 and 4.02% (mean 1.06, n 17). The Fe HR /Fe T ratios vary between 0.20 and 1.00 (mean 0.81, n 17). Fe py /Fe HR ratios vary between 0.07 and 0.90 (mean 0.75, n 17).

DISCUSSION
A numerical model was developed to simulate the syndepositional pyrite formation in sediment porewater (Lang et al., 2020). The syndepositional pyrite formation involves dissimilatory sulfate reduction (DSR) in sediment porewater with sulfate supply from seawater diffusion. H 2 S, the product of DSR, either converts to pyrite by reacting with reactive Fe in sediment or is oxidized back to sulfate. The modeling results indicate that both δ 34 S py and pyrite contents in sediments or sedimentary rocks are sensitive to some environmental factors, including the seawater sulfate concentration, the isotopic composition of seawater sulfate (δ 34 S SW ), reaction rate constant of DSR (R DSR ), biological S isotope fractionation in DSR (Δ DSR ), organic matter contents in sediments (Org0), reactive Fe contents in sediments (Fe0), reaction rate constant of pyrite formation (R py ), and seafloor redox conditions. The latter controls the fraction of H 2 S being oxidized by sulfur species with higher valence states (Lang et al., 2020). It should be noted that, for the simplicity, H 2 S oxidation exclusively generates sulfate that returns to the seawater sulfate inventory, while intermediate sulfur species, such as elemental S or thiosulfate (Canfield and Thamdrup, 1994;Habicht et al., 1998), are not considered.
The purpose of Lang et al. (2020) is to display how δ 34 S py and pyrite contents are affected by the process of syndepositional pyrite formation in sediment porewater. The parameters used in the model do not represent the reactions that occurred in the natural environment. Thus, the original model only provides a theoretical framework but cannot be applied directly to quantify the marine sulfur cycle in deep time. In this study, we refine Lang's model to develop an applicable method to quantify the seawater sulfate concentration by using δ 34 S py and Fe speciation data. The differential equations that describe porewater sulfate, organic matter, reactive Fe, and H 2 S profiles have been described in detail in Lang et al. (2020) and thus will not be illustrated here (for the detailed mathematics, please see Supplementary Material).
In this study, we made the following modifications. 1) We calibrated the key parameters, including R DSR and Δ DSR , based on the porewater geochemical profiles of modern marine sediments. 2) We applied the Monte Carlo simulation to eliminate uncertainties raised by other loosely constrained parameters, including sedimentation rate, initial organic matter contents in sediments, R py , and redox conditions at the seafloor. The inputs from the sample measurements include δ 34 S py , Fe speciation data, and δ 34 S SW (from CAS data).

Parameter Calibration
With the inputs of sedimentation rate, concentration, and isotopic composition of seawater sulfate, the initial organic carbon and reactive Fe contents in sediments, and the redox condition in seawater/seafloor, both δ 34 S py and pyrite contents in sediments/sedimentary rocks can be calculated (Lang et al., 2020). There are two key parameters that would significantly affect the modeling outputs: the reaction rate constant of DSR (R DSR ) and isotope fractionation in DSR (Δ DSR ). To apply this model to quantify the deep time seawater sulfate concentration, we need to calibrate these two parameters by porewater geochemical profiles of modern sediments in Santa Barbara Basin (Soutar and Crill, 1977;Raven et al., 2016). The porewater profile of sulfate concentration is controlled by the seawater sulfate concentration (Org0) and R DSR . The former two values are already known for modern sediments. Thus, R DSR can be determined by fitting the porewater sulfate concentration profile. The modeling results are illustrated in Figure 3, showing that D DSR 0.0033 ((mM/L)*ka) −1 can simulate the porewater sulfate concentration profile.
In addition, the porewater profile of sulfur isotope composition of porewater sulfate (δ 34 S pw ) is controlled by Δ DSR , R DSR , (Org0), and seawater sulfate concentration. It is proposed that Δ DSR is 46‰ if sulfate concentration is greater   Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 767857 5 than 0.2 mM (Habicht and Canfield, 1997;Habicht et al., 2002). Recent studies indicate that the isotopic difference between seawater sulfate and sediment pyrite can be up to 70‰ in the past 500 ky (Pasquier et al., 2017). Thus, we allow Δ DSR varying between 46 and 70‰ ( Figure 4). The modeling results indicate that Δ DSR of 70‰ can best simulate the porewater δ 34 S pw profile of modern sediments. In fact, Δ DSR of 70% is consistent with the cultural DSR experiments, which indicate a similar magnitude (up to 66‰) of Δ DSR , if there is no sulfur disproportionation (Sim et al., 2011).
Another key issue is the correlation between seawater/seafloor redox condition and pyrite formation reaction constant (R py ). R py is a constant, but it is a variant in the model that is related to the seawater redox condition. Oxidation of H 2 S is a complex process, involving the generation of different sulfur species of various valence states. Since we have ignored the formation of sulfur species with intermediate valence states and assumed sulfate as the only product of H 2 S oxidation (Lang et al., 2020), recycling of H 2 S-derived sulfate in DSR is excluded. In addition, H 2 S oxidation could be further complicated by its possible connections with Fe/Mn oxidation, nitrate oxidation, and O 2 oxidation . In order to simplify the simulation, we replace the seawater/seafloor redox condition with R py . It is proposed that a larger R py implies less H 2 S oxidation, i.e., more reduced conditions. Furthermore, the oxidation of H 2 S is simulated by the H 2 S upward diffusion to seawater instead of the downward diffusion of oxidants into sediments. For these reasons, the profiles of H 2 S concentration and isotope profiles can't be well simulated in our model ( Figure 5). In the simulation, we have R py 0.3 [(mM/L)*ka] −1 as a cutoff value for euxinic bottom water, which is closest to the actual profiles with the seafloor oxygen fugacity (f O2 ) < 10 uM ( Figure 5). It should be noted that the linkage between R py and seafloor f O2 is not fixed due to the lack of modern observation.

The Monte Carlo Simulation
Because the number of parameters is much larger than the number of equations, there are multiple solutions in quantifying the seawater sulfate concentration with pyrite sulfur isotope and pyrite contents. To reduce the uncertainties raised by loosely constrained parameters, here, we apply the Monte Carlo simulation. In this method, each parameter is allowed a range of variation ( Table 2). All possible solutions are calculated based on each assemblage To save the computation resource, we take the following assumptions. 1) Sedimentation rate of similar lithology in the same section is limited to a narrow range (0.05-0.15 m/ky for our samples). This assumption is generally consistent with the estimation of sedimentation rate based on the available biostratigraphic data. In detail, the Yurtus Formation with a total thickness of ∼60 m was deposited within 21 million years, including black shale of ∼30 m, dolomite of 3.3 m, and carbonate dominated interval of ∼27 m (Zhu et al., 2021). Considering the faster precipitation of carbonate rocks, it is estimated that the black shale was deposited between 10.5 and 21 million years. Assuming the compaction rate of 0.8, the sedimentation rate of muddy sediment (black shale) in the Yurtus Formation is between 0.075 and 0.15 m/ky. 2) Pyrite is the only major sink for reactive Fe and precipitates in sediment porewater. The second assumption allows the replacement of pyrite content by Fe PY /Fe HR ratio. In addition, this also implies (Fe0) equals to Fe HR . For each sample, if Fe PY /Fe HR and δ 34 S py fall between the two nearby input points, the input points are included in the solution.
In addition, to limit the multiplicity of solution from multiple samples, we have the following assumptions. 1) The seawater sulfate concentration is invariant for samples from the same section (non-sulfidic seawater). This is likely the case given the residence time of seawater sulfate is longer than the duration of sample collections (the Geological Background section). 2) The seawater/seafloor redox condition was the same for all samples of the same lithology. For each set of seawater sulfate concentration and seawater/seafloor redox condition values, we calculate the frequency that the solution set of samples contains the seawater sulfate concentration range and the seawater redox condition range (Table 3). Because sulfate concentration is homogeneous in non-sulfidic seawater, the seawater sulfate concentration range with the frequency of 1 (indicating possible for all samples) is the plausible range of seawater sulfate concentration.
Technically, we justified the two parameters, R DSR and Δ DSR , based on the simulations of porewater geochemical profiles of modern marine sediments. Furthermore, to eliminate uncertainties raised by other loosely constrained or unconstrained parameters, such as sedimentation rate and redox condition, we run the Monte Carlo simulation to calculate all possible outcomes. With certain constraints from geological observations and by assuming homogeneous seawater sulfate concentration, the seawater sulfate concentration of the early Cambrian ocean can be quantified by δ 34 S py and Fe speciation data. Below, we will calculate the early Cambrian seawater sulfate concentration by this methodology.

Quantifying the Early Cambrian Seawater Sulfate Concentration
The sulfate concentration should be homogeneous in non-sulfidic seawater because the ocean mixing time is four orders of magnitude shorter than the residence time of seawater sulfate (1000s of years vs. 10s million years). Even the seawater sulfate concentration was an order of magnitude lower (Loyd et al., 2012), a homogeneous seawater sulfate concentration in the nonsulfidic seawater is still valid in a million-year time scale. In contrast, because of active DSR in sulfidic seawater, the sulfate concentration in sulfidic seawater could be different from nonsulfidic seawater (Lyons, 1997). Thus, this model can only be applied to the non-sulfidic ocean.
It is noticed that some Yurtus samples have Fe PY /Fe HR ratios >0.8, suggesting the deposition under sulfidic conditions (Poulton and Raiswell, 2002;Raiswell et al., 2018). However, deposition under sulfidic conditions is inconsistent with abundant fossils from the Yurtus Formation, including typical early Cambrian Asteridium-Heliosphaeridium-Comasphaeridium acritarch assemblage, and tubular fossil Megathrix, as well as small shelly fossils including Shabaktiella multiformis, Parcaconus xinjiangensis, Eoyochelcionella aksuensis (Qian, 1999;Qian et al., 2001;Yao et al., 2005;Dong et al., 2009), the latter of which implies the seafloor might be habitable for animals. In fact, recent studies indicate that the terrestrial input of reactive Fe is strongly affected by the intensity of continental weathering (Wei et al., 2021), implying that the Fe HR /Fe T ratio is not only affected by the redox condition of seawater. Similarly, pyrite precipitation within sediment porewater is affected by many factors, and a high fraction (high Fe PY /Fe HR ratio) of syndepositional pyrite could be generated in non-sulfidic conditions (Lang et al., 2020). Deposition in a non-sulfidic water column is also consistent with abundant diagenetic euhedral pyrites from the Yurtus Formation. Therefore, we argue that the above simulations can be applied to the Yurtus samples.  To validate the calculated seawater sulfate concentration from the Tarim samples, we also choose samples from the Yangtze Block, South China. The early Cambrian successions in the Yangtze Block have been extensively studied, and several sections have both δ 34 S py and Fe speciation data reported, including the Xiaotan section (Och et al., 2016), Shatan and Songtao sections (Goldberg et al., 2007), Jinsha and Weng'an sections (Jin et al., 2016), Dingtai section (Xu et al., 2012), Yangjiaping section (Feng et al., 2014), and Longbizui section . As discussed above, this model can only simulate syndepositional pyrite formation in non-sulfidic seawater. Although there are many samples from the Yangtze Block that might have been deposited under non-sulfidic conditions, the sampling interval from the Yangtze Block should overlap with that of the Yurtus Formation in the Tarim Block. Thus, we choose data from the Songtao section for analysis (Goldberg et al., 2007).
In the simulation, we use Fe speciation data and δ 34 S py data from the Yurtus Formation in the X1 drill core in the Tarim Block and from the Niutitang Formation (∼529−515 Ma) (Jin et al., 2020;Na and Kiessling, 2015) in the Songtao section in the Yangtze Block. We also apply the calibrated R DSR and Δ DSR values (Figures 3,4). The early Cambrian δ 34 S SW is set to +30‰ (Fike et al., 2015). Although the validity of δ 34 S CAS has been challenged (Marenco et al., 2008;Peng et al., 2014), the assigned δ 34 S SW value is consistent with the δ 34 S CAS data from the Yurtus Formation ( Figure 2). Other parameters are allowed a range of variation, which are listed in Table 1.
Based on the outputs of the simulation, the frequency map of both seawater sulfate concentration and redox condition for all samples can be created ( Table 3). The seawater sulfate concentration is bracketed between 8.9 and 14.0 mM. These values are in agreement with a rough estimate of ∼10 mM based on the marine sulfur isotope mass balance calculation (Canfield and Farquhar, 2009), but at least 4 times higher than the calculation of ∼2 mM based on the stratigraphic variation of CAS sulfur isotope values (Loyd et al., 2012).
Our modeling result indicates that the seawater sulfate concentration was already high during the early Cambrian, equivalent to ∼30-50% of the present level of 28 mM. The seawater sulfate concentration is controlled by both terrestrial input of sulfate and burial in the ocean. Evaporate deposition (mainly gypsum) and pyrite precipitation and burial are the two major sinks of seawater sulfate (Fike et al., 2015). It is unclear how seawater sulfate concentration could increase to 30-50% of the present level. Neither is known about when the seawater sulfate concentration increased during the EdiacaranCambrian transition. We suggest that high seawater sulfate concentration during the early Cambrian might be favored due to the following reasons. First, the terrestrial sulfate input might have dramatically increased due to extensive evaporate dissolution during the formation of "Great Unconformity" (Peters and Gaines, 2012;Shields et al., 2019). Second, although the redox landscape in the early Cambrian ocean might be highly heterogeneous with the possible development of sulfidic wedge in continental margins (Jin et al., 2014;Jin et al., 2016;Li et al., 2018), the proportion of sulfidic seafloor might be smaller than that of early Neoproterozoic and Mesoproterozoic (Lyons et al., 2014). Limited seafloor euxinia also indicates less efficient pyrite burial, reducing the pyrite sink of the seawater sulfate. This argument is consistent with high Mo content in the lower Cambrian black shales (Sahoo et al., 2012;Reinhard et al., 2013;Sahoo et al., 2016). Finally, our model provides a new approach to quantify the seawater sulfate concentration in paleoceans. In addition to justifying some key parameters based on the modern sediment observation, the Monte Carlo simulation could reduce the uncertainties raised by loosed constrained parameters. It should be noted that the assumption of the Monte Carlo simulation is invariant seawater sulfate concentration during the interval of simulation. Thus, high-resolution sampling from non-sulfidic deposits is required. If samples were collected from multiple sections, the chrono-and/or biostratigraphic framework is required to justify the coeval deposition.

CONCLUSION
In this study, we justified the key parameters (reaction rate constant of DSR and sulfur isotopic fractionation in DSR) of the syndepositional pyrite formation model. We also develop the Monte Carlo simulation approach to avoid uncertainties raised by loosely constrained parameters, such as sedimentation rate and the initial organic matter content in sediment. The new model allows the quantification of seawater sulfate concentration in deep time by using pyrite sulfur isotope and Fe speciation data. Based on the study of the lower Cambrian Yurtus Formation in the Tarim Block, combining with the data of the coeval Niutitang Formation in the Yangtze Block, the early Cambrian seawater sulfate concentration is bracketed between 8.9 and 14.0 mM, approaching to 30-50% of the present level. The relatively high seawater sulfate concentration might be attributed to enhanced terrestrial sulfate input in the context of "the Great Unconformity" and the reduced sulfidic seafloor in the second rise of atmospheric O 2 level. Our model provides a new approach to quantify the seawater sulfate concentration in paleoceans.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
BS, TH and XL contributed to conception and design of the study. GZ, TL and KZ collected and analysed samples. TH, BS, XL, WT, and RW contributed to modeling and data analyzing. TH and BS wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.