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).

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 barite-bearing 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)₃O₈] and cymrite BaAl₂Si₂(O, OH)₈ ∙ H₂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₃) 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).

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.

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.

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 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.

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 (δ³⁴S _V-CDT = 0.0 ± 0.2‰; Liseroudi et al., 2021) and barite S0327 (δ³⁴S_V-CDT = 11.0 ± 0.5‰; δ¹⁸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 ¹³³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 ³²S and ³⁴S or ¹⁶O and ¹⁸O were measured simultaneously in two Faraday cups.

A total of 1,032 sulfur and oxygen isotope measurements (δ³⁴S_pyrite = 200; δ³⁴S_barite = 485, δ¹⁸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 δ³⁴S_barite and δ¹⁸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 ³⁴S/³²S ratios are reported relative to Vienna Canyon Diablo Troilite (V-CDT) and ¹⁸O/¹⁶O ratios relative to Vienna Standard Mean Ocean Water (V-SMOW), using conventional delta notation: δ 34 S  ( ‰ ,  V − CDT ) = [ ( (   34 S /   32 S )   sample (   34 S /   32 S )  V − CDT ) − 1 ] × 1000 δ 18 O  ( ‰ ,  V − SMOW ) = [ ( (   18 O /   16 O )  sample (   18 O /   16 O )  V − SMOW ) − 1 ] × 1000

External analytical reproducibility (1 σ) was typically ±0.04‰ δ³⁴S for pyrite, ± 0.15‰ δ³⁴S, and ±0.12‰ δ¹⁸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).

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.

Stage 2: The second stage, pyrite (Py-2), forms porous, subhedral to idiomorphic grains, primarily concentrated in laminae and nodules of the mudstones (Figures 8A,B). In the lower part of the Harp section, where the mudstones contain no barite laminae or nodules, Py-2 either occurs as disseminated crystals (<70 µm) or forms aggregates that are up to 2.7 mm (Figure 4J). The Py-2 grains sometimes form overgrowths on Py-1 and may contain inclusions of framboids, quartz, and barite (Figures 8C,D). Barite (Brt-2) forms subhedral to anhedral stratiform barite crystals (Brt-2a, Figure 9A) and spherical to irregular barite nodules (Brt-2b, Figure 9B). Together with other barium phases, this barite replaces the earlier formed Brt-1 (Figures 6D, 7B). The Brt-2a mainly occurs as well-formed horizontal (Figure 4D) and wavy to irregular laminae with intercalations of clay-rich laminae (Figure 9A), nodular Brt-2b is the dominant form of barite in the samples (Figure 9C). The Brt-2a laminae tend to wrap around rosette and nodules of Brt-2b when they occur together, associated with witherite, quartz, hyalophane, and cymrite (Figures 9D–F). Importantly, Brt-2b is intergrown with Py-2 with clear planar crystal boundaries, suggesting coprecipitation (e.g., Figures 8B,D). Notably, Brt-2 and Brt-1 are replaced by witherite, mainly in the dark grey mudstones of Harp, Cowan, and Axe sections. The vein and pore-filling barite (Brt-2c) comprises precipitation in the pore, intergranular spaces, and fractures and is most commonly observed in the Axe, NAFCAC-1, and Harp sections. Barite-2c tends to form replacement and feeder (conduits) textures, including veins, veinlets, lensoids, and idiomorphic barite crystals overgrowing or replacing earlier formed barite and pyrite and cymrite grains (Figures 7B, 10).

Stage 3: The third stage in the paragenesis comprises Ba-bearing minerals, including cymrite and hyalophane, precipitating as pseudomorphs of earlier formed mineral phases (Figure 7C). In essence, the formation of cymrite and hyalophane, together with witherite and quartz, continues after the precipitation of Py-2 and Brt-2 (Figure 10C).

The results of the SIMS analyses of the isotopic compositions of pyrite and barite are provided in Grema et al. (2021) and presented in Figures 11, 12. Pyrite has a large range of δ³⁴S values (Figures 11, 13). Framboidal pyrite (Py-1) preserves mostly negative δ³⁴S values (Figure 13A), between −29.4‰ and −9.4‰, with a mean δ³⁴S value of −21.0 ± 6.3‰ (1σ, n = 23). Pyrite-2 preserves more positive δ³⁴S values (Figure 13B) between +9.4‰ and +44.5‰ and the mean is +31.0 ± 7.7‰ (1σ, n = 177).

The δ³⁴S_barite values in this study overlap (Figure 12A) and range between + 37.1‰ and + 67.9‰ (+51.1 ± 7.7‰, 1σ) with an outlier of + 15.3‰ from the mudstones of the NAFCAC-2 section. Measured δ¹⁸O_barite values from the sections also overlap and are between + 8.8‰ and + 23.9‰ (+21.0 ± 1.7‰, 1σ). The Brt-2a δ³⁴S values are between + 37.1‰ and + 60.1‰ with a mean δ³⁴S value of +49.2 ± 4.7‰ (1σ, n = 63). Corresponding δ¹⁸O values of +19.5‰ to + 23.7‰ with mean δ³⁴S value of +21.8 ± 1.5‰ (1σ, n = 24) are preserved in the Brt-2a barite (Figure 12A). The δ³⁴S and δ¹⁸O values of Brt-2b (+ 51.4 ± 8.1‰, 1σ) are between + 38.8‰ and +67.9‰ (n = 391) and +8.8‰ to +23.9‰ (+21.0 ± 1.7‰, 1σ, n = 284) respectively, with an outlier of δ¹⁸O value of −5.5‰ (corresponding δ³⁴S value of + 40.1‰). The Brt-2b barite from the Cowan and NAFCAC-2 sections record the highest δ³⁴S values between + 60.8‰ and + 67.9‰, with corresponding δ¹⁸O values between +18.4‰ and +23.5‰ (Figure 12A). Analyses of barite veins produced δ³⁴S values between + 44.4 and + 51.1‰ (+48.5 ± 3.6‰, 1σ) and δ¹⁸O values between + 16.6‰ and + 21.5‰ (+19.5 ± 2.6‰, 1σ) respectively. Overall, there is observed covariation (Figures 13C,D) between δ³⁴S_barite and δ¹⁸O_barite values, as indicated by the coefficient of determination (r² = 0.52; Figure 12A).

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 open-system conditions, biogenic barite preserves δ³⁴S and δ¹⁸O values that represent unmodified seawater sulfate (Paytan and Griffith, 2007; Griffith and Paytan, 2012). The constraints for unmodified Late Devonian seawater sulfate (δ³⁴S values = 20 and 25‰; δ¹⁸O values = 10–16‰) are provided by analyses of carbonate associated sulfate (CAS) from the Frasnian-Famennian boundary (John et al., 2010; Chen et al., 2013). However, the barite from all stratigraphic sections in this study preserves higher δ³⁴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₄-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 δ³⁴S_barite and δ¹⁸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 δ³⁴S_barite and δ¹⁸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 δ³⁴S_barite and δ¹⁸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 δ³⁴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 δ³⁴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 δ¹⁸O and δ³⁴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 δ¹⁸O and δ³⁴S values, has been linked to sulfate reduction rate (SRR; Böttcher et al., 1998; Böttcher et al., 1999; Aharon and Fu, 2000; Brunner et al., 2005). 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₂O (Fritz et al., 1989; Brunner and Bernasconi, 2005), resulting in high SALP values. In contrast, where δ³⁴S values increase without a corresponding increase in δ¹⁸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 δ¹⁸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 δ¹⁸O and δ³⁴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).

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₄ overlaps with conditions in which reduced sulfur is the dominant sulfur species, rather than sulfate (Figure 14). The dissolution and 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).

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₂ 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 in the Canol Formation samples preserves end-member δ³⁴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 δ³⁴S values that are significantly lower than Late Devonian seawater (Figure 11), which represents a large isotopic fractionation (≤75‰). The δ³⁴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 δ³⁴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 δ³⁴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 δ³⁴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 δ³⁴S_pyrite values in the Canol Formation, which may imply that depositional environment was not the primary control on δ³⁴S_pyrite values.

The combination of highly positive δ³⁴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 δ³⁴S_pyrite values do not exceed coeval δ³⁴S_barite values (Figure 11), meaning they are not strictly “superheavy” (δ³⁴S_pyrite > δ³⁴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 (Δ³⁴S = 43.3‰), the Δ³⁴S for Brt-2 and Py-2 is smaller (Figure 11). This reduced Δ³⁴S value could represent a smaller isotopic fractionation (ε³⁴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 δ³⁴S_pyrite values in the Canol Formation.

Previous studies on diagenetic pyrite and barite linked highly positive δ³⁴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 δ³⁴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 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 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 δ³⁴S_pyrite values in this study and other similar settings. We would suggest that highly positive δ³⁴S_pyrite values and barite formation, dissolution, and replacement by other Ba-bearing phases might have been a more general feature of methane diagenesis in the lower Paleozoic.