The Formation of Highly Positive δ34S Values in Late Devonian Mudstones: Microscale Analysis of Pyrite (δ34S) and Barite (δ34S, δ18O) in the Canol Formation (Selwyn Basin, Canada)

The sulfur isotope composition of pyrite in marine sedimentary rocks is often difficult to interpret due to a lack of precise isotopic constraints for coeval sulfate. This study examines pyrite and barite in the Late Devonian Canol Formation (Selwyn Basin, Canada), which provides an archive of δ34S and δ18O values during diagenesis. Scanning electron microscopy (SEM) has been combined with microscale secondary ion mass spectrometry (SIMS) analysis (n = 1,032) of pyrite (δ34S) and barite (δ34S and δ18O) on samples collected from nine stratigraphic sections of the Canol Formation. Two paragenetic stages of pyrite and barite formation have been distinguished, both replaced by barium carbonate and feldspar. The δ34Sbarite and δ18Obarite values from all sections overlap, between +37.1‰ and +67.9‰ (median = +45.7‰) and +8.8‰ and +23.9‰ (median = +20.0‰), respectively. Barite morphologies and isotopic values are consistent with precipitation from diagenetically modified porewater sulfate (sulfate resupply << sulfate depletion) during early diagenesis. The two pyrite generations (Py-1 and Py-2) preserve distinct textures and end-member isotopic records. There is a large offset from coeval Late Devonian seawater sulfate in the δ34Spyrite values of framboidal pyrite (-29.4‰ to -9.3‰), consistent with dissimilatory microbial sulfate reduction (MSR) during early diagenesis. The Py-2 is in textural equilibrium with barite generation 2 (Brt-2) and records a broad range of more positive δ34SPy-2 values (+9.4‰ to + 44.5‰). The distinctive highly positive δ34Spyrite values developed from sulfate limited conditions around the sulfate methane transition zone (SMTZ). We propose that a combination of factors, including low sulfate concentrations, MSR, and sulfate reduction coupled to anaerobic oxidation of methane (SR-AOM), led to the formation of highly positive δ34Spyrite and δ34Sbarite values in the Canol Formation. The presence of highly positive δ34Spyrite values in other Late Devonian sedimentary units indicate that diagenetic pyrite formation at the SMTZ may be a more general feature of other Lower Paleozoic basins.


