The Cambrian strata in the Arksu region, western Tarim Block, consist of, in ascending order, the Yurtus (YF), Xiaoerblak (XF), Wusongger (WF), Shayilike (SF), Awatage (AF), and Qiulitage (QF) formations. The Yurtus Formation unconformably overlies the Ediacaran Qigeblak Formation and conformably underlies the Xiaoerblak Formation (Figure 1) (He et al., 2010, 2018; Zhu et al., 2018). The basal Yurtus Formation is restricted to the Tommotian (equivalent to Stage 2 of Terreneuvian Series in the Geological Time Scale) by the appearances of Asteridium-Heliosphaeridium-Comasphaeridium acritarch assemblage (Yao et al., 2005; Dong et al., 2009). Univalve mollusk fossils (Shabaktiella multiformis, Parcaconus xinjiangensis, Eoyochelcionella aksuensis, etc.) discovered from the upper Yurtus Formation might be correlated with the small shelly fossils from the Qiongzhusian strata in South China or the Atdabanian Stage in the Eastern European Platform (equivalent to Stage 3 of Cambrian Series 2 in the Geological Time Scale) (Qian, 1999; Qian et al., 2000; Gradstein et al., 2012). Furthermore, the absence of trilobite suggests that the Yurtus Formation might belong to Cambrian Stage 2 and probably predate Cambrian Stage 3, i.e., between 529 and 521 Ma (Dong et al., 2009). Samples were collected from the drill core (X1). The Yurtus Formation in the X1 core is composed of, in stratigraphic order, bedded chert (∼3.7 m), black shale (∼23 m), intercalated mudstone and dolostone (∼13 m), and dolostone (∼13 m) lithological units. In this study, only the black shale samples from the lower Yurtus Formation were analyzed.

Pyrite sulfur isotope ratios were determined at the State Key Laboratory of Biogeology and Environmental Geology, China University of Geosciences (Wuhan). The purified Ag₂S precipitate (after chromium reduction) was mixed with an excessive amount of V₂O₅ 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.

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 δ³⁴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)⁻¹ can simulate the porewater sulfate concentration profile.

In addition, the porewater profile of sulfur isotope composition of porewater sulfate (δ³⁴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 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 δ³⁴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₂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₂S oxidation (Lang et al., 2020), recycling of H₂S-derived sulfate in DSR is excluded. In addition, H₂S oxidation could be further complicated by its possible connections with Fe/Mn oxidation, nitrate oxidation, and O₂ oxidation (Li et al., 2010). 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₂S oxidation, i.e., more reduced conditions. Furthermore, the oxidation of H₂S is simulated by the H₂S upward diffusion to seawater instead of the downward diffusion of oxidants into sediments. For these reasons, the profiles of H₂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]⁻¹ 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.

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 of parameters. If there are n parameters, each of which has m_i possible values (i ranges from 1 to n). There are ∏ i n m i assemblages or outputs.

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 δ³⁴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 δ³⁴S_py and Fe speciation data. Below, we will calculate the early Cambrian seawater sulfate concentration by this methodology.

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 non-sulfidic 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 non-sulfidic 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 δ³⁴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 (Wang et al., 2012). 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 δ³⁴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 δ³⁴S_SW is set to +30‰ (Fike et al., 2015). Although the validity of δ³⁴S_CAS has been challenged (Marenco et al., 2008; Peng et al., 2014), the assigned δ³⁴S_SW value is consistent with the δ³⁴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.