Sedimentary Pyrite Formation in a Seasonally Oxygen‐Stressed Estuary: Potential Imprints of Microbial Ecology and Position‐Specific Isotope Fractionation

Sedimentary pyrite records are essential for reconstructing paleoenvironmental conditions, but these records may be affected by seasonal fluctuations in oxygen concentration and temperature, which can impact bioturbation, sulfide fluxes, and distributions of sulfide oxidizing microbes (SOMs). To investigate how seasonal oxygen stress influences surficial (<2 cm) pyrite formation, we measured time‐series concentrations and sulfur isotope (δ34S) compositions of pyrite sulfur along with those of potential precursor compounds at a bioturbated shoal site and an oxygen‐deficient channel site in Chesapeake Bay. We also measured radioisotope depth profiles to estimate sedimentation rates and bioturbation intensities. Results show that net pyrite precipitation was restricted to summer and early autumn at both sites. Pyrite concentration was higher and apparently more responsive to precursor compound concentration at the mildly bioturbated site than at the non‐bioturbated site. This disparity may be driven by differences in the dominant SOM communities between the two sites. Despite this, the sites' similar pyrite δ34S values imply that changes in SOM communities have limited effects on surficial pyrite δ34S values here. However, we found that pyrite δ34S values are consistently and anomalously lower than coeval precursor compounds at both sites. A steady‐state model demonstrates that equilibrium position‐specific isotope fractionation (PSIF) effects in the S8‐polysulfide pool can create a 4.3–7.3‰ gap between δ34S values of pyrite and zero‐valent sulfur. This study suggests that SOM communities may have distinct effects on pyrite accumulation in seasonally dynamic systems, and that PSIF in the polysulfide pool may leave an imprint in pyrite isotope records.