INTRODUCTION
Methane is a powerful greenhouse gas that is produced during the final stage of organic matter fermentation (Knittel and Boetius, 2009). Changes in the flux of methane (from sediment to oceans) have been linked with major climatic impacts at particular stages of earth history (e.g., Dickens et al., 1995). In modern ocean sediments, sulfate reduction coupled with the anaerobic oxidation of methane (SR-AOM) accounts for approximately 80% of methane oxidation, thereby regulating the release of methane into the atmosphere (Egger et al., 2018). Authigenic pyrite (FeS 2 ) and barite (BaSO 4 ) can both form as by-products of SR-AOM, meaning these phases provide a potential archive of methane oxidation (e.g., Borowski et al., 2013;Wood et al., 2021).
Pyrite and barite also provide an important archive for sulfur isotopes in marine environments, which can be used to reconstruct biogeochemical processes that link the sulfur, carbon, and iron cycles (Bottrell and Newton, 2006;Fike et al., 2015). For example, pyrite forms as a by-product of microbial sulfate reduction (MSR) during organoclastic sulfate reduction (OSR) and SR-AOM. There is a large isotopic fractionation associated with MSR, due to the differential reaction rates of the sulfate isotopologues ( 32 S 16 O 4 > 34 S 18 O 4 ; Kaplan and Rittenberg, 1964;Seal et al., 2000;Canfield, 2001a). As a result, pyrite often preserves δ 34 S values that are offset relative to coeval seawater sulfate (δ 34 S pyrite << δ 34 S seawater ; Fike et al., 2015). Stratigraphic variability in δ 34 S values have been used to infer regional to global-scale changes in the sulfur cycle that reflect enhanced pyrite burial (e.g., Goodfellow and Jonasson, 1984) and the size of the marine sulfate reservoir (e.g., Kah et al., 2004).
More recently, studies have shown how sulfur isotope variability may instead be controlled by sedimentary facies and diagenetic processes (Magnall et al., 2016;Pasquier et al., 2017;Marin-Carbonne et al., 2018;Bryant et al., 2019;Richardson et al., 2019;Bryant et al., 2020;Pasquier et al., 2021). For example, the progressive modification of pore fluid sulfate by MSR during diagenesis can result in strong isotopic gradients and a range of δ 34 S pyrite values, depending on the location of pyrite formation in the sediment (Canfield, 2001b;Canfield et al., 2010). In particular, studies have benefited from the generation of isotopic data using microscale techniques (e.g., secondary ion mass spectrometry; SIMS), which enable the determination of paragenetically constrained phase specific δ 34 S values. However, there are relatively few examples where δ 34 S pyrite and proxies for coeval δ 34 S sulfate values have been paired at high spatial resolutions (e.g., Magnall et al., 2016). As a result, the origin of highly positive or 'superheavy' δ 34 S pyrite values (δ 34 S pyrite > δ 34 S SO4 ; Ries et al., 2009) can be particularly difficult to constrain. For example, highly positive δ 34 S values in pyrite in sedimentary units from the Late Devonian period have been linked with MSR in low sulfate water/ sediment columns (Goodfellow and Jonasson, 1984;Sim et al., 2015;Zhang et al., 2020) or formation through hydrothermal or thermochemical sulfate reduction (TSR; Yan et al., 2020).
The Late Devonian was one of the major periods of organic carbon burial in earth history (Klemme and Ulmishek, 1991). In the Selwyn Basin, Canada, clastic-dominated type (CD-type) Zn-Pb ± Ba mineralization is hosted by Late Devonian mudstones (Goodfellow and Jonasson, 1984;Goodfellow, 1987;Hanor, 2000;Johnson et al., 2009;Farquhar et al., 2010;Magnall et al., 2020b). Bedded barite is hosted by unmineralized Late Devonian mudstones in the Selwyn Basin (Fernandes et al., 2017), although the spatial association with the CD-type massive sulfide deposits meant the barite was considered to be a distal expression of sedimentary exhalative (SEDEX) hydrothermal activity (Large et al., 2005;Leach et al., 2005;Goodfellow and Lydon, 2007;Leach et al., 2010). In this SEDEX model, highly enriched δ 34 S values in pyrite and barite have been interpreted to indicate nearly complete sulfate reduction during MSR in a stagnant and stratified anoxic water column (Goodfellow and Jonasson, 1984).
Recent studies, however, have highlighted how barite may have been formed by diagenetic processes before being subsequently replaced during hydrothermal sulfide mineralization (Johnson et al., 2009;Magnall et al., 2016;Johnson et al., 2018;Magnall et al., 2020a;Reynolds et al., 2021). In the diagenetic model, it is proposed that the barite formed in a setting analogous to cold seep environments where methane is oxidized by sulfate through microbial metabolic processes (Greinert et al., 2002;Paytan et al., 2002;Torres et al., 2003;Canet et al., 2014). The sulfate methane transition zone (SMTZ) is a diagenetic redox boundary in organic carbon-bearing sediments that is an important habitat for a consortium of sulfate-reducing and methanotrophic microorganisms (Torres et al., 1996;Machel, 2001;Arning et al., 2015). Organic matter is converted during early diagenesis in such sediments, and the water-soluble products (e.g., acetic acid, CO 2 , and CH 4 ) change the porewater composition (e.g., pH), resulting in a series of hydrogeochemical reactions of dissolution (e.g., feldspar) and precipitation (e.g., barite). Diverse δ 34 S barite values have been recorded at modern cold seeps (up to ∼40‰; Greinert et al., 2002;Wood et al., 2021) that are correspondingly similar to highly positive δ 34 S values in Paleozoic bedded barite (Johnson et al., 2009;Canet et al., 2014). Such highly positive δ 34 S barite values are interpreted to result from the residual sulfate pool at the SMTZ (Canet et al., 2014;Clark et al., 2015). When methane oxidation is coupled to sulfate reduction at the SMTZ and in the presence of an iron source, pyrite is formed at the SMTZ, and this pyrite can have positive δ 34 S pyrite values (e.g., Borowski et al., 2013).
This study integrates high-resolution scanning electron (SEM) microscopy petrography of barite (+ associated barium phases) and pyrite, together with microscale isotopic microanalyses of δ 34 S pyrite , δ 34 S barite, and δ 18 O barite of selected samples from the Late Devonian Canol Formation of the Selwyn Basin. We have targeted samples containing both barite and pyrite to develop paired isotopic constraints on the evolution of sulfur during diagenesis. In particular, we have focused on the precise mechanism by which highly positive δ 34 S pyrite developed in the Canol Formation and discussed the implications for interpreting sulfur isotopes in similar settings. margin of the ancestral North American continent (Mair et al., 2006). Formation of the basin stems from widespread protracted extensional tectonics of the Rodinian supercontinent that led to the re-emergence of the Laurentian craton between 825 and 740 Ma (Martel et al., 2011) with resulting development of the epicontinental margin and Selwyn Basin during the Ediacaran-Cambrian periods (Gordey, 1993). The oldest strata in the basin consist of the syn-rift Neoproterozoic-Terreneuvian Windermere Supergroup overlain by basinal post-rift Paleozoic sedimentary rocks with a combined total thickness of around 7,500 m (Ootes et al., 2013). Abrupt episodic extensions and volcanism during the Early Cambrian, Middle Ordovician, and Devonian are characterized by mafic volcanics within the platform carbonate and the basinal strata (Goodfellow and Lydon, 2007). Collision with an island arc during the Late Devonian is suggested to have led to deformation and subsequent incorporation of the Selwyn Basin strata into the fold and thrust belt of the North American Cordillera (Mair et al., 2006). A sudden change in depositional regime during the Late Devonian led to the deposition of siliciclastic sediments that spread from the margin of Laurentia toward the interior of the craton (Martel et al., 2011). Together with an overlying turbidite unit, the siliciclastic sediments are collectively known as the Earn Group ( Figure 1B; Gordey, 2013;Ootes et al., 2013).
The Earn Group is subdivided into an upper unit of coarsegrained siliciclastic turbidites and sandstones of the Imperial Formation and a lower Canol Formation (Martel et al., 2011). The Canol Formation is of Upper-Devonian (Frasnian-early Famennian) age and consists primarily of dark grey to black mudstones locally calcareous or siliceous with occasional carbonate concretions of variable sizes (Mair et al., 2006;Martel et al., 2011). The formation is extensively widespread in the Mackenzie Mountains region with variable thickness that reaches 400 m (Cecile et al., 1983;Martel et al., 2011;Gadd et al., 2016). The depositional setting of the Canol Formation has been interpreted to be of deep-water marine, formed during the early stages of foredeep basin development (Cecile et al., 1983).