10.1029/2022JG007324
2 of 20 et al., 2019), and the rate of its formation and/or the formation of its precursor compounds may be enhanced by the availability of microbial biomass as a template or nucleation site for crystallization (Donald & Southam, 1999;Ferris et al., 1987;Picard et al., 2018). The sulfide oxidizing microbe (SOM) community in a sediment column may also affect the rate and chemical pathway of pyrite formation because of the distinct effects that different SOM have on sediment geochemistry. For example, cable bacteria tend to dissolve FeS by lowering pore water pH, while Beggiatoa colonizing the same sediments can induce FeS accumulation under a local pH maximum (Seitaj et al., 2015). These changes in pH affect the relative amounts of H 2 S (a pyrite reactant), HS − (not a pyrite reactant), and polysulfide-thereby altering the relative rates of the two primary pyrite precipitation reactions, that of FeS with H 2 S (D. Rickard and Luther, 1997) and that of FeS with polysulfide (Luther, 1991).
The influence of SOM on pyrite precipitation pathways is relevant to interpretation of the sedimentary pyrite record. Distributions of stable sulfur isotopes in pyrite are used to reconstruct geochemical conditions in ancient environments (Fike et al., 2015), but recent work has emphasized the role of localized physical and geochemical forcings in affecting pyrite sulfur isotope ratios (Pasquier, Bryant, et al., 2021;. Because the two pyrite reaction pathways can have distinct effects on the sulfur isotope composition of the resulting pyrite (Butler et al., 2004), interpretation of the pyrite isotopic record is further obfuscated by uncertainty over which pyrite precipitation reaction predominates in sediments. Relative rates of these reaction pathways in sediments are often difficult to ascertain (J. Liu et al., 2020;D. Rickard & Luther, 2007) in part because reactant concentrations are strongly influenced by oxidative sulfur cycling by microbes (Canfield et al., 1998) and in part because the rate of pyrite formation via polysulfide may depend on the surface areas of elemental sulfur and FeS (D. T. Rickard, 1975), which are difficult to measure and may also be influenced by microbial activity.
Pyrite formation and SOM ecology can also be influenced by bioturbation, the physical mixing and ventilation of sediments by burrowing benthic animals (Kristensen et al., 2012). This occurs by physicochemical alteration, for example, vertical mixing and oxidative dissolution of pyrite (Thamdrup et al., 1994), by altering sedimentary depth profiles of sulfide and electron donors (Aller, 1982;Bonaglia et al., 2019), and by influencing SOM community composition. Bioturbation can increase or decrease the metabolic rates of sulfate reduction and sulfide oxidation, depending on the type of bioturbation that occurs (Bonaglia et al., 2013;Reichardt, 1988;Vasquez-Cardenas et al., 2016). Low rates of bioturbation under low-oxygen conditions promote formation of microbial mats dominated by SOM such as Beggiatoa or Thioploca, which exploit the oxygen-sulfide chemical gradient around the sediment-water interface (Levin, 2003). Cable bacteria can colonize non-burrowed sediment to exploit the same chemical gradient (Malkin et al., 2014), though they have also been documented around burrow linings in bioturbated sediment (Aller et al., 2019).
These findings point to a complex relationship between bioturbation, microbial sulfur cycling, and pyrite formation, especially under low-oxygen conditions. The geological pyrite record may reflect sensitivity to variations in the type and intensity of bioturbation and to the distinct sulfur-cycling microbial communities that may have existed in a given sediment column. This is of particular relevance to intervals of Earth history such as the Early Paleozoic, when global oxygen concentrations were relatively low and oxygen minimum zones may have impinged onto continental shelves (Tostevin & Mills, 2020). In such environments, seasonal deoxygenation may have been a common feature of shallow marine and estuarine settings. Therefore, to assess the effects of contrasting bioturbation intensities and distinct SOM communities on pyrite formation rates and pyrite sulfur isotope composition, we analyzed surficial sediments collected at monthly to seasonal intervals from two sites in Chesapeake Bay: site CB4. 3C (38.55505°N, 76.42794°W) in the central channel of the Bay at a depth of 26 m and site CB4.3W (38.55728°N, 76.49402°W) on the Bay's western shoal at a depth of 9 m (Figure 1).
Chesapeake Bay is a large estuary formed by the post-glacial inundation of the lower Susquehanna River basin (Colman et al., 1990) and it is subject to strong seasonal oscillations in geochemical conditions (Breitburg, 1990;Kemp et al., 2005), including summertime oxygen depletion that affects bioturbation rates (Schaffner et al., 1992). Each site has been monitored for water column biogeochemical parameters since the mid-1980s, and compiled monitoring data are available from the Chesapeake Bay Program Water Quality Database (https://www. chesapeakebay.net/what/downloads/cbp-water-quality-database-1984-present). While both sites experience strong seasonal oscillations in temperature and oxygen concentration, the channel site experiences more severe oxygen depletion and bears higher sulfide concentrations near the sediment surface in the summer than the shoal site (Table 1). We used aliquots of samples that had previously been used to document that the channel site tends to be dominated by cable bacteria over Beggiatoa, while the shoal site tends to be dominated by Beggiatoa over cable 10.1029/2022JG007324 3 of 20 bacteria (Malkin et al., 2022). As such, these two sites serve as useful locations to study the interplay between bioturbation, oxidative sulfur cycling, and pyrite formation.

Sampling Methods
Sediment cores were collected with an 8.6-cm diameter Uwitec gravity corer in March-August and October 2017, and again in February, May, and August 2018. Cores were sealed, stored aboard ship in water baths at ambient bottom water temperature (ranging from 5 to 27°C, depending on season) and sectioned under a nitrogen (N 2 ) atmosphere within 18 hr at the University of Maryland Center for Environmental Science, Horn Point Laboratory  (Cambridge, MD). Sectioned sediment samples were placed in centrifuge tubes, centrifuged to remove pore water, sealed in gas-tight bags, and stored at −20°C. Additional core collection was conducted in August 2021 for the purpose of collecting pore water sulfide data and measuring radiotracer depth profiles; initial processing of these cores was conducted in the same way.

Sulfur Extraction and Geochemical Analysis
We used a sedimentary sulfur extraction method (J. Liu et al., 2020) to extract zero-valent sulfur (ZVS), acid volatile sulfur (AVS, predominantly composed of iron monosulfides; Luther, 2005; cf. D. Rickard & Morse, 2005), and chromium reducible sulfur (CRS, predominantly composed of pyrite; Canfield et al., 1986). Samples were thawed in gas-tight bags in a refrigerator at 4°C for 24 hr and then processed inside a N 2 -purged glove bag. Two aliquots of equal volume were taken from each sample. The first aliquot was dried in an oven at 70°C for 24 hr to measure dry bulk mass, while the second aliquot was mixed into a centrifuge tube with 1.0 mL of a N 2 -purged 3% Zn acetate/0.1 M acetic acid solution. The Zn acetate solution reacted with the sediment sample for 30 min in the glove bag to stabilize the AVS fraction, after which 5.0 mL of a 3:1 volumetric mixture of N 2 -purged methanol:toluene was pipetted into the tube. Tubes were then placed inside gas-tight barrier bags and placed on a shaker table at 150 revolutions per minute (RPM) for 18 hr. After shaking, samples were centrifuged for 25 min at 4,200 RPM; supernatant fluids were decanted into round-bottom flasks in a N 2 -purged glove bag. The sediment fraction in each tube was placed into a separate round-bottom flask which was subjected to sequential AVS-CRS extraction (described below), while the flasks containing solvent were left sealed in the glove bag until the first AVS-CRS extractions were complete.
Sequential AVS-CRS extractions were conducted at 70°C for 2 hr under N 2 after Canfield et al. (1986). AVS was extracted with 40 mL of N 2 -purged 6 N HCl, and CRS was extracted with an additional 40 mL of 1 M CrCl 2 /0.5 M HCl. Sulfide was transported via N 2 carrier gas to a trap vessel containing 3% Zn acetate/0.1 M acetic acid. The resulting zinc sulfide was transformed to silver sulfide (Ag 2 S) via reaction with excess silver nitrate. Solvent-extractable sulfur (i.e., the ZVS fraction) was extracted from methanol:toluene solvent via hot CrCl 2 reduction in the same manner as the CRS fraction, except that this extraction injected 20 mL of N 2 -purged 6 N HCl and 30 mL of 1 M CrCl 2 /0.5 M HCl solution at the same time. Ag 2 S samples produced from the AVS, CRS, and ZVS extractions were purified with 1 M NH 4 OH (Firsching, 1961), triple-rinsed with deionized water, centrifuged, and dried at 70°C. Dried Ag 2 S was weighed on a microbalance to determine sulfur species concentrations. Concentrations of AVS and CRS for the same sample set were separately measured in Malkin et al. (2022), but the aliquots measured in that study were not reacted with organic solvent to separate ZVS from CRS and their sulfur isotope compositions were not analyzed. The standard error (σ/√2) in solid-phase sulfur concentration measurement in this study's data, based on seven replicates, was 0.14% for pyrite, 0.02% for ZVS, and <0.01% for AVS.
For δ 34 S analysis, 200-400 μg aliquots of dried Ag 2 S were weighed into tin capsules along with 2.0-4.0 mg of vanadium pentoxide. Capsules were combusted at 1020°C in a Thermo Scientific EA Isolink elemental analyzer. Produced SO 2 was analyzed in a Thermo Scientific Delta V Plus Isotope Ratio Mass Spectrometer via a Conflo IV in continuous flow mode. Sulfur isotope ratios were determined based on the international stand ards IAEA-S1, IAEA-S2, and IAEA-S3 along with an in-house Ag 2 S standard and are reported in delta notation (δ 34 S = [( 34 S/ 32 S) sample /( 34 S/ 32 S) VCDT − 1] * 1,000) relative to Vienna Canyon Diablo Troilite (VCDT). Standard error of Ag 2 S replicates (n = 11 pairs) was 0.08‰ across all sample sets. A second source of error in measuring δ 34 S ZVS values was temperature variability in the heating mantles during CrCl 2 extraction of the ZVS-bearing organic solvents. This is a source of greater imprecision in ZVS extraction than in pyrite extraction because of the temperature-dependent solubility of elemental sulfur . Replicate extractions of sediment samples from several other shallow sediment cores from Chesapeake Bay were performed concurrently with this set of extractions, and the temperature of the heating mantles was varied between 45 and 70°C for a subset of those samples. Resulting measurements of replicate ZVS δ 34 S values from samples subjected to variable extraction temperatures (n = 16 pairs) had a standard error of 2.2‰, which includes the instrumental error of the mass spectrometer.
In the 2021 sampling season, we also measured dissolved pore water sulfide concentrations. Aliquots of 1.40 mL of pore water were pipetted from centrifuged samples and placed into 1.5 mL microcentrifuge tubes, each of which contained 0.10 mL of a 3% Zn-acetate solution. These were stored at 4°C until AVS extractions were 5 of 20 conducted on them. After AVS extraction, resulting Ag 2 S was weighed for abundance and samples with sufficient mass were analyzed for their δ 34 S values in the same manner as the solid phase sulfur samples. For samples from the 2017-2018 sampling dates, we refer to the pore water sulfide concentrations determined via microsensor by Malkin et al. (2022).
After AVS-ZVS-CRS extractions were completed, the remaining solid fraction of sediment was triple-rinsed with DI water via centrifuging and vortexing, then dried at 70°C. These samples were weighed into tin capsules and analyzed for organic carbon concentration and δ 13 C composition on a Thermo Scientific EA Isolink elemental analyzer coupled to a Thermo Scientific Delta V Plus Isotope Ratio Mass Spectrometer via a Conflo IV in continuous flow mode. The AVS-CRS extraction process is not expected to affect organic carbon δ 13 C values (Muller et al., 2017). The standard error of δ 13 C sample replicates (n = 7 pairs) was 0.13‰ across all sample sets. All δ 13 C values are reported relative to the Vienna Pee Dee Belemnite (VPDB) isotopic standard. After centrifuging, drying and homogenization with an agate mortar and pestle, 1-2 g aliquots from the channel site collected in July 2019 were also sent to SGS Canada, Inc. (Lakefield, ON, Canada) and analyzed for solid phase trace metal concentration via inductively coupled plasma mass spectroscopy with standard reference materials OREAS-524, OREAS-135, and OREAS-920. Results for arsenic (As) and iron (Fe) are included in Supporting Information S1.

Radioisotope Analysis
Sediment aliquots were prepared for radioisotope analysis after Schelske et al. (1994) (Schelske et al., 1994). Unsupported 210 Pb activity in each sample was determined by subtracting 226 Ra activity (i.e., supported 210 Pb activity) from total 210 Pb activity at each level. We determined 7 Be (477.6 keV) by gamma spectroscopy. Analyses for 7 Be were completed in less than one half-life (53 days) from the date of core collection. All reported activities were corrected to the date of core collection. Detector efficiencies were established by repeated analyses of a Baltic Sea Sediment standard (IAEA-300; see Ballestra et al., 1994).
We estimated biodiffusive mixing with the stratigraphic profile of excess 210 Pb using the approach of Nittrouer et al. (1984). Sediment mixed zone depth was determined from the depth of penetration of 7 Be. We modeled values for biodiffusive mixing (D b , cm 2 yr −1 ), linear sedimentation rate (cm yr −1 ), and excess 210 Pb activity at time/depth zero (A 0 , dpm g −1 ). These parameters were optimized for mixed zone sediments by minimizing the sum of the negative log-likelihood of modeled and observed pairs for excess 210 Pb activity at the mid-depth of each sediment section. Similarly, values of sedimentation rate for the accumulation zone were optimized by minimizing the sum of the negative log-likelihood of modeled and observed pairs for excess 210 Pb activity at the mid-depth of each sediment section.
We determined recent sediment chronology by 210 Pb dating. We calculated sediment ages by applying the constant-flux, constant-sedimentation (Constant Flux Constant Sedimentation [CFCS]) model (Krishnaswamy et al., 1971) to the accumulation zone samples with unsupported 210 Pb. We estimated the average sedimentation rate as the quotient of the decay constant of 210 Pb (0.031) divided by the slope of the regression of the natural log of unsupported 210 Pb against the mid-depth of each sediment interval. We calculated age errors using the 95% confidence interval of the slope of the regression.

Results
All sediment intervals in this study were 0.5-cm thick increments between 0 and 2 cm: 0.0-0.5 0.5-1.0, 1.0-1.5, and 1.5-2.0 cm. The operationally defined sedimentary sulfur pools analyzed include CRS, AVS, and ZVS (see definitions in Section 2.2). Tests of the solvent extraction-CRS extraction procedure yielded >90% recovery for flower sulfur samples (Acros Organics) below 1.8 mg of total elemental sulfur; therefore, we assume that at least 90% of ZVS was dissolved into organic solvents before CRS extractions, and that the CRS fractions are thus nearly equivalent to the pyrite fractions. The ZVS pool is assumed to comprise primarily S 8 as well as small amounts of polysulfide (J. Liu et al., 2020). The AVS pool is assumed to comprise primarily solid FeS with a negligible amount of pore water sulfide after centrifuging (Luther, 2005).
Solid phase concentrations are expressed as dry weight percent of sulfur or carbon. We describe time series concentration and δ 34 S value trends both as depth-averaged values (i.e., the average value of all depth intervals at one given sampling date at each site) and time-averaged values (i.e., the average of all samples extracted over time from one given depth interval at each site). All data are available in Supporting Information S1.

Sulfur Species Abundance
All solid-phase sulfur species increased in concentration during summer and fall and decreased in concentration during winter and spring. Pyrite predominated over AVS, ZVS, and pore water sulfide at all sampling dates and depths ( Figure 2), and pyrite was more abundant at the shoal site than at the channel site at every date and sampling depth. Between autumn 2017 and February 2018, depth-averaged pyrite sulfur abundance declined 20% at the channel and 21% at the shoal, then increased again through late summer 2018. This pattern appeared to repeat on an annual basis, with lower abundances in winter and spring followed by higher abundances in summer and autumn. Depth-averaged pyrite abundance in the top 2 cm of the sediment at both sites varied over time between 0.52%-0.65% at the channel site and 0.69%-1.01% at the shoal site.
AVS reached a maximum depth-averaged concentration of 0.17% at the shoal site in August 2018, and 0.25% at the channel site in July 2017. Time-averaged AVS concentration at the shoal increased with depth from 0.04% to 0.08%, but at the channel site it declined with depth from 0.10% to 0.07%. The maximum depth-averaged ZVS concentration was 0.09% in June 2017 at the channel site and 0.06% in August 2018 at the shoal site. Time-averaged ZVS concentration was constant with depth at the shoal (0.03%-0.04%), but at the channel site it declined with depth from 0.06% to 0.02%.  (Table 1).

Sulfur Species Isotopic Composition
At both sites, δ 34 S pyr values were more negative and less variable with depth and time than δ 34 S AVS and δ 34 S ZVS values ( Figure 3). The δ 34 S H2S value obtained from August 2021 was also higher than the coeval δ 34 S pyr value (−9.4‰ vs. −18.9‰). Depth-averaged δ 34 S pyr values at both sites featured slightly lower values in late winter to early spring and higher values in late summer (a range of −21.7‰ to −20.3‰ at the shoal site and −23.0‰ to −20.0‰ at the channel site).
At the shoal site, depth-averaged δ 34 S AVS values ranged from −9.5 to −18.9‰ and generally reached higher values in the summer. In contrast, the channel site had a larger range of δ 34 S AVS values (−1.8 to −18.5‰) but did not feature any clear seasonal trends in δ 34 S AVS values. Likewise, depth-averaged δ 34 S ZVS values showed a seasonal trend with lower winter values and higher summer values at the shoal site (total range of −21.3 to −11.4‰) but not at the channel site (total range of −5.7 to −15.2‰).
Time-averaged δ 34 S pyr values at the shoal and channel sites varied by less than 0.5‰ with depth (−21.6 to −21.4‰ at the shoal and −22.6 to −22.2‰ at the channel). Time-averaged δ 34 S AVS values declined with depth at both of the sites, from −14.6 to −16.5‰ at the shoal site and from −9.5 to −11.8‰ at the channel site. Time-averaged δ 34 S ZVS values at both sites decreased with depth: from −16.9 to −17.8‰ at the shoal site, and from −9.3 to −13.6‰ at the channel site. Of the cores collected in 2021, only one sample from the top 2 cm of either core returned a sufficient amount of dissolved pore water sulfide to yield a δ 34 S measurement. This sample was from the interval 1.5-2.0 cm at the shoal site, and its δ 34 S value was −9.4‰.

TOC Abundance and Carbon Isotope Composition
Total organic carbon (TOC) concentration at the shoal site generally exceeded that at the channel site ( Figure S1 in the Supporting Information S1). At the shoal, depth-averaged TOC concentration increased from 2.8% (March 2017) to 3.2% (April 2017), decreased over summer 2017 to 2.9% (August 2017), and later increased to a maximum value of 3.5% in August 2018. At the channel site, average %TOC reached maxima of 2.8% in July 2017 and 2.8% in May 2018. Depth-averaged values of δ 13 C at the channel site ranged from −25.2 to −25.8‰, with no apparent seasonality. Depth-averaged values of δ 13 C at the shoal site ranged from −25.1 to −25.5‰, also with no apparent seasonality.

Radioisotope Profiles
At the shoal site, radiometric data indicated a 4 cm deep mixed layer. Excess 210 Pb at the shoal site varied between 6.6 and 8.7 disintegrations per minute per gram (dpm g −1 ) in the upper 4 cm, then decreased below 4 cm to a minimum value of 1.2 dpm g −1 at 15 cm before dropping below detection ( Figure 4). The depth-integrated excess 210 Pb inventory (25.1 dpm cm −2 ) found at the shoal site was similar to the expected inventory (27.7 dpm cm −2 ) from atmospheric fallout alone (Graustein & Turekian, 1986). The activity of 7 Be at the shoal was greatest (4.4 dpm g −1 ) in the uppermost sample and decreased with depth to a minimum value (1.1 dpm g −1 ) at 4 cm before dropping below detection-likewise implying a mixed layer depth of 4 cm. The 7 Be inventories for both sites (5.9 dpm cm −2 at the shoal and 3.5 dpm cm −2 at the channel) were similar to those previously measured in Chesapeake Bay (Dibb & Rice, 1989).
The constant-flux, constant-sedimentation (CFCS) model of the shoal site yielded a linear sedimentation rate of 0.215 cm yr −1 and a likely range of 0.182-0.245 cm yr −1 , resulting in a date of 1947 at 15 cm depth. The shoal showed a 137 Cs maximum of 1.4 dpm g −1 at 11 cm; the CFCS model estimate for the date of the 137 Cs maximum in this core (11 cm = 1966) showed good agreement with the expected date (1964) for peak 137 Cs fallout in the northern hemisphere. The biodiffusive mixing model ( Figure S2 in Supporting Information S1) indicated that biodiffusion (D b = 5.1 cm 2 yr −1 , and a likely range of 4-7 cm 2 yr −1 ) was similar to the lower range of values found by Dellapenna et al. (1998) in Chesapeake Bay. Our estimate of D b was not sensitive to variability in sedimentation rate within the range of values for sedimentation rate predicted by the CFCS model. The optimization approach employed in the biodiffusive mixing model estimated a sedimentation rate of 0.234 cm yr −1 , similar to the CFCS model estimate.
Radioisotope profiles at the channel site were more ambiguous than those at the shoal site. Excess 210 Pb was only present in the upper 2.5 cm and reached a maximum value of only 0.6 dpm g −1 , while 7 Be showed a random stratigraphic pattern down to 16 cm. Thus, it is possible that the channel site is non-depositional or that the site was dredged at some point. However, dredging is unlikely at a water depth of 24 m because Chesapeake Bay shipping  channels are not dredged below a depth of 16 m (U.S. Army Corps of Engineers, 2005), and the geochemical data we collected does not point to any abrupt changes in sediment geochemistry that would result from dredging within the primary period of sample collection. We consider it unlikely that non-depositional conditions would compromise our interpretations of surficial (<2 cm) pyrite precipitation processes during the course of our sampling.

Effects of Seasonal Dissolution and Bioturbation on the Interpretation of the Pyrite Record
In this study, we report seasonal oscillations in pyrite precipitation and dissolution or burial in sediments below mesohaline waters of Chesapeake Bay, at both a shallow shoal site and at a deep channel site which each experience different levels of bottom water oxygen depletion. Both sites undergo an annual cycle of pyrite precipitation in summer to early autumn followed by dissolution and/or burial in winter and spring (Figure 2), and both sites lose nearly identical fractions of pyrite from summer to winter: pyrite concentration declines by 20% at the channel site and 21% at the shoal site. This matches the 21% loss of pyrite sulfur observed in upper sediments at a point between these two sites at intermediate water depth (12 m) from June 1994 to March 1995 (Cooper & Morse, 1998). Here we discuss the causes of seasonal precipitation-dissolution-burial cycles and how they might affect the pyrite geochemical record.
Pyrite can be removed from surficial sediments through biodiffusive burial into deeper layers and/or through oxidative dissolution. Biodiffusion tends to homogenize surficial signals (Hülse et al., 2022), while oxidative dissolution may partially reset them on an annual basis. Oxygen penetration depths never exceeded 0.4 cm at either site during the sampling interval (Malkin et al., 2022), but other electron acceptors may still contribute to oxidative dissolution of pyrite below this depth (Schippers & Jørgensen, 2002). Biodiffusion can also contribute to oxidative dissolution by mixing electron acceptors deeper down, but at low biodiffusion rates this may be counteracted by downward mixing of organic matter (S. van de Velde and Meysman, 2016). On the other hand, if strong biodiffusion is accompanied by strong bioirrigation, the addition of oxygenated water may be expected to promote further oxidative dissolution. Thus, strong biodiffusion paired with strong bioirrigation will probably promote oxidative dissolution throughout the sedimentary mixed layer, but weak biodiffusion is less likely to do so.
At the shoal site, both long-lived excess 210 Pb and short-lived 7 Be indicated moderate biodiffusion rates compared to other sites in Chesapeake Bay (Dellapenna et al., 1998) and globally compiled sites at similar water depths (Lecroart et al., 2010), as expected for a seasonally hypoxic site. With regard to oxidative dissolution, we note that the shallowest interval at the shoal site had more pyrite than underlying layers in summer but less pyrite in winter. This suggests that oxidative dissolution contributed to pyrite removal at the shoal site because biodiffusion by itself would only homogenize pyrite concentrations but would not bring surficial pyrite concentration below that of the underlying layers. Therefore, oxidative dissolution is likely the primary process promoting pyrite removal at the shoal site in the winter, but there is probably a minor contribution from biodiffusive burial.
The channel site is more difficult to characterize because radioisotope profiles did not enable a robust estimate of biodiffusion rates. However, the channel site is severely oxygen-stressed and even under oxic winter conditions it features very limited infaunal colonization, primarily by small oligochaetes (Marvin-DiPasquale & Capone, 1998). Therefore, pyrite removal at the channel site in winter is more likely to result from oxidative dissolution than from biodiffusive burial. The fact that both sites lose nearly identical fractions of surficial pyrite from summer to winter is interesting and somewhat unexpected in light of the contrasts in bioturbation intensity between the two sites. Because the shoal site accumulates more total pyrite in summer, it also loses more total pyrite in winter; therefore, it may be reasonable to conclude that oxidative dissolution under oxic wintertime conditions at both sites (Table 1) removes similar absolute amounts of surficial pyrite, but that additional pyrite removal at the shoal results from biodiffusive burial.
Both of the studied sites illustrate the axiom that the sedimentary pyrite record must primarily reflect net-precipitation intervals. In other words, trace metal compositions and δ 34 S values of pyrite at these sites will reflect environmental conditions that prevail in summer and autumn, not during winter. By this mechanism, pyrites from ancient settings that experienced seasonal changes in geochemistry may reflect a "summertime signal" imparted during warm, sulfidic net-precipitation periods. Intervals of net dissolution may also alter the existing geochemical record if metals or sulfur isotopes are distributed unevenly in the pyrite-for example, by removing certain trace metals such as Cu and Ni that can be more concentrated in overgrowths of pyrite framboids than in the framboids themselves (Gregory et al., 2022). Oxidative pyrite dissolution has been shown to fractionate pyrite sulfur and product sulfate, typically by <1‰ (Heidel & Tichomirowa, 2011;Pisapia et al., 2007;Taylor et al., 1984;Toran & Harris, 1989), which may contribute to the relatively minor variations in pyrite δ 34 S values across seasons.

Apparent Disparity in Pyrite Formation Rates Between Sites
Pyrite concentration at the shoal site consistently exceeded pyrite concentration at the channel site in spite of the higher summertime concentration of shallow pore water sulfide at the channel site (Table 1). In fact, the correlation between pyrite concentration and precursor concentration was more positive at the shoal site for all three precursor pools (AVS, ZVS, and H 2 S) as well as for TOC (Figure 5a-5d). Pearson correlation coefficients and p-values for concentrations of pyrite versus its precursors at the shoal site were 0.69 for AVS (p = 0.0003, n = 23), 0.76 for ZVS (p = 0.00003, n = 23), and 0.66 for log 10 (H 2 S) (p = 0.0006, n = 23), versus 0.34 (p = 0.05, n = 35), 0.13 (p = 0.46, n = 35), and 0.30 (p = 0.08, n = 35) at the channel site, respectively. Pyrite sulfur concentrations typically increase with TOC concentration in marine sediments in a C/S ratio of 2.8 ± 0.8, which has led to the suggestion that TOC concentration is a key determinant of pyrite sulfur concentration when sulfate is not limiting (Morse & Berner, 1995). However, this C/S relationship did not hold true at the channel site, where TOC and pyrite had a slightly negative Pearson correlation of −0.17, versus 0.58 at the shoal site. While high TOC concentration (7% or more) may cause organic compounds to scavenge reduced sulfur more quickly than FeS (Shawar et al., 2018), the C/S slope remained positive at the shoal site even though bulk TOC concentration was generally higher there. It is possible that the organic matter of channel sediments is more sulfide-reactive than that of shoal sediments, but both sites have lower TOC concentrations than the sites studied in Shawar et al. (2018). We infer that the lower Pearson coefficients for pyrite concentration versus precursor concentration at the channel site result from slower pyrite precipitation in these sediments. In this section we explore the possible drivers of this disparity, including SOM community dynamics.
Certain factors influence the concentrations of pyrite and its reactant species, but these factors do not explain the apparent difference in reaction rates of FeS with H 2 S/ZVS between the two sites. For example, finer sediment grain size can increase sulfate reduction rates (Hemingway et al., 2019;Lin et al., 2022) and sediments near the shoal site tend to be finer than channel sediments-clays/silty clays versus silty clays, respectively (Kerhin et al., 1988). Despite this, H 2 S accumulates to much higher concentrations in channel sediments, yet this is not correlated with higher pyrite concentrations there ( Figure 5). It is unknown whether sediment grain size can affect reaction rates between FeS, ZVS, and H 2 S. Considering pyrite dissolution rather than precipitation, it is possible that the more permeable channel sediments may encourage pyrite dissolution by slowing the downward diffusion of oxidants. However, both sites lose equal fractions of pyrite in winter, meaning that greater oxidative dissolution in channel sediments resulting from larger grain size would need to occur concurrently with net precipitation in summer and early fall, when the channel site's water column is strongly hypoxic.
Similarly, while Fe/Mn-oxides are likely to attain higher concentrations at the shoal site (Sinex & Helz, 1981), it is unknown whether Fe/Mn-oxides can affect the direct reaction between FeS and H 2 S/ZVS. MnO 2 can oxidize pyrite sulfur to sulfate and FeS sulfur to elemental sulfur (Schippers & Jørgensen, 2001), and Fe-oxides can oxidize H 2 S to elemental sulfur (Poulton et al., 2004;Pyzik & Sommer, 1981) but cannot directly oxidize pyrite (Schippers & Jørgensen, 2002). Nitrate may also oxidize FeS or dissolved sulfide (Schippers & Jørgensen, 2002). The lower benthic salinity of the shoal site (∼13 vs. ∼19) may increase the rate of H 2 S oxidation by Fe-oxides at the shoal site compared to the channel site (Poulton et al., 2004), which could partially explain the lower H 2 S concentrations there. Although these reactions affect the concentrations of FeS, ZVS, and H 2 S, they do not account for why pyrite concentration is more positively correlated with FeS, ZVS, and H 2 S concentrations at the shoal site. Direct pyrite oxidation by MnO 2 is more likely at the more oxic shoal site, but this reaction should lower the pyrite/reactant slope, not increase it.
What might drive the apparent contrast in the kinetics of pyrite precipitation between channel and shoal? One possibility is that new pyrite nucleated more easily at the shoal site; another is that post-nucleation crystal growth was more favorable there. Under abiotic laboratory conditions, pyrite nucleation requires strong supersaturation (D. Rickard, 2012) and in some cases is regarded as rate-limiting compared to subsequent crystal growth (Schoonen & Barnes, 1991). Other studies, however, have posited that the presence of biological material (Grimes et al., 2001) or preexisting pyrite crystals (Harmandas et al., 1998) may substantially lower the nucleation barrier. Rates of pyrite nucleation and subsequent growth can also be affected by trace metal inventories. For example, at a trace element to Fe molar ratio of 0.1%-0.5%, Mn and Ni decrease the time required for pyrite nucleation by 60% and 80%-99%, respectively (Baya et al., 2021(Baya et al., , 2022. The western shoal in this area generally has more Mn and more Ni (∼1.5% Mn/Fe and ∼0.09% Ni/Fe) than the channel does (∼0.6% Mn/Fe and ∼0.04% Ni/Fe), which accordingly may induce faster pyrite nucleation in shoal sediments (Sinex & Helz, 1981). On the other hand, as previously noted, MnO 2 can dissolve pyrite (Schippers & Jørgensen, 2001), which means that without knowing the net oxidation state of manganese in these sediments, it is difficult to say whether Mn at the shoal site increased the pyrite precipitation rate, the pyrite dissolution rate, or both.
The activity of SOM may alter the concentrations of certain trace metals that affect pyrite precipitation rates. For example, arsenic slows down pyrite nucleation and crystal growth; addition of As to 0.1%-0.5% of Fe increases pyrite nucleation time by 70%-140%, and As also increases the time required for pyrite to become the dominant Fe-S mineral in its batch from 14 days to over 4 months (Baya et al., 2021(Baya et al., , 2022). An As-induced lag time of 4 months could strongly influence our study sites, where net surficial pyrite accumulation is limited to seasonal timescales. Although comparative shoal-channel data are not available for As, cable bacteria may accumulate As in surficial sediments of the channel site. Active cable bacteria colonization promotes adsorption of As onto Fe-oxides in surficial sediments, attaining As/Fe molar ratios as high as ∼0.2% until anoxic conditions release the As into bottom waters (S. J. van de Velde et al., 2022). Surface sediments from the channel site in July 2019 had an As/Fe ratio of 0.03% (Supporting Information S1), and although this ratio is lower than the As/Fe ratios that were experimentally demonstrated to impede pyrite nucleation, this measurement was taken during the apex of summertime anoxia in the Chesapeake Bay channel, by which point most of the As accumulated by cable bacteria colonization already would have dissolved. The higher density of cable bacteria documented at the channel site in springtime relative to summer-and at the channel relative to the shoal-suggests that As may accumulate to higher concentrations at the channel site than at the shoal site in spring and early summer. In this way, cable bacteria would delay seasonal pyrite nucleation and impede its growth in the channel by concentrating As in surficial sediments in spring and early summer. This process would be analogous to the seasonal "firewall" against sulfide release from sediments that cable bacteria provide by generating shallow sedimentary Fe-oxides (Seitaj et al., 2015).
Beggiatoa may alternately quicken pyrite precipitation by making more reactive ZVS compounds. Elemental sulfur (S 8 ) can exist in dissolved, colloidal, or crystalline form, and the timescale of its reaction with the sulfide-polysulfide pool can range from seconds in dissolved form (Kamyshny et al., 2003) to minutes or days in colloidal form (Fossing & Jørgensen, 1990;Kamyshny & Ferdelman, 2010) to more than a year in crystalline form when polysulfide is absent (Avetisyan et al., 2019). The form that S 8 takes is strongly influenced by Beggiatoa, which store large amounts of intracellular S 8 mixed with polysulfide (Berg et al., 2014). The S 8 stored in Beggiatoa is also frequently encased in sulfur globule proteins (Schmidt et al., 1986;Strohl et al., 1981) that render the biogenic S 8 more hydrophilic and make it more sulfide-reactive than inorganic S 8 (Dahl, 2020;Kleinjan et al., 2005a). Because Beggiatoa predominate over cable bacteria at the shoal site but not at the channel site (Malkin et al., 2022), they may generate soluble, pyrite-reactive S 8 in larger quantities at the shoal site.
The relationship between pyrite concentrations and the predominance of Beggiatoa over cable bacteria (quantified as the difference between their biovolume; Figure 5e) highlights the possibility that contrasting SOM communities might influence rates of pyrite precipitation. It shows that samples with a greater predominance of Beggiatoa over cable bacteria (quantified as Beggiatoa biovolume minus cable bacteria biovolume) tend to have higher concentrations of pyrite. It should be noted, however, that many of the Beggiatoa-dominated samples are from the shoal site while many of the cable bacteria-dominated samples are from the channel site, so it is difficult to state with certainty whether a predominance of Beggiatoa over cable bacteria actually promotes pyrite precipitation or whether other disparities between the two sites are influencing both SOM biovolume and pyrite precipitation rates. Nevertheless, the correlations documented in Figure 5 point to an avenue for further research into the roles of SOM in influencing rates of sedimentary pyrite formation directly from pyrite-forming sulfur compounds, in addition to the broader impacts of microbes in general (Picard et al., 2016(Picard et al., , 2018Wacey et al., 2015).

Sulfur Isotope Signatures of Pyrite Compared to Reactants
Despite much wider variation in the coeval δ 34 S values of AVS and ZVS, depth-averaged δ 34 S values of pyrite varied by only 3‰ over time at the channel site and only 2‰ over time at the shoal site (Figure 3). Surficial pyrite δ 34 S values at these sites thus appear insensitive to seasonal fluctuations in biogeochemical conditions. It is known that sedimentation rates and physical reworking can influence pyrite sulfur isotope compositions (Aller et al., 2010;X. Liu et al., 2019;Pasquier, Bryant, et al., 2021), in part by setting the ratio of advective to diffusive sulfate flux (Claypool, 2004). Radioisotope results indicate that the shoal site has a sedimentation rate of ∼0.2 cm yr −1 while the channel site may be net non-depositional, yet the difference between average δ 34 S pyr values at the two sites is only 1.0‰ (−22.4‰ in the channel vs. −21.4‰ on the shoal for all samples). Part of this insensitivity likely results from the lack of net pyrite formation in winter and spring, which decreases the likelihood that differences between the sites in winter and spring play a role in determining δ 34 S values. It may also be the case that different sedimentation rates have not yet induced significant Rayleigh distillation effects in surficial sediments that are still in relatively close contact with the water column.
While pyrite is the most thermodynamically stable iron sulfide mineral at Earth surface conditions, its formation is kinetically limited by nucleation (D. Rickard, 2012). Variations in microbial sulfate reduction rate, reoxidative sulfur cycling, system openness, and sedimentary redox conditions will strongly influence the isotopic composition of all of the pyrite-forming sulfur compounds discussed here. However, it has typically been assumed that there is no net fractionation of sulfur isotopes associated with the direct reactions of FeS with H 2 S and FeS with S n 2− ; these reactions have instead been modeled in terms of simple isotopic mass balance based solely on the stoichiometry of combination or substitution reactions (Butler et al., 2004). Yet even during net-precipitation intervals, δ 34 S values of pyrite at our study sites reveal another unusual feature: not only do they remain stable over time, they also remain lower than those of coeval precursor compounds. For example, in Aug. 2017, when pyrite concentration reached its annual peak, mean δ 34 S values of AVS, ZVS, and pyrite at the shoal site were −14.3, −15.7, and −21.4‰ respectively. The next summer (Aug. 2018), AVS and ZVS at the shoal site reached higher δ 34 S values of −9.7 and −10.8‰, but the δ 34 S value of pyrite remained low and stable at −21.1‰. In fact, pyrite δ 34 S values were lower than any of the analyzed pyrite precursor compound δ 34 S values (i.e., AVS, ZVS, and one recovered value for H 2 S) at each sampling date and depth. This raises a question posed by Raven et al. (2016), who measured 34 S-depleted pyrite in marine sediments of the Santa Barbara Basin: how does sedimentary pyrite precipitate with δ 34 S values that are lower than those of all coexisting precursor compounds?
This apparent deviation from mass balance has been observed in other shallow marine and estuarine sediments (Kaplan et al., 1963;Raven et al., 2016;Strauss et al., 2012;Zopfi et al., 2008). For example, in Beggiatoa-colonized surface sediments of the Bay of Concepción, Chile, pyrite δ 34 S values were 2‰-7‰ lower than coeval δ 34 S ZVS values and 2‰-13‰ lower than coeval δ 34 S AVS values, while at a nearby Thioploca-colonized site, surficial δ 34 S pyr values were 2‰-4‰ lower than δ 34 S ZVS values and 4‰-6‰ lower than δ 34 S AVS values (Zopfi et al., 2008). Anomalous 34 S depletion generally does not exceed a few per mil in shallow (<10 cm) sediments, but it sometimes increases in deeper sediments (e.g., Brüchert et al., 1995;Raven et al., 2016). This increase with depth been attributed to the δ 34 S value that pyrite inherits from rapid formation by the polysulfide reaction in shallow sediments (J. Liu et al., 2020). However, this explanation still does not account for very shallow sedimentary pyrites that are 34 S-depleted relative to all coexisting reactants.
An alternate explanation invokes pyrite formation in sedimentary microenvironments and isotope exchange between pore water sulfide and other sulfur compounds (Raven et al., 2016). Slow infilling of pyrite in microenvironments followed by later precipitation of pyrite reactant compounds may create a reactant-product δ 34 S offset, especially in Fe-limited settings, though it is difficult to conclusively document these micro-scale processes in sediments. HCl-extractable iron persists in mid-Bay surficial sediments at concentrations of 100-200 μmol g −1 in early summer and the degree of pyritization of Fe above 2 cm depth in the mid-Bay is between 25% and 50% (Cornwell & Sampou, 1995). This suggests that Fe is not limiting in surficial mid-Bay sediments, although the wide range of reaction timescales of HCl-extractable Fe (days to >10 4 yr) makes it challenging to determine the extent of Fe-limitation (Raiswell et al., 1994).
Because δ 34 S values for H 2 S are not present for most of our sampling period, we cannot fully rule out the possibility that this apparent mass imbalance is actually resolved by low δ 34 S values in surficial H 2 S. However, pyrite formation via the H 2 S pathway results from a 1:1 combination of FeS sulfur and H 2 S sulfur (Butler et al., 2004), which would require summer δ 34 S H2S values between −24 and −32‰ to balance the relatively high δ 34 S values of AVS compared to produced pyrite. By comparison, pore water sulfide at the shoal site had a measured δ 34 S value of −9.4‰ in August 2021. Summertime δ 34 S values as low as −24 to −32‰ for H 2 S are therefore unlikely, especially considering that the solid AVS pool is reset each spring by dissolution ( Figure 2) and that the summertime δ 34 S value of surficial FeS is likely to be similar to that of the coexisting dissolved sulfide pool from which it precipitates (Böttcher et al., 1998). This makes the H 2 S pathway improbable from an isotope mass balance perspective.
Even if the H 2 S pathway appears to be an unlikely source of the pyrite in these sediments, the polysulfide pathway must still account for a 5-10‰ gap between pyrite and coeval ZVS. We propose that the anomalous δ 34 S offset between reactant ZVS and product pyrite in shallow sediments can result from position-specific isotope fractionation (PSIF) in the polysulfide pool. Although PSIF effects have been implicated in polysulfide equilibrium isotope distributions (Amrani et al., 2006), it is not clear how PSIF would impact sulfur isotope distributions in natural sediments. Here we develop a simple model of equilibrium isotopic fractionation between S 8 , polysulfide, and pyrite that arises from position-specific effects. This steady-state model solves a system of algebraic equations based on observational data of S 8 and polysulfide δ 34 S values from two previous studies (Amrani & Aizenshtat, 2004;Amrani et al., 2006). The first modeled system component is the equilibrium between S 8 , sulfide, and polysulfide. Polysulfides are highly reactive chains of sulfur atoms with the formula S n 2− , where 2 ≤ n ≤ 9 (D. Rickard & Luther, 2007). Multiple chain lengths exist in equilibrium with each other, with n = 4-6 forming the bulk of the polysulfide pool (Kamyshny et al., 2004). Polysulfides equilibrate with each other and with dissolved sulfide and dissolved S 8 on a timescale of seconds (Kamyshny et al., 2003). This equilibrium involves the disproportionation of polysulfide into elemental sulfur and sulfide: Generally, at lower pH, the system's equilibrium favors polysulfide disproportionation to S 8 and H 2 S, while at higher pH the system's equilibrium favors polysulfide and HS − . In Amrani and Aizenshtat (2004), S 8 produced from complete disproportionation of the polysulfide pool had a δ 34 S value 4.4‰ ± 1.0‰ higher than that of the produced sulfide; in the model, we refer to this isotopic offset as the variable ε-S8-HS .
The second component of the model is based on equilibrium fractionation within the polysulfide pool observed by Amrani et al. (2006), who determined δ 34 S values for S 4 2− , S 5 2− , S 6 2− , and S 7 2− in a solution that was undersaturated with respect to S 8 . The authors found that long-chain polysulfides had higher δ 34 S values than short-chain polysulfides and concluded that this effect resulted from fractionation of 34 S toward the middle of each chain, which made long chains with internal 34 S more stable than those without it. The δ 34 S values recorded for S 4 2− through S 7 2− in Amrani et al. (2006) were, in ascending order: 16.9‰, 18.1‰, 19.2‰, and 20.6‰ VCDT (1σ = 0.3‰ for all values) in a system whose total dissolved reduced sulfur (TDRS, defined as the sum of dissolved S 8 , dissolved sulfide, and polysulfide) had a δ 34 S value of 16.5‰. The remaining HS − in the system was calculated to have a δ 34 S value of 13.8‰, and linear extrapolation of the four measured points predicts δ 34 S values of 21.8‰ and 23.0‰ for S 8 2− and S 9 2− . The bulk δ 34 S value of TDRS does not affect the model solution; as such, all model input and output values are expressed relative to TDRS, not on the VCDT scale.
It should be noted that the two published data sources used for this model had very different S 8 saturation states from those of surficial estuarine sediments, and by definition S 8 would not precipitate from the undersaturated solution used in Amrani et al. (2006). However, under ambient sediment conditions at circumneutral pH, S 8 will be much closer to saturation than it was in the experiments that informed our model inputs (Avetisyan et al., 2019). At near-saturation conditions, biogenic S 8 -which is likely to be in colloidal form )may remain stable long enough to constitute a significant part of the ZVS pool while also approaching isotopic equilibrium with pyrite-forming polysulfides on a timescale of minutes to days (Fossing & Jørgensen, 1990;Kamyshny & Ferdelman, 2010).
For each chain length n, we model the δ 34 S values of a pool of S n 2− molecules as the average of the δ 34 S values of all the positions in that chain length. For example, a pool of S 4 2− whose terminal sulfurs have a δ 34 S value of 5.0‰ and whose central sulfurs have a δ 34 S value of 8.0‰ will have a bulk δ 34 S value of 6.5‰. The equations for the δ 34 S values of each polysulfide and its component positions follow the form where m = 1 for odd values of n and m = 2 for even values of n.
The terminal position in all polysulfides is notated as a, the second position inward in all polysulfides is notated as b, the third inward as c, the fourth inward as d, and the fifth inward (i.e., the central sulfur in S 9 2− ) is notated as e (Figure 6a). The variable q is the amount by which the δ 34 S value of the terminal position in S 4 2− is offset from the δ 34 S value of the system's TDRS; we assume that q has a range of values between that of the bisulfide pool (−2.7‰ relative to TDRS) and that of the measured S 4 2− (+0.3‰ relative to TDRS) from Amrani et al. (2006) ( Table 2). The variable p is the per mil amount by which the δ 34 S value of a sulfur in a given position in a chain of length n + 1 exceeds that of the sulfur in the same position in chain length n; for example, the c position in S 8 2− is p per mil higher than the c position in S 7 2− . We assume that p increases linearly with chain length, which simplifies the model and fits the linear increase in δ 34 S values versus chain length (R 2 = 0.998) measured by Amrani et al. (2006). The variable ϵ b is the per mil δ 34 S offset between position a and b in all chain lengths, while f c , f d , and f e are coefficients that modify ϵ b between positions b-c, c-d, and d-e. Because polysulfide disproportionation (Equation 1) involves expulsion of the terminal sulfur (Steudel, 1996), it is likely that ϵ b has a similar value to ε-S8-HS . Any pair of input values q and p generates values of ϵ b , f c , f d , and f e that can be used to calculate the δ 34 S value of each position in each polysulfide (see Equations S11-S14 in the Supporting Information S1).
The final model component is the polysulfide-pyrite reaction, an exchange reaction in which two adjacent sulfurs at one end of a polysulfide chain are incorporated into the new pyrite molecule while the sulfur from the reacting FeS is exchanged and added back to the polysulfide (Butler et al., 2004;Luther, 1991). The reaction for the dominant polysulfide chain length, S 5 2− , is Our baseline model calculates the total δ 34 S offset between pyrite and S 8 as the difference between (a) the δ 34 S value of S 8 formed from S 9 2− disproportionation and (b) the δ 34 S value of pyrite formed from the most abundant polysulfide, S 5 2− (Kamyshny et al., 2004), with the pyrite and the S 8 thus derived from the same polysulfide pool. The model with all p-q combinations specified above yielded δ 34 S S6 values that match the measured values in Amrani et al. (2006) within instrumental precision (see Supporting Information S1 for further details on model validation and sensitivity tests). We compared the p-q parameter space consistent with the results of Amrani et al. (2006) and used the zone of overlap to determine that the mostly likely values of the S 8 -pyrite δ 34 S offset range from 5.5‰ to 6.4‰, but can vary between 4.3‰ and 7.3‰ ( Figures S3 and S4 in the Supporting Information S1).
Can the modeled S 8 -pyrite offsets of 4.4‰-7.3‰ explain the observed gaps between δ 34 S ZVS and δ 34 S pyr values in Chesapeake Bay sediments? At the shoal site, the only time when the offset between δ 34 S ZVS and δ 34 S pyr values exceeded this range was in August 2018, when it reached a mean value of 10.2‰. Mass balance calculation of the δ 34 S value of "new pyrite" added from Febuary 2018 to August 2018 indicates that the freshly formed pyrite would have had a higher δ 34 S value than the preexisting pool, requiring a smaller offset of 8.5‰. The record at the channel site is somewhat more ambiguous. At sampling dates of net pyrite precipitation-i.e., dates when depth-averaged pyrite concentration increased relative to the previous sampling-the difference between depth-averaged δ 34 S ZVS and δ 34 S pyr values at the channel site ranged from 5.9‰ (August 2018) to 11.9‰ (July 2017). The August 2018 value represents the only net-precipitation interval at that site in 2018, and it lies within the model's predicted S 8 -pyrite offset range. In 2017 there is one net-precipitation interval when the offset is nearly 5‰ larger than the largest modeled offset of 7.0‰, but ZVS-pyrite offsets from the other three net-precipitation  -FeS exchange reaction that forms pyrite, with the same two sulfur atoms again colored red. (b) For the shoal site, average δ 34 S values of acid volatile sulfur (AVS), H 2 S, zero-valent sulfur (ZVS), and pyrite measured during intervals of net pyrite precipitation are plotted on the left side. On the right side are pyrite δ 34 S values generated by modeling three different reaction pathways: green squares represent the H 2 S pathway, blue diamonds represent the ZVS pathway without position-specific isotope fractionation (PSIF) effects, and red gradients represent the histogram distribution (see Figure S4 in the Supporting Information S1) of S 8 -pyrite δ 34 S offsets ranging from 4.3‰ to 7.3‰. (c) This panel contains the same comparisons for the channel site; the H 2 S pathway model assumes that δ 34 S H2S values equal δ 34 S AVS values. intervals in 2017 fall in the range of 6.5-8.5‰. Interpretation of these four intervals at the channel site in 2017 is somewhat obfuscated by the small successive increases in depth-averaged pyrite concentration from 1 month to the next, each of which is much smaller than instrument precision. Despite the limitations imposed by measurement precision, we conclude that the PSIF model can account for much, if not all, of the apparent isotopic mass imbalance in surficial pyrite at these two sites (Figures 6b and 6c). This hypothesis does not rule out other processes such as precipitation in biomass-associated microenvironments (Raven et al., 2016) that may also contribute to 34 S depletion in pyrite relative to precursor compounds.
Our model assumes equilibrium between polysulfides and S 8 , but the timescale of equilibration between these species varies between seconds for dissolved S 8 (Kamyshny et al., 2003) and minutes to days for colloidal S 8 (Fossing & Jørgensen, 1990;Kamyshny & Ferdelman, 2010). The model's ability to predict ZVS-pyrite δ 34 S offsets in natural sediments depends on whether the sedimentary S 8 pool equilibrates with the polysulfide pool quickly enough to impact the pyrite formation reaction (Equation 3). Equilibration rates depend on S 8 saturation and S 8 solubility, which depend on the form that S 8 takes (dissolved, colloidal, or crystalline), the ambient pH, the concentration of sulfide and polysulfide, and the presence or absence of S 8 -encasing proteins (Kafantaris & Druschel, 2020;Kamyshny et al., 2004;Kleinjan et al., 2005b;Maki, 2013). The presence of polysulfide in intracellular S 8 globules in Beggiatoa (Berg et al., 2014), the hydrophilicity of protein membranes that typically coat these globules (Kleinjan et al., 2005a), and the high colloidal fraction of intracellular S 8 in Beggiatoa ) all point to a soluble and reactive pool of S 8 that is likely to equilibrate quickly with the ambient sulfide-polysulfide pool.
Finally, we also consider the possibility that intracellular cycling between S 8 and polysulfide could increase ZVS-pyrite S isotope fractionation beyond the modeled range via repeated S 8 cleavage and S 9 2− disproportionation as Beggiatoa consume stored intracellular S 8 . For example, if two Beggiatoa cells initially store intracellular ZVS with the same δ 34 S value but the first cell ruptures, its ZVS can begin forming pyrite and will move toward an equilibrium isotopic distribution with the ambient polysulfide pool at a rate determined by its solubility. Meanwhile, if the second cell migrates toward a sulfide-poor region and begins oxidizing part of its S 8 stock (Berg et al., 2014), it will cleave its S 8 to use polysulfide as a labile redox intermediary. If some of this polysulfide is reconstituted back into S 8 via Equation 1 at a later point, the reprecipitation may concentrate more 34 S into the new S 8 molecules. This mixture of extracellular ZVS and 34 S-enriched intracellular S 8 could widen the gap between measured δ 34 S ZVS and δ 34 S pyr .

Conclusions
In this study we compared seasonal changes in sulfur compound concentrations and sulfur isotope distributions between two sites in Chesapeake Bay: a shallow, bioturbated site on the Bay's western shoal (CB4.3W) and a deeper, defaunated site in the Bay's central channel (CB4.3C). The strong seasonality, low bioturbation rates, and oxygen-stressed conditions at these sites may offer an analog for environmental conditions in Earth's pastparticularly the Early Paleozoic, when global oxygen concentrations and bioturbation rates were lower than today (Tarhan et al., 2015;Tostevin & Mills, 2020). For this reason, seasonally dynamic, oxygen-stressed benthic modern settings may be useful settings to explore the benthic biogeochemistry of the Earth's past.
Radioisotope profiles indicated a mixed layer of 4 cm and a biodiffusion rate of 5.1 cm 2 yr −1 at the shoal site, but also suggested that the channel site may be a non-depositional setting. Both sites experienced net pyrite precipitation in summer and the loss of ∼20% of surficial pyrite in winter. The bioturbated shoal site consistently had higher concentrations of pyrite than the more oxygen-depleted and depauperate channel site, and changes in pyrite concentration at the shoal site appeared to be more responsive to the concentrations of measured pyrite precursors (FeS, ZVS, and pore water sulfide) than at the channel site. We hypothesize that this may stem from the different SOM communities that dominate each site (Figure 7). The Beggiatoa that predominate at the shoal site may provide ZVS in a more soluble, reactive form than the ZVS at the channel site; alternately, the cable bacteria that predominate at the channel site may delay the seasonal growth of pyrite by accumulating arsenic in surface sediments in spring and early summer (Baya et al., 2021; in a process analogous to the seasonal delay of sulfide release from sediments (Seitaj et al., 2015).
We determined that the polysulfide reaction pathway was probably forming much of the pyrite in the surficial sediments. Future work may better differentiate the H 2 S and polysulfide reaction pathways by targeting the distinct isotopic fingerprints of these reactions. For example, sediment geochemical models could include a more fully resolved polysulfide pool and the sulfur exchange reaction between polysulfide and FeS, although the influence of microbial sulfur metabolisms still poses a challenge in modeling reaction rates (Canfield et al., 1998). Alternately, microanalytical techniques may yield useful data by combining high-resolution morphological analysis with targeted isotopic analysis of different pyrite morphologies (Bryant et al., 2020;Cui et al., 2018;Wacey et al., 2010).
Surficial pyrite δ 34 S values at both sites remained stable over time despite oscillating environmental conditions and despite much more variability in the δ 34 S values of coexisting pyrite precursor compounds. Pyrite δ 34 S values were consistently more negative than those of their coexisting precursors, which has been observed in other settings (Kaplan et al., 1963;Raven et al., 2016;Strauss et al., 2012;Zopfi et al., 2008). We used a steady-state model to evaluate whether PSIF effects in polysulfides can explain the observed gap between δ 34 S values of pyrite and the ZVS pool. Based on previously published measurements of δ 34 S values in the S 8 -polysulfide-sulfide system, we estimate that this effect can generate pyrite δ 34 S values that are up to 7‰ lower than ambient dissolved/colloidal S 8 . This effect would explain much of the anomalous 34 S-depletion of surficial pyrite relative to ZVS at the studied sites ( Figure 6). The uncertain expression of PSIF effects in systems close to S 8 saturation and the potentially wide range of S 8 solubility in SOM-colonized sediments suggest that further investigation of PSIF effects in polysulfides is needed to corroborate these results. Nevertheless, the very rapid (∼10 s) equilibration time of polysulfides raises the possibility that the isotopic imprint of PSIF in polysulfides is a common feature in pyrites that have formed via the polysulfide reaction pathway.
Finally, we note that the similarity in pyrite δ 34 S values between these two sites offers a useful benchmark in interpreting pyrite δ 34 S signals in ancient sediments, as it appears that changes in SOM ecology-at least in regard to the balance of cable bacteria versus Beggiatoa-do not appear to affect the δ 34 S record of surficial pyrite. On the other hand, our data point to the likelihood that shifts in benthic SOM communities can influence rates of shallow pyrite precipitation. Additionally, the "summertime signal" inherent to seasonal pyrite precipitation-dissolution cycles may need to be taken into consideration when analyzing ancient sediments that were subject to periodic seasonal variation in geochemical conditions.

Data Availability Statement
The geochemical data and the PSIF model code used for interpretation of pyrite formation and sediment dynamics in this study are available at the Johns Hopkins Research Data Repository (https://archive.data.jhu.edu) via https://doi.org/10.7281/T1/EHSEV2. are schematics that illustrate how annual cycles of pyrite precipitation and dissolution compare to increasing and decreasing activity of Beggiatoa, cable bacteria, and bioturbating fauna. Hypothesized SOM-induced mechanisms for promoting faster or slower pyrite precipitation (reactive S 8 and arsenic buildup, respectively) are indicated in these panels. In all three panels, the white zone represents periods of net pyrite precipitation and the gray zone represents intervals of net dissolution.