Local Geology
This study builds on an earlier study by Fernandes et al. (2017) located northeast of the Macmillan Pass district ( Figure 1A). The sedimentary rocks of the Canol Formation in the study area consist of gently open-folded, bedded coarser-grained siliciclastics and organic-rich mudstones that are moderate to steeply dipping and weathering to silver color (Fernandes, 2011;Fernandes et al., 2017). These lithologic units host several stratiform barite beds that form barite horizons in equivalent Devonian stratigraphic intervals on a regional scale (e.g., Goodfellow and Jonasson, 1984;Smith et al., 1993;Magnall et al., 2016;Fernandes et al., 2017).  Goodfellow & Lydon, 2007). Major barite and clastic-dominated (CD-type) Zn-Pb deposits are indicated by the yellow cycle and red triangle, respectively. The red area in B depicts the current study location. (B) Devonian stratigraphic sequence of the Selwyn Basin from south to north (modified from Martel et al., 2011;Ootes et al., 2013;Morrow, 2018 In the Northwest Territories' part of the Selwyn Basin, 22-72 m thick barite sequences occur within stratigraphic sections at Bunk-2, Bunk-2, Cowan, NAFCAC-1, NAFCAC-2, Anita, Axe, Harp, and Wise locations (Figures 2, 3;Fernandes et al., 2017). These baritebearing sections occur as topographic highs, in a broad 200 km long NW-SE trend, are all confined to the upper parts of the Canol Formation, below the unconformity with the fine-grained siliciclastics of the Imperial Formation and are considered to be Frasnian in age (Martel et al., 2011;Fernandes et al., 2017). Mineralogically, the mudstones consist of varying amounts of clay minerals, quartz, and organic matter, with barite, pyrite, and Ba-bearing feldspar (hyalophane (K, Ba) [Al (Si, Al) 3 O 8 ] and cymrite BaAl 2 Si 2 (O, OH) 8 · H 2 O) constituting the major accessory minerals (Fernandes, 2011). Fernandes et al. (2017) further describe the intricate relationship between barite and Ba-bearing carbonates and silicates; witherite (BaCO 3 ) forms a textural association with barite, replacing both laminated and nodular barite grains, while hyalophane and cymrite are 0.5-2 mm crystals found replacing barite and witherite in both the laminae and nodules (Fernandes et al., 2017).

Bedded Barite Mineralization Style in the Canol Formation
There are two types of barite (laminated and nodular) identified in the upper Canol Formation that occur interlaminated together; nodular barite was dominant over laminated form (Figure 3; Fernandes et al., 2017). Brief descriptions of the barite forms from Fernandes et al. (2017) are highlighted below and shown in Figure 4.
The barite nodules in mudstones are spherical, ellipsoidal, or irregular ( Figures 4A-C). The nodules often contain barite crystals that are rosette or tabular with a size range from less than 100 μm to 1.2 cm ( Figure 4E). Witherite, hyalophane, cymrite, and quartz are associated primarily with the barite crystals ( Figures 4G,H). The laminated barite has been described to primarily occur as clusters of intergrown anhedral crystals within 50-100 µm laminae (e.g., Figure 4D). Associated minerals include pyrite, cymrite, quartz, and hyalophane with intercalated clay-rich laminae ( Figure 4H). Grey mudstone with about a meter of laminated and nodular barite occurs in the Bunk-2 section between dark grey siltstone hosting nodular barite ( Figure 4I). The barite crystals range from <10 to 70 μm, either concordant with the lamina or irregular with no specific direction.

METHODOLOGY Sampling
Sixty-five (65) mudstone samples from the nine barite-bearing sections from the suite of samples investigated by Fernandes (2011) were resampled. Samples containing pyrite and barite were targeted, and their mineralogical and paragenetic relationships were examined using binocular microscopy as a first step.

Petrography
Detailed petrographic examination of the mineralogical and textural relationships was carried out on thin sections using an Olympus BX51 microscope, equipped with a Sc50 camera, in transmitted and reflected light settings. Polished thin sections were additionally prepared and carbon-coated to a thickness of 20 nm for further examination and imaging using an electron Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 784824 probe microanalyzer (EPMA) and SEM. Backscatter electron images (BSE) were obtained using Japan Electron Optics Limited (JEOL) JXA-8530F Hyperprobe, equipped with a wavelength and energy dispersive spectrometry combined system. The EPMA was operated using a beam diameter between 1 and 3 μm, beam current of 15 nA, and accelerating potential of 15 kV, in secondary electron (SE) and BSE modes. Organic petrography using a standard reflection microscope showed that randomly distributed organoclasts occur as < 2 µm sized particles. Due to the thermal overmaturity, the particles are inertinite and prevent a broader reconstruction of organic matter type.

Secondary Ion Mass Spectrometry (SIMS) Analyses of Sulfur and Oxygen Isotopes
Microdrills of regions of interest (n 54) were made on polished sections to obtain suitable subsamples, using a 4 mm diameter diamond core drill, from the whole sample suite in Fernandes (2011). Several representative subsamples were cast into 25 mm epoxy pucks, together with reference materials (RMs) of pyrite S0302A (δ 34 S V-CDT 0.0 ± 0.2‰; Liseroudi et al., 2021) and barite S0327 (δ 34 S V-CDT 11.0 ± 0.5‰; δ 18 O V-SMOW 21.3 ± 0.2‰; Magnall et al., 2016). Re-examination and further BSE imaging were carried out with EPMA after carbon-coating, with subsequent 30 nm gold coating added to the mounts before isotope measurement. Microscale isotopic analyses were carried out using Cameca IMS1280 large-geometry secondary ion mass spectrometer (SIMS) operated in multi-collector mode at the NordSIMS laboratory, Stockholm, Sweden. For the measurements, a 133 Cs + , 20 kV impact energy primary beam was utilized. The beam current was 1 nA for all barite analyses, and the larger pyrite targets, yielding a ca. 10 μm spot; a 400-500 pA beam was used for smaller pyrite targets yielding a ca. 6 μm spot. A normal incidence low-energy electron flooding gun was utilized for charge compensation. Sulfur and oxygen isotopes were determined in separate analytical sessions in which secondary ion signals of 32 S and 34 S or 16 O and 18 O were measured simultaneously in two Faraday cups. Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 784824 A total of 1,032 sulfur and oxygen isotope measurements (δ 34 S pyrite 200; δ 34 S barite 485, δ 18 O barite 338) on pyrite and barite grains were obtained using automated analytical sequences in which every 6 to 7 measurements are followed by 1-2 RM analyses. The δ 34 S barite and δ 18 O barite values measurements were carried out on the same barite grains to capture the covariation between the sulfur and oxygen isotopic systems. Within-session drift and instrumental mass fractionation (IMF) were corrected using the regularly interspersed analyses of the RMs in each session. The IMF-corrected 34 S/ 32 S ratios are reported relative to Vienna Canyon Diablo Troilite (V-CDT) and  Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 784824 7 O/ 16 O ratios relative to Vienna Standard Mean Ocean Water (V-SMOW), using conventional delta notation: External analytical reproducibility (1 σ) was typically ±0.04‰ δ 34 S for pyrite, ± 0.15‰ δ 34 S, and ±0.12‰ δ 18 O for barite. The reproducibility for each session is the standard deviation of the reference material measurements in that session. For any individual measurement, the external uncertainty is propagated together with the within-run uncertainty for an overall value. Post SIMS SEM imaging was carried out on each barite and pyrite spot for confirmation of target integrity; measurements on pyrite-barite grain boundaries were discarded (n 14).

Mineralogy and Paragenesis
Different pyrite and barite formation stages have been defined in terms of shape, size, and distribution, broadly comparable across the different stratigraphic intervals ( Figure 5). The different stages of pyrite and barite are described in terms of their relative timing of formation (paragenesis). Detailed descriptions are given below.
Stage 1: The earliest stage of pyrite (Py-1) comprises framboids present in all the sections (apart from NAFCAC-2 and Anita). Framboidal pyrite is particularly enriched in the mudstones of Harp and Axe sections. Py-1 is commonly located in the inter-and intra-granular pore space of the mudstone laminae and nodules ( Figures 6A,B). The framboids have a broad size distribution, mostly <45 µm ( Figure 6C), although spherical to irregular clusters of polyframboids reach up to 80 µm in diameter ( Figure 6D). The individual framboid microcrystals are <5 μm, mostly equant, cubic, or pyritohedral ( Figure 6C). The first stage of barite (Brt-1) is microcrystalline, <35 µm with subhedral to euhedral disseminated grains, and is found in mudstones in all sections. The Brt-1 is intergrown with quartz and cymrite in clay-rich or quartz-dominated matrices ( Figures 7A-D). The Brt-1 and Py-1 grains rarely occur together, and establishing their paragenetic relationship was consequently difficult.

DISCUSSION
The mineralogical paragenesis and microscale isotopic constraints for pyrite and barite enable the reconstruction of the sulfur cycle during the deposition of the Canol Formation in the Late Devonian. The paired isotope data (δ 34 S pyrite , δ 34 S barite , and δ 18 O barite ) can be used to unravel the fate and behavior of the archived sulfate and to interpret the end-member δ 34 S pyrite values. In particular, the highly positive values occur in samples from several regionally correlative sections; this suggests the formation of these distinctive isotopic values could represent a regionally important (> 10 s km) process during the Late Devonian Selwyn Basin.

Barite Formation
The formation of pelagic barite in marine environments is initially associated with sinking particulate organic matter (Paytan et al., 2002;Gonzalez-Muñoz et al., 2012;Martínez Ruíz et al., 2020), which may then be recycled under reducing conditions during diagenesis (Torres et al., 1996). Under opensystem conditions, biogenic barite preserves δ 34 S and δ 18 O values that represent unmodified seawater sulfate (Paytan and Griffith, 2007;Griffith and Paytan, 2012 (John et al., 2010;Chen et al., 2013). However, the barite from all stratigraphic sections in this study preserves higher δ 34 S barite values (median 45‰), representing a substantial offset from Late Devonian seawater sulfate (Figure 11). Barite crystals with variable isotopic ratios are not restricted to a single section, depth, or lithology; the analyses reveal isotopic heterogeneity within and between grains of the same barite generation that are only a few microns apart ( Figures 13B-D). This lack of isotopic distinction between the different barite types and the distribution of the isotopically heterogenous barite within the lithologies and sections may indicate a similar environment of formation across all stratigraphic sections (e.g., Magnall et al., 2016; Figure 12B).
The size (>5 µm) and morphology of the barite crystals are also consistent with barite precipitation in diagenetic pore fluids (Paytan et al., 2002;Paytan and Griffith, 2007). The two barite generations (Brt-1 and Brt-2) are isotopically indistinct ( Figure 12A), despite paragenetic relationships that suggest barite precipitation during two stages (Figures 6,9). Both generations of barite are also surrounded by the siliciclastic and organic-rich host mudstones suggesting formation in early diagenesis near the seafloor (Aplin and Macquaker, 2011). The formation of the barite nodules (Brt-2b) might have coincided with increased compaction of the sediments, with clay minerals and organic matter dehydration possibly causing the nodules to form into ellipsoidal and irregular shapes (e.g., Figure 9C), parallel to the bedding plane (Goldberg et al., 2006;Paytan and Griffith, 2007;Zan et al., 2020). Diagenetic barite formation is generally associated with the sulfate methane transition zone (SMTZ), where opposing diffusional fluxes of methane and sulfate interact (Barnes and Goldberg, 1976;Reeburgh, 1976;Jørgensen and Kasten, 2006). Barium is soluble in the strongly reducing CH 4 -bearing fluids but will form diagenetic barite upon mixing with downward diffusing sulfate (Torres et al., 1996;Dickens, 2001). The depth of diagenetic barite formation and the degree to which pore fluid sulfate has been modified via MSR will significantly influence δ 34 S barite and δ 18 O barite values (Bottrell and Newton, 2006;Goldberg et al., 2006;Antler et al., 2013). Compared to the mineral separate analyses from the previous study (Fernandes et al., 2017), the microscale δ 34 S barite and δ 18 O barite values in this study are biased towards an end member that represents strongly modified seawater ( Figure 12B). Samples containing both pyrite and barite were targeted in this study, meaning the barite only samples that contained lower δ 34 S barite and δ 18 O barite values are underrepresented in the microscale dataset. Nevertheless, when interpreted together, the bulk rock and microscale data plot along a consistent trend ( Figure 12B) that is typical of the progressive modification of seawater sulfate via MSR (e.g., Antler et al., 2013;Pellerin et al., 2019).
The overall trend between unmodified Late Devonian seawater and higher δ 34 S barite values indicates that barite formed under progressively sulfate limited conditions in which the rate of sulfate depletion exceeds that of sulfate resupply (Torres et al., 1996;Fike et al., 2015). The δ 34 S barite values of Brt-2b represent the most evolved isotope signatures ( Figure 12A), consistent with precipitation under sulfate limitation during later stages of diagenesis (Turchyn and Schrag, 2006;Aller et al., 2010;Gomes and Johnston, 2017). Later changes in the barite front may have resulted in the formation of the vein and pseudomorphic barite (Brt-2c), which formed within pore spaces of existing pyrite and barite (Figures 7, 10) but appears to have precipitated from a similar sulfate pool (Figure 11) during the second stage of paragenesis.
Covariation between δ 18 O and δ 34 S values can provide further information on the diagenetic environment of barite formation. In sediment pore fluids, the slope of the apparent linear phase (SALP; Antler et al., 2013), which describes covariation between δ 18 O and δ 34 S values, has been linked to sulfate reduction rate (SRR; Böttcher et al., 1998;Böttcher et al., 1999;Aharon and Fu, 2000;. At low SRR, a high degree of reversibility in the enzymatic pathway of sulfate reduction is thought to promote oxygen isotope exchange between intermediate sulfur phases (e.g., sulfite) and H 2 O (Fritz et al., 1989;, resulting in high SALP values. In contrast, where δ 34 S values increase without a corresponding increase in δ 18 O values (low SALP), it is thought that high SRR results in a lower degree of reversibility (Antler and Pellerin, 2018). At high SRR, kinetic isotope effects will provide the primary control on δ 18 O values, resulting in a SALP value of ∼0.25 that represents an oxygen fractionation factor that is approximately 25% of sulfur (Mizutani and Rafter, 1973). Importantly, diagenetic mineral phases that contain sulfate (e.g., barite, celestine and carbonates) have the potential to preserve the SALP signature (Antler and Pellerin, 2018). There is a strong positive correlation between δ 18 O and δ 34 S values in barite from the Canol Formation ( Figure 12B), which corresponds with a low SALP value and high SRR. Similar low SALP values have been described in authigenic carbonate formed in ancient cold seep environments (e.g., Feng et al., 2016).

Barite Replacement
Barium concentrations in basinal waters are controlled by sulfate concentration and reduction (Torres et al., 2003). The low solubility of barite means the stability field of BaSO 4 overlaps with conditions in which reduced sulfur is the dominant sulfur species, rather than sulfate (Figure 14). The dissolution and FIGURE 11 | Box and whisker plot of barite and pyrite δ 34 S value from this study. Grey area indicates the extent of Late Devonian seawater (John et al., 2010). The colors are intended as visual aids, highlighting the barite and pyrite types.
Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 784824 replacement of diagenetic barite may then proceed under conditions of extreme sulfate limitation that typically develop in carbonaceous sediments (Hanor, 2000). The diagenetic barite in the Canol Formation has been replaced by Ba-carbonates (witherite and barytocalcite), followed by Ba-feldspars (hyalophane and cymrite; Figures 9, 13D).  (Fernandes et al., 2017) and SIMS analyses of barite from Macmillan Pass (Magnall et al., 2016). The regression analysis includes the SIMS and mineral separate data for Canol Formation barite but excludes the outliers from either dataset (highlighted by grey circles). The blue box represents the range of constraints for Late Devonian seawater (John et al., 2010;Chen et al., 2013). The presence of replacive Ba-bearing phases may imply the progressive diagenetic replacement of barite arising from continuous depletion of the sulfate concentration (e.g., Hanor, 2000). High Ba concentrations (around six orders of magnitude relative to seawater) coupled with extreme sulfate depletion have been shown to result in witherite formation (Maynard and Okita (1991). Barium concentrations in modern seawater are between ∼5 and 20 μg/L in the open ocean (Hanor, 2000) and up to 60 μg/ L occurs in anoxic deep marine environments (Falkner et al., 1993). Thus, Ba concentrations may increase with depth (Paytan and Griffith, 2007;Carter et al., 2020), mediated by microbial activities in the presence of organic matter. Two conditions are suggested by Maynard and Okita (1991) as a prerequisite for witherite replacement of barite; a) a closed or restricted system with sulfate resupply << sulfate reduction; b) high organic matter contents for diagenetic barite conversion to witherite. However, barite dissolution and witherite formation have also been shown as a function of temperature, CO 2 dissolution (Busenberg and Plummer, 1986), and pH (Melero-García et al., 2009;Hill et al., 2014), with dissolved carbonate, likely supplied from organic matter degradation (Maynard and Okita, 1991;Hanor, 2000).

Pyrite Formation
Pyrite in the Canol Formation samples preserves end-member δ 34 S values ( Figure 11) that are associated with two morphologically distinct stages of pyrite formation ( Figures 13A,B). Framboidal pyrite (Py-1) in the Canol Formation preserves δ 34 S values that are significantly lower than Late Devonian seawater (Figure 11), which represents a large isotopic fractionation (≤75‰). The δ 34 S py-2 values are more positive (mean +31.0‰) than upper constraints for Late Devonian seawater (+25‰; John et al., 2010;Chen et al., 2013). The formation of highly positive δ 34 S pyrite values has been attributed to various processes that include MSR (Drake et al., 2018), sulfide reoxidation (Kah et al., 2016), TSR (Cui et al., 2018;Yan et al., 2020), rapid sedimentation (Pasquier et al., 2017), and sulfate limitation (Ries et al., 2009). Highly positive δ 34 S values from mineral separate analyses of pyrite and barite in siliciclastic units from the Selwyn Basin have been interpreted to have developed from water mass restriction and euxinic conditions in the Late Devonian (Goodfellow and Jonasson, 1984). However, the pyrite paragenesis indicates Py-2 precipitated below the SWI during diagenesis (e.g., Figure 8) when sulfate limitation would have developed on a much smaller pore fluid scale. Importantly, there is no observable alteration relating to hydrothermal fluids (e.g., Cui et al., 2018;Yan et al., 2020), which rules out any high-temperature origin for Py-2. Highly positive δ 34 S pyrite values could have developed due to high sedimentation rates reducing the diffusional exchange between diagenetic pore fluids and overlying seawater (e.g., Pasquier et al., 2017). However, there does not appear to be any lithological control on the development of highly positive δ 34 S pyrite values in the Canol Formation, which may imply that depositional environment was not the primary control on δ 34 S pyrite values.
The combination of highly positive δ 34 S pyrite values and evidence of barite solubility represents a trend of progressive sulfate depletion, which is typical of the diagenetic evolution of anoxic sediments (e.g., Torres et al., 1996;Aloisi et al., 2004). Importantly, there is evidence of textural equilibrium between Py-2 and Brt-2 (e.g., Figure 8) and it is apparent that δ 34 S pyrite values do not exceed coeval δ 34 S barite values (Figure 11), meaning they are not strictly "superheavy" (δ 34 S pyrite > δ 34 S SO4 ; Ries et al., 2009;Cui et al., 2018). A similar pyrite paragenesis and assemblage between Py-2 and Brt-2 has also been described at Macmillan Pass, where it was linked with two separate stages of diagenetic pyrite formation associated with OSR and AOM-SR (Magnall et al., 2016). Compared to the offset between Py-1 and Late Devonian seawater (Δ 34 S 43.3‰), the Δ 34 S for Brt-2 and Py-2 is smaller ( Figure 11). This reduced Δ 34 S value could represent a smaller isotopic fractionation (ε 34 S), possibly linked with higher SRR that are typical of methane and gas seeps (Deusner et al., 2014). Indeed, the low SALP value that is recorded in barite ( Figure 12A) would support this interpretation, although it is also possible that Py-2 and Brt-2 did not precipitate precisely at the same time from the same porewater fluids. Nevertheless, the overall evidence of high SRR and increasing sulfate depletion still provides the most plausible explanation for the development of highly positive δ 34 S pyrite values in the Canol Formation.

Implications
Previous studies on diagenetic pyrite and barite linked highly positive δ 34 S values to sulfate limitation in a euxinic water column (Goodfellow and Jonasson, 1984). More recent studies have now shown that these isotopic values developed as part of a diagenetic assemblage associated with the SMTZ, where sulfate limitation occurred at the pore fluid scale (e.g., Magnall et al., 2016;Johnson et al., 2018). Notably, the same 2-stage paragenesis comprising framboidal pyrite followed by an assemblage of euhedral pyrite and barite has been identified in Late Devonian strata that host CD-type deposits in the Macmillan Pass district (Magnall et al., 2016). Thus, the formation of highly positive δ 34 S pyrite values and bedded barite in the Canol Formation is evidence of a similar diagenetic assemblage that has been preserved in unmineralized Late Devonian stratigraphy. Moreover, the preservation of this  Maynard and Okita (1991).
Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 784824 assemblage in multiple stratigraphic sections in the Canol Formation implies periodical stability in the SMTZ on a large scale (> 10 km) in order to allow pyrite and barite to accumulate (e.g., Lin et al., 2017;Liu et al., 2019). In the modern oceans, shallow SMTZ and high methane fluxes are located along productive continental margins with high sedimentation rates (Egger et al., 2018). Other factors controlling the depth of the SMTZ include organic matter content and decomposition and sulfate concentrations (Torres et al., 2003;Jørgensen et al., 2004;Borowski et al., 2013;Lin et al., 2016;Lin et al., 2017;Liu et al., 2019;Liu et al., 2020;Liu et al., 2021). In the Canol Formation, high total organic carbon (1.5-9.2 wt%; Kabanov and Gouwy, 2017) is associated with facies characterized by dynamic sedimentation, in which bioturbation has been linked with partial oxygenation (Biddle et al., 2021). The stabilization of the SMTZ is likely to have required constant fluxes of methane and sulfate coupled with a hiatus in sedimentation, and changes to any of these parameters could potentially have resulted in the alteration of authigenic mineral phases (Arning et al., 2015). The relatively minor replacement of Py-2 by veined Brt-2c ( Figure 10) may therefore suggest abrupt albeit short-lived fluctuation of the SMTZ due to changes in either methane or sulfate fluxes within the sediments. Barite solubility is drastically increased below the SMTZ, and where sulfate concentrations diminish, Ba diffuses upward and reacts with other available cations (Maynard and Okita, 1991;Hanor, 2000). Considering that SR-AOM in the SMTZ also leads to the generation of alkalinity, this may explain the origin of barium carbonate in these samples (e.g., Figure 13D). The origin of carbonate in these samples could be evaluated using carbon isotopes, but this is beyond the scope of this study.
Similar barite-pyrite-Ba-carbonate/silicate assemblages have also been reported in other Paleozoic strata (e.g., Maynard and Okita, 1991;Jewell, 2000;Koski and Hein, 2004), including Qinling-Daba region, southern China (Wang and Li, 1991;Xu et al., 2016). The low seawater concentrations in the Lower Paleozoic (ca. <1-10 mM; Horita et al., 2002;Brennan et al., 2004;Gill et al., 2007), relative to modern ocean seawater sulfate that is about 28 mM (Rickard, 2012;Fike et al., 2015; Jorgensen Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 784824 et al., 2019), would have restricted the amount of sulfate resupply below the seafloor in organic-rich sediments (Greinert et al., 2002;Paytan et al., 2002). Continuous fractionation of the porewater sulfate through OSR and SR-AOM likely led to the formation of highly positive δ 34 S pyrite values in this study and other similar settings. We would suggest that highly positive δ 34 S pyrite values and barite formation, dissolution, and replacement by other Babearing phases might have been a more general feature of methane diagenesis in the lower Paleozoic.

CONCLUSION
The diagenetic sulfur cycle in the Selwyn Basin has been reconstructed using microscale (SIMS) paired isotope constraints on pyrite and barite in samples from multiple stratigraphic sections of Late Devonian sedimentary rocks. Two distinct stages of pyrite formation formed via microbial sulfate reduction (MSR) under contrasting levels of sulfate availability. The initial stage of pyrite formation developed during early diagenesis under relatively opensystem conditions, which resulted in the precipitation of the framboidal pyrite (Py-1) and preservation of negative δ 34 S values. Deeper in the sediment profile, MSR resulted in progressive sulfate depletion and the development of the sulfate methane transition zone (SMTZ), where sulfate reduction was coupled with the anaerobic oxidation of methane (SR-AOM).
Highly positive δ 34 S values in pyrite developed as a result of high sulfate reduction rates and progressive depletion of sulfate. Sulfate limited, methane-rich diagenetic fluids beneath the SMTZ provided conditions under which barium was soluble. Barite formed when opposing diffusional fluxes of barium and sulfate bearing fluids mixed at the SMTZ. Importantly, the paragenetically constrained analyses indicate that highly positive δ 34 S pyrite values formed during diagenesis in multiple correlated stratigraphic sections from the Late Devonian in the Selwyn Basin, indicating this diagenetic assemblage was a regional feature. We propose that the formation of highly positive δ 34 S values in pyrite, dissolution of barite, and replacement by other Ba-bearing phases, could have been a more general features of methane diagenesis in the lower Paleozoic.

DATA AVAILABILITY STATEMENT
The dataset presented in this study can be found online at the data repository of the GFZ German Research Centre for Geosciences, Potsdam, Germany. GFZ Data Services. https://doi.org/10.5880/ GFZ.3.1.2021.006.

AUTHOR CONTRIBUTIONS
SG and JM designed the study. HG prepared the samples and performed the petrographic analysis (optical microscopy, EPMA-EDS, and SEM). MW conducted the SIMS analyses in coordination with HG and JM. HG interpreted the data and wrote the manuscript with contributions from JM, SG, H-MS, and MW.