Fresh specimens for geochemical analysis were collected from sediment cores at the Alberta Energy Regulator Core Research Centre (Calgary, Alberta, Canada). Cores were described to identify sedimentary facies and formation unit boundaries, which were differentiated based on ichnology, bioturbation intensity, grain size, fossil prevalence, reactivity to 10% HCl, fracture frequency, sedimentary structures, and color (Frucci, 2021). The cores were then depth-corrected to match log depth by comparing core-observed rock attributes and gamma-ray, induction, and neutron/density log response (Visy, 2022). Additional samples from two outcrop exposures of the Exshaw Formation at Jura Creek and Crowsnest Pass (Locations 21 and 22, Figure 1) were also analyzed in this study. For each outcrop sample, exposed/weathered surfaces were removed using a grinder before being broken into small chips using a rock hammer and powdered in a shatterbox (SPEX Industries Inc., Catalog No. 8500). For C and N isotopic analysis, the powdered samples were decarbonated using 10% HCl (v/v) for 48 h and centrifuged at 3000 RCF for 4 min. The supernatant was then decanted, and acid was added to test for effervescence that indicates residual inorganic carbon. After decarbonation, samples were rinsed three times or until pH neutral with deionized water, centrifuged, and the supernatant discarded. The solid residue was dried using a freeze drier for 24–48 h.

For determination of bulk δ¹³C_org and δ¹⁵N_bulk, 10–30 mg decarbonated samples were weighed into tin capsules and loaded into a zero-blank auto-sampler on an elemental analyzer (EA; Costech, ECS 4010) coupled with an isotope ratio mass spectrometer (IRMS; Thermo-Electron Delta V Advantage) in the Baylor Geosciences Stable Isotope Lab. We analyzed the N isotopic composition of a set of 8 non-acidified Exshaw Shale samples to examine the effect of acidification on δ¹⁵N values. On average the acidified samples were ¹⁵N-enriched by 0.68 ± 0.69‰. While this difference is not trivial, the geographic and stratigraphic trends are captured by the δ¹⁵N values of acidified samples. The N₂ and CO₂ produced in the EA passed through a thermal conductivity detector for total nitrogen and organic carbon determination prior to transfer to the IRMS. Samples were analyzed in batches of 10 bracketed by blanks and an internal standard acetanilide calibrated to international isotope standards (USGS 40: δ¹³C = −26.39‰, δ¹⁵N = −4.52‰; USGS 41: δ¹³C = +37.63‰, δ¹⁵N = +47.57‰). Isotopic composition is reported in conventional δ notation relative to Vienna Pee Dee Belemnite (VPDB) for carbon and atmospheric N₂ for nitrogen. Standard deviations range between 0.01 and 0.1‰ for δ¹³C_org and 0.12–0.24‰ _for δ¹⁵N.

Trace and major element concentrations were measured using a portable Bruker TRACER 5i X-Ray Fluorescence (XRF) spectrometer. To assess redox conditions, vanadium, molybdenum, uranium, nickel, chromium, and cobalt concentrations are reported and normalized to aluminum to account for terrigenous inputs of the trace elements (TE). Four standards, Zentrales Geolosches Institut Black Shale (ZGI TS), Mintek Carbonaceous Shale (SARM 41), Geological Survey of Japan Black forest soil (JSO-1), and Geological Survey of Japan Porites sp. Coral (Jcp-1), were analyzed before and after each batch of samples for a total of 10 replicates. These values were used to determine accuracy and precision of analyses and calculate response factors to assess the stability of the internal “Mudrock” calibration. For ZGI TS, the relative uncertainty (1σ standard deviation) calculated during analysis was 10.5% for Al, 3.8% for V, 3% for Mo, 18.7% for U, 5.8% for Ni, 13% for Co, and 5.3% for Cr. Trace element enrichment factors (EF) were calculated relative to the reference elemental concentration values reported for the global average Post-Archean Australian Shale (PAAS) (Taylor and McLennan, 1985; Tribovillard et al., 2006; Algeo and Liu, 2020).

The Peace River Embayment was a prominent bathymetric feature of the WCSB at the end of the Devonian and served as a relatively deep-water connection between the open ocean and Prophet Trough to the basin’s interior (Robison, 1995). At Loc. 3 and 4, the C_org content increased from ∼0.4% in the Big Valley (BV) Fm. to 6.1% and 11.7%, respectively, at the base of the Exshaw Fm. (Figures 2A,B). The average C_org content at Loc. 4 decreased up-core from 5.7% ± 1.6% in the Exshaw Shale to 3.3% ± 0.2% in the Banff Fm. The δ¹³C_org values decreased from −28.4‰ in the BV to −28.6 ± 0.2‰ in the Exshaw Shale, −29.2 ± 0.4‰ in the Siltstone, and −29.1 ± 0.1‰ in the Banff Fm. At Loc. 3, a consistent C_org content of 5.2% ± 0.7% and δ¹³C_org of −28.8‰ was maintained throughout the shale and siltstone units. A similar published profile for the Peace River Embayment near Loc. 2 (Figure 1) reported δ¹³C_org values between −28.2 and −28.5‰ across the Exshaw interval (Caplan and Bustin, 1998).

The Loc. 11 core includes the Big Valley carbonate at the base, a very thin layer of Exshaw Shale (∼1.25 m) overlain by siltstone and Banff carbonate (Figure 2D). In comparison, the Exshaw Shale at Loc. 9 is 10.75 m thick without a siltstone member above it (Figure 2C). These two locations occur on the Madison Shelf, which was a shallow marine setting compared to the Peace River Embayment in the north and the Prophet Trough in the west (Smith and Bustin, 2000; Vernon, 2001). At Loc. 11, C_org increased from 0.4% just below the Exshaw Fm. to a maximum of 8.2% in the Exshaw shale. The average C_org content decreased rapidly thereafter, to < 2% in both the Exshaw Siltstone and Banff Fm. Average δ¹³C_org value increased from −29.3 ± 0.1‰ in the Exshaw Shale to −27.7 ± 0.4‰ to −27.0 ± 0.2‰ in the Banff Fm. At Loc.9, C_org increased from 0.5% in the BV Fm. to a maximum of 10.5% near the base of the Exshaw Shale. It then decreased to 0.8% in the Banff Fm. The δ¹³C_org values of Exshaw Shale at Loc.9 were similar to Loc. 11 despite differences in unit thickness. At Loc. 9, average δ¹³C_org values of Exshaw Shale was −29.1 ± 0.1‰ and rose to −28.5‰ in the Banff Shale.

The Banff Fm. sediments at Jura Creek (Loc. 21) and Crowsnest Pass (Loc. 22) were deposited along the eastern margin of the Prophet Trough. Here, the Palliser Fm. underlies the Exshaw Fm., and the Exshaw Shale is formally divided into Lower Shale and Upper Shale members based on their C_org richness, microfossil abundance and mineral content (Vernon, 2001). Both shale members consist of black laminated mudrocks, but the Lower Shale had a higher C_org content. At Loc. 21, the average C_org concentration for the Lower Shale was 3.9% ± 1.1% and the Upper Shale was 1.9% ± 0.7%. There was a drop in δ¹³C_org values from −28.3‰ to −28.9‰ in the Lower Shale, followed by a¹³C enrichment to −26.5‰ in the Upper Shale (Figure 2E). At Loc. 22, the C_org content for the Lower Shale ranged between 0.5% – 9.6%, with a mean C_org of 2.8% ± 2.4%, whereas for the Upper Shale, it ranged between 1.4% and 4.3%, and a mean C_org of 2.6% ± 1.1% was recorded. There was a notable maximum of 13.4% at the base of the siltstone unit (Figure 2F). δ¹³C_org values ranged from −29.4‰ at the base of the Exshaw to −27.2‰ in the siltstone member. The occurrence of a ∼3‰¹³C-enrichment at the transition from the Lower to Upper Shale at both Loc. 21 and 22 is notable, as this is the most significant carbon isotopic excursion detected in the WCSB.

Banff Assemblage δ¹⁵N profiles were broadly similar at all locations, exhibiting relatively low values (0.1–2.9‰) at or near the base of Exshaw Shale associated with high C_org content (Figure 2). The maximum δ¹⁵N values ranged between 4.5‰ and 7.2‰, corresponding with lower C_org near the top of the cored intervals. In the Peace River Embayment at Loc. 4, the δ¹⁵N minimum of 1.5‰ was at the base of the Exshaw Shale, and the average was 2.8 ± 0.6‰ for the Shale and 4.0 ± 0.9‰ for the Exshaw Siltstone. The overlying Banff Fm. had an average δ¹⁵N value of 4.9 ± 0.2‰. A previous investigation by Caplan and Bustin (1998) at a site close to Loc. 2 reported δ¹⁵N with a stratigraphic trend similar to that at Loc. 4, with a minimum of 0.5‰ at the base of Exshaw and increasing upward to a maximum of 5.1‰. At Loc. 3, there was a gradual increase in δ¹⁵N from 2.9‰ at the Exshaw base to a maximum of 4.2‰ in the siltstone unit. The average δ¹⁵N value was 3.2 ± 0.6‰ for Exshaw Shale and 3.7 ± 0.2‰ for Siltstone.

The two Madison Shelf locations recorded 3–4‰ negative δ¹⁵N excursions across the Big Valley/Exshaw Boundary and positive excursions of similar magnitude at the top of the Exshaw Shale (Figures 2C,D). Loc. 9 had a minimum δ¹⁵N value of 1.6‰ near the base of the Exshaw Fm., an up-core increase to 3.0‰ near the Exshaw top, and a maximum of 5.9‰ at the bottom of Banff Fm. The average δ¹⁵N for the Exshaw Shale was 2.4 ± 0.7‰. There was a positive trend in δ¹⁵N values up-core at Loc. 11 from 1.9 ± 0.3‰ in the Exshaw Shale to 7.0 ± 0.3‰ in the Banff Fm.

Locations 21 and 22 also had positive up-core trends in δ¹⁵N (Figures 2E,F), recording a ∼2‰ negative excursion at the Exshaw base, with a minimum of 2.1 and 2.2‰, respectively. Loc. 21 had an average δ¹⁵N of 3.7 ± 1.1‰ for the Lower Shale and 4.5 ± 0.8‰ for the Upper Shale. At Loc. 22, average δ¹⁵N values increased from 4.6 ± 0.3‰ for the Lower Shale to 6.2‰ for one sample in the Banff Fm.

Sediments from the base of the Exshaw Shale had an average C_org concentration of 7.9% ± 5.9% across 20 locations. There was significant variability among sites, with the greatest C_org content at Loc. 2 (19.9%) in the Peace River Embayment, Loc. 5 (18.6%) on the northern edge of the Rundle Shelf, and Loc. 17 (18.4%) on the Madison Shelf (Figure 3). Low C_org content was found at Loc. 6 (0.7%) on the Rundle Shelf and Locs. 15 (0.3%), 16 (0.2%), and 18 (2.0%) on the southwestern edge of the Madison Shelf adjacent to the Prophet Trough. The δ¹³C_org values ranged between −29.5 and −26.8‰ across the basin, with low C_org locations having relatively ¹³C-enriched (−27.8 ± 0.7‰) OM compared with the moderate (−28.9 ± 0.4‰) and high (−28.6 ± 0.6‰) C_org locations (Figure 4). The δ¹⁵N values at the base of the Exshaw Shale ranged from 0.0 to 6.3‰ and were inversely correlated with C_org content (Figures 3, 4). Locations with low C_org content were most ¹⁵N-enriched (δ¹⁵N = 6.0 ± 0.5‰) compared with locations with moderate (2.4 ± 1.1‰) and high (0.5 ± 0.1‰) C_org content.

The Exshaw Shale member has consistently higher concentrations of redox-sensitive TEs relative to the underlying and overlying strata (Supplementary Figure S1; Frucci, 2021). Trace elements are relatively immobile irrespective of the level of catagenesis (Ross and Bustin, 2008), so spatial trends in their enrichments reflect geographic variability in redox conditions during deposition and early burial (Pi et al., 2014; Algeo and Liu, 2020). We focused our analysis on samples from the base of the Exshaw Shale to compare spatial patterns in redox conditions with C and N stable isotope distributions during a single time interval when transgression created sediment accommodation space throughout the depositional basin (Figure 3). We found greater than two orders of magnitude difference in EFs across the basin (Table 1). The ratios of Mo, U, and Ni to Al for all sampled locations were greater than those of the reference PAAS, with maximum EFs of 3.0 × 10³ for Mo, 4.5 × 10² for U, and 1.3 × 10³ for Ni at Loc. 14 and similarly high EFs at Loc. 19 (Figure 3; Table 1). The Zn_EF and Co_EF were also elevated at Locs. 14 and 19, but at many locations Zn_EF and Co_EF was <1.0 indicating content similar to and below that of PAAS (Table 1). The V_EF was >1.0 at all locations except one, and it was greatest at Locs. 14 and 17. The Cr_EF ranged only between 0.52 and 3.14, and the Cr concentration was below the level of detection at Locs. 13, 14, 16, 18, and 19 in the south of the WCSB.

There is a general pattern of covariance among the TEs at the base of the Exshaw Shale as shown by the positive Spearman rank order correlation coefficients (r; Figure 5). Elevated Mo_EF and U_EF in black shales are commonly used to infer bottom water anoxia (Algeo and Li, 2020), and are strongly correlated at the base of the Exshaw Shale (r = 0.84, p < 0.001). Strong correlation with Mo_EF was also observed for Ni_EF, Zn_EF, and V_EF (Figure 5; Supplementary Figure S2), all with p < 0.001). The Cr_EF was weakly correlated with Mo_EF (r = 0.57, p = 0.033), and Co_EF was not significantly correlated with Mo_EF (p > 0.05). Devonian-Carboniferous black shales deposited in other North American basins are also highly enriched in Mo and U, with molar ratios of Mo to U ranging from 0.4 to 2 times the seawater molar ratios (Algeo and Tribovillard, 2009), and variable sedimentation rate and dilution by detrital minerals can explain variable enrichment of TEs deposited at that time (Algeo and Liu, 2020). The thickness of the Exshaw Shale interval in the WCSB, however, was not significantly correlated with any TE enrichment factors (Figure 5). Considering shale member thickness as a proxy for sedimentation rate, with a greater rate inferred for thicker intervals, the spatial variation in the degree of enrichment likely reflects differences in ocean redox conditions or circulation patterns across the basin. Paleogeographic reconstructions indicate that shallower regions on the eastern boundary of the Peace River Embayment, the eastern Rundle Shelf, and the southwestern Madison Shelf generally experienced lower TE enrichment (Figure 3D). Locations with greater enrichment of TEs, especially V and Mo (Figure 5), typically corresponded to greater C_org content (Figures 3C,D), though the relationship was relatively weak suggesting other controls such as primary productivity were more influential on C_org.

The onset of Exshaw Shale deposition was tied to the first of three transgressive events recorded in the Famennian-Tournaisian Banff assemblage (Smith and Bustin, 2000), which spans the D-C boundary in the WCSB (Macqueen and Sandberg, 1970). Crustal deformation associated with the Antler Orogeny, active during the Late Devonian and Carboniferous, produced the relatively complex seafloor bathymetry of the WCSB (Jewell, 1994; Robison, 1995). The intense tectonic activity affected the entire western extent of the North American craton, and in the WCSB it caused deepening of the Prophet Trough and parts of the cratonic platform and shelf margin (Poole, 1973; Savoy, 1990; Barclay et al., 1994; Smith et al., 1995). Crustal extension and vertical fault block movement also resulted in deepening of the Peace River Embayment (Caplan and Bustin, 2001; O'Connell et al., 1990). Thicker deposits of the Exshaw Shale accumulated on descending blocks as deeper locations had greater depositional accommodation (Smith et al., 1995; Smith and Bustin, 2000; Zaitlin, 2011; Visy, 2022). The difference in shale thickness and therefore sediment accumulation rate, for example, at Loc. 9 compared with Loc. 11, did not result in a dilution effect on geochemical characteristics (Figures 2, 5).

The complex bathymetry of the WCSB, which extended from the paleo-equator to ca. 10° N, likely affected the location of bottom water currents and regions of upwelling (Caplan and Bustin, 1998). On modern tropical western continental margins, interactions among wind stress, Ekman transport, and basin geometry can produce focused regions of surface water divergence and upwelling (Brandt et al., 2023). The WCSB likely included regions within the Peace River Embayment where upwelling was favored. As a modern geographical analog to the WCSB, the tropical seas between Indonesia and Australia have similar complex bathymetry. A study of the Arafura Sea modeled subsurface currents and related the complex bathymetry with the upwelling of nutrient-rich water near the shelf break, where there is elevated primary productivity (Kämpf, 2016). Subsequently, the newly upwelled waters become relatively nutrient-depleted, and these residual waters advect via surface currents to less productive regions.

Similar effects of bathymetry, currents, and upwelling on primary productivity can explain the spatial variability in C_org content of the Exshaw Shale. Greater C_org content is associated with eutrophic locations, and relatively oligotrophic (or less eutrophic) locations have lower C_org content (O’Connell et al., 1990; Richards, 1989). Based on Rock-Eval pyrolysis hydrogen index (Supplementary Table S1), OM in the Exshaw Fm. at all locations throughout the WCSB was principally derived from marine organisms. Short chain alkanes basin-wide also indicate the dominance of marine-derived OM (Frucci, 2021). The greater C_org content at the base of the Exshaw Fm. reflects increased primary productivity and export production coinciding with the transgressive stratal succession at the onset of black shale deposition (Caplan and Bustin, 2001; Visy, 2022). The broad distribution of marine sourced OM agrees with the previous analysis of biomarkers in oils derived from the Exshaw Fm. (Allan and Creaney, 1991). Further, accumulation of C_org-rich sediments in the Peace River Embayment has been previously connected to enhanced marine productivity in an upwelling zone (Caplan and Bustin, 1998; 2001). Previous geochemical analyses, however, focused on locations where the Exshaw Shale was most C_org rich, and our results expand this finding to show that sedimentary OM at all locations including those with relatively low C_org content has a marine origin.

The extent of bottom water anoxia during the deposition of Exshaw Shale at each location was assessed based on sedimentary facies classification and relative concentrations of trace elements. Trace fossil burrows in the Exshaw shale at Locs. 5, 15, 16, and 18 indicate greater oxygen exposure at the sediment surface during or soon after deposition on the Rundle Shelf adjacent to the Peace River Embayment and on the Madison Shelf adjacent to Prophet Trough (Visy, 2022). These locations also have elevated oxygen index values consistent with organic matter oxidation during early burial (Supplementary Table S1). Burrowing is interpreted to occur in shallower water environments where oxygenated conditions occur more frequently throughout the water column (Barclay et al., 1994; Richards et al., 1994; Smith et al., 1995; Smith and Bustin, 2000; Aretz and Corradini, 2021). Similarly, the TE enrichment factors were lowest at Locs. 3, 6, 8, 9, 16, 18, and 20, also indicating bottom waters in regions with burrows were more frequently oxygenated (Figures 3D, 5). While the entire basin was relatively shallow, the occurrence of subsidence by extensional faulting formed bathymetric lows that experienced more persistent stratification, thus those regions behaved as silled basins and were less prone to mixing between surface and bottom waters, similar to the modern Cariaco Basin (Thunell et al., 2004). Redox sensitive TE enrichment factors indicate that the deeper waters of the Peace River Embayment and Prophet Trough experienced oxygen depletion, whereas the paleobathymetric highs of the Rundle and Madison Shelf were comparatively oxygenated (Figure 5).

The nitrogen isotopic composition of marine sediments depends on the isotopic composition of phytoplankton and microbial transformations during transport to the ocean floor. Sinking particles are ¹⁵N-enriched by ca. 0.15‰/200 m, so the N isotopic composition of sediments in shallow basins with anoxic bottom water broadly reflects average phytoplankton (Uveges et al., 2019). The inverse relationship between δ¹⁵N and C_org at the base of the Exshaw Shale is intriguing (Figure 4), and the natural logarithmic fit (R ² = 0.83) suggests the possible influence of Rayleigh fractionation of N isotopes during phytoplankton uptake (Okhouchi and Takano, 2014). The data can be explained by advection of surface waters from upwelling locations with high nutrient concentration, productivity, and OM accumulation to locations with low nutrient concentration, productivity, and OM accumulation. Thus, the low δ¹⁵N values associated with higher sedimentary C_org would result from isotopic fractionation during incomplete assimilation of upwelled nutrients, leaving the residual nutrient pool ¹⁵N-enriched. The higher δ¹⁵N values associated with lower sedimentary C_org would derive from phytoplankton growth on the advected residual nutrient pool. The ¹⁵N enrichment of the residual nutrient pool might also be affected by dissimilatory nitrification/denitrification and anammox, with the cumulative effect of assimilatory and dissimilatory processes contributing to the shape of the Rayleigh fractionation curve (Waser et al., 1998; Sigman et al., 2001; Altabet, 2005).

Collectively, phytoplankton populations preferentially draw on available NH₄ ⁺ or NO₃ pools before N₂ fixation becomes widespread (Howarth et al., 1988). During oceanic anoxic events (OAEs), the global “Nitrostat” augments bioavailable N with increased N₂ fixation to replace N lost from the system due to increased denitrification and anammox (Kuypers et al., 2004). The concentrations of NH₄ ⁺, NO₃ ⁻, and NO₂ ⁻ (used as an electron acceptor in the anammox process) vary in the marine water column, depending on circulation, uptake by phytoplankton, ammonification during OM decomposition, and redox controls on microbial transformations (Sigman et al., 2001; Ohkouchi and Takano, 2014). In redox-stratified waters, nitrification, denitrification, and anammox are associated with the interface between reducing and oxidizing conditions and collectively diminish the nutrient-N pool. In deep waters below the photic zone, nutrient-N accumulates as NH₄ ⁺ in anoxic waters and NO₃ ⁻ where O₂ ventilation supports nitrification (Junium et al., 2018). The anoxic/suboxic conditions that predominated in the WCSB favored bottom-water NH₄ ⁺, which would constitute an internal N source that could support primary productivity when upwelled to the photic zone (Domingues et al., 2011; Peng, 2015; Naafs et al., 2019; Du et al., 2023).

In modern anoxic basins, the isotopic composition of bottom water NH₄ ⁺ is similar to that of surface sediments. In the Cariaco Basin, deep NH₄ ⁺ (δ¹⁵N = 4.6‰) is ¹⁵N enriched by 1.2‰ compared with surface sediments (δ¹⁵N = 3.4‰; Thunell et al., 2004). In the Black Sea, deep NH₄ ⁺ (δ¹⁵N = 1.7 ± 0.2‰; Velinsky et al., 1991) is ¹⁵N depleted by ∼1.2‰ compared with surface sediments (δ¹⁵N = 2.9 ± 0.6‰; Fulton et al., 2012). Smaller stratified basins like Framvaren Fjord (δ¹⁵N = 1.7‰ for deep NH₄ ⁺ and 2.2‰ for surface sediments; Velinsky and Fogel, 1999) and Fayetteville Green Lake (δ¹⁵N = 1.7‰ for deep NH₄ ⁺ and 3.2‰ for surface sediments; Fulton et al., 2018) have similar relationships. The average isotopic difference (δ¹⁵N_NH4+ - δ¹⁵N_sed) for these four representative stratified anoxic basins is −0.5 ± 1.2‰, and we apply this value in the calculation of δ¹⁵N of deep NH₄ ⁺ for the WCSB during the onset of Exshaw Shale deposition. For 20 locations, the average δ¹⁵N value for the base of the Exshaw Shale was 1.7‰ (weighted by C_org content). Therefore, the calculated average deep water δ¹⁵N_NH4+ value was 1.2 ± 1.2‰.

Unlike nutrient-N that can be fixed in situ from dissolved N₂, new phosphorus is mainly delivered to the ocean via continental weathering (Arthur and Sageman, 1994; Algeo and Scheckler, 2010). Sinking particulate matter transfers P along with N to deep waters and sediments. In the presence of O₂, N:P is typically close to the Redfield Ratio (16:1) for both particulate and dissolved components (Redfield, 1958; Falkowski et al., 1998; Tyrrell, 1999). In changing redox conditions, dissolved PO₄ ^3- concentration is regulated by adsorption and desorption from mineral phases (Furumai et al., 1989; Van Cappellen and Ingall, 1996; Ingalls et al., 2022). Under anoxic bottom waters, sedimentary P burial efficiency is reduced, leading to increased PO₄ ^3- concentration and decreased nutrient N:P, as has been observed in the Black Sea (Brewer and Murray, 1973; Fulton et al., 2012). Denitrification and anammox further contribute to lower N:P nutrient ratios. In such systems, upwelling of water with low N:P to the photic zone supports cyanobacterial N₂ fixation after nutrient-N has been exhausted.

For OAEs, modeled elevated PO₄ ^3- concentrations and enhanced denitrification contribute to elevated N₂ fixation. With a 2-fold increase in PO₄ ³⁻concentration associated with OAE 2, the modeled global contribution to export includes 55% N₂ fixation, 35% NO₃ ⁻ assimilation, and 10% NH₄ ⁺ assimilation (Naafs et al., 2019). Regions of high primary productivity were focused in upwelling zones where NH₄ ⁺ assimilation dominated, whereas N₂ fixation was prominent in locations receiving a low supply of bioavailable N. Martinez et al. (2019) reported ¹⁵N-enrichment throughout the Cleveland Shale interval in the Appalachian Basin for the D-C boundary, which they associated with increased denitrification rates during a transgressive event. Other D-C boundary sections reported negative δ¹⁵N excursions interpreted as indicative of increased N₂ fixation (Caplan and Bustin, 1998; Liu et al., 2016), following the common interpretation for similar ¹⁵N-depleted black shales in Cretaceous marine sediments (Kuypers et al., 2004; Ohkouchi et al., 2006; Junium and Arthur, 2007; Monteiro et al., 2012; Fujisaki et al., 2016; Ruebsam and Schwark, 2023).We propose spatial variability in nutrient availability within the WCSB could drive the observed differences in the N isotopic composition of Exshaw Shale sediments (Figure 3). Similar variability in N cycling has also been proposed for the D-C boundary associated with complex paleogeography and location-specific differences in the extent of denitrification across a shallow marine basin in South China (Liu et al., 2016). We consider the role of surface water advection as the primary driver of divergent sedimentary δ¹⁵N values.

The variability of δ¹⁵N within the WCSB (Figure 3), and its inverse relationship with C_org (Figure 4) may be most directly related to isotopic fractionation during incomplete assimilation of upwelled nutrients and subsequent surface water advection. This differs from the common interpretation that isotopic fractionation during denitrification and N₂ fixation is the direct cause of higher and lower δ¹⁵N values, respectively. Kinetic isotopic fractionation is quantified by the isotope enrichment factor, denoted by epsilon (ε). For N, ¹⁵ε reflects the ratio of the reaction rate constants (k) of light (¹⁴N) and heavy (¹⁵N) isotopes. ε 15 = k 14 / k 15 − 1 × 1000 ‰ (1)

This calculation yields a positive ¹⁵ε value for a normal kinetic isotope effect with ¹⁴N reacting faster than ¹⁵N. Isotopic fractionation between nutrient-N and biomass depends on the N source, which is affected by redox conditions. The ¹⁵ε for NH₄ ⁺ assimilation can vary widely from 4 to 27‰ depending on its concentration. In NH₄ ⁺ replete conditions as would be expected for OAEs, partial uptake with ¹⁵ε ≥ 20‰ can produce phytoplankton biomass with exceptionally low δ¹⁵N, leaving behind an exceptionally ¹⁵N-enriched NH₄ ⁺ pool (Waser et al., 1998; Higgins et al., 2012). The ¹⁵ε value for NO₃ ⁻ assimilation in modern sea surface water is ∼5‰ (Sigman et al., 1999).

The isotopic composition of phytoplankton biomass can be calculated by applying the Rayleigh equation for isotopic fractionation: δ org = δ init + f / 1 − f · ε · ⁡ ln f (2) where δ_org is the δ¹⁵N value for biomass N, δ_init is the δ¹⁵N value for the nutrient source, ¹⁵ε is the isotope enrichment factor, and f is the fraction of unutilized nutrient-N in the system after assimilation, i.e., the ratio of NH₄ ⁺/NH₄ ⁺-initial, and NO₃ ^-/NO₃ ^--initial. For high C_org (δ¹⁵N = 6.0‰), moderate C_org (δ¹⁵N = 2.4‰), and low C_org (δ¹⁵N = 0.5‰) locations (Figure 4), we assessed the range of potential δ_init values as a function of residual nutrient-N (0 ≤ f ≤ 1) for ¹⁵ε = 20‰ (a relatively conservative value for NH₄ ⁺ assimilation, Figure 6A) and ¹⁵ε = 5‰ (NO₃ ⁻ assimilation; Figure 6D).

Using Eq. 2 and isotope mass balance, δ¹⁵N of the residual nutrient pool (δ_resid) as a function of δ_init, δ_org, and f is determined as follows: δ resid = δ init − 1 − f · δ org / f (3)

The δ_resid values increase dramatically with greater assimilation of the initial nutrient pool, resulting in δ¹⁵N values much greater than any observed in the Exshaw sediments (Figures 6B,E). Assuming that the residual nutrient pool was advected away from the upwelling source and supported phytoplankton growth in less productive locations, the residual nutrient-N concentrations would decrease and the δ¹⁵N would increase even more. However, if the initial upwelled nutrient pool had an N:P ratio lower than the Redfield ratio as expected for anoxic bottom waters, and PO₄ ^3- was advected in the surface waters along with nutrient-N, then the nutrient N:P ratio of advected waters would become progressively lower (assuming uptake close to the Redfield ratio). Thus, the advected water mass would have excess P that could ultimately support N₂ fixation distant from the initial upwelling region. Given the maximum δ¹⁵N value near 6.0‰ for the base of the Exshaw Fm., cyanobacterial N₂ fixation in locations distal to upwelling would be required to lower the very high phytoplankton δ¹⁵N values predicted for assimilation of advected nutrient-N (δ_resid).

The ratio of N₂ fixation to total N assimilation (r) can be defined as: r = δ resid − δ org / δ resid − δ fix (4) where 0 ≤ r ≤ 1 and δ_fix = −1‰ for N₂ fixation. (Modified from Eqs 2 and 3; Ohkouchi and Takano, 2014). The range of values for r required to produce observed δ_org values for the Exshaw Fm. are shown in Figures 6C,F.

We consider three example scenarios to evaluate nutrient utilization based on the advection model for the WCSB (Figure 7). All three scenarios feature deep NH₄ ⁺ as the primary N-nutrient source in the WCSB and allow for N₂ fixation after NH₄ ⁺ has been consumed.

1. Deep-water with δ¹⁵N-NH₄ ⁺ = 1.2‰ is upwelled or mixed directly to the surface, where it is assimilated by phytoplankton. Residual NH₄ ⁺ is advected away from the region of focused upwelling to less productive locations (Figure 7A).

With the three scenarios, we estimated the initial δ¹⁵N values of nutrient-N that would support the production of biomass. We considered three characteristic locations in the WCSB with high (19.0% ± 0.8%), moderate (7.7% ± 2.6%), and low (0.8% ± 0.8%) C_org and average δ¹⁵N values of 0.5‰, 2.4‰, and 6.0‰, respectively, based on Figure 4. Using the advection model parameters, our goal was to determine which scenario fits best with our observed δ¹⁵N and C_org data and calculate a range of expected values for proportion of N fixation: assimilation (r). The lowest δ¹⁵N values corresponded to productive upwelling locations with minimal N₂ fixation.

Scenario 1: Using the weighted average δ¹⁵N_sed value of 1.7‰, the calculated δ¹⁵N of deep water NH₄ ⁺ in the WCSB is 1.2‰. In this scenario deep-water NH₄ ⁺ is upwelled to the surface with minimal isotopic fractionation and supplied a nutrient-N source to phytoplankton with a low initial δ¹⁵N value (Figure 7A). This scenario is analogous to the whole water-column mixing model of Uveges and Pearson (2023), which shows that low sedimentary δ¹⁵N values are a result of periodic overturning of a stratified water column in an anoxic ocean. Assuming a high concentration of upwelled NH₄ ⁺ was partially utilized, the resulting ¹⁵ε would be high. The most productive locations with highest C_org would require f < 0.01 to produce δ_org = 0.5‰ with δ¹⁵N_init = 1.2‰ and ¹⁵ε = 20‰ (Point X in Figure 6A), implying near complete consumption of NH₄ ⁺ at the site of upwelling. Thus, the advected residual NH₄ ⁺ pool would be very small and ¹⁵N-enriched (Point X in Figure 6B). To produce a maximum δ¹⁵N_sed value of 6.0‰, phytoplankton production would require 94% of N to be derived from N₂ fixation (Point X in Figure 6C).

Scenario 2: In a strongly stratified basin that inhibits deep mixing, NH₄ ⁺ diffuses to the chemocline and is utilized by annamox and nitrifying/denitrifying bacteria (Konovalov et al., 2008; Romaniello and Derry, 2010). This would increase the NH₄ ⁺ δ¹⁵N value of the residual pool. With strong stratification, the concentration of NH₄ ⁺ decreases upward through the chemocline, and NH₄ ⁺ δ¹⁵N values can range up 10–20‰ at and above the chemocline (Fulton et al., 2018). In scenario 2, we assume a δ¹⁵N value of 15‰ for NH₄ ⁺ that is assimilated by phytoplankton at the top of the chemocline (Figure 7B). At the most productive locations with the highest C_org, this would require f = 0.55 to produce δ_org = 0.5‰ with ¹⁵ε = 20‰ (Point Y in Figure 6A). The δ¹⁵N value of the residual NH₄ ⁺ pool would be predicted to be 27.1‰ (Point Y in Figure 6B). Using Eq. 4, complete assimilation of the residual NH₄ ⁺ along with additional 79% N₂ fixation would account for δ_resid = 27.1‰ and δ_org = 6‰ (Point Y in Figure 6C). This means that the advected NH₄ ⁺ residual pool would provide 21% of the N requirement in the most oligotrophic locations, and 79% N₂ fixation could make up for the N-deficit.

Scenario 3: During upward mixing of chemocline water, NH₄ ⁺ might be nitrified and partially denitrified as it is mixed into oxygenated surface waters, thereby producing a NO₃ ⁻ pool with relatively high δ¹⁵N (Figure 7C). We apply a value of 5‰, which is similar to observed surface water NO₃ ⁻ δ¹⁵N values in modern oxygen-limited basins (Brandes et al., 1998; Voss et al., 2001). In the WCSB, with δ_init = 5‰, the most productive locations with highest C_org would require f = 0.18 (Point Z in Figure 6D), i.e., 82% of the initial NO₃ ⁻ would be transported away from the upwelling/mixing region by advection. Consequently, the advected nutrients could support greater primary productivity, as most of the NO₃ ⁻ (82%) would be laterally advected. This could explain progressively higher phytoplankton δ¹⁵N values away from the source region, but not decreasing C_org. This scenario does not permit the introduction of significant N₂ fixation with r = 0.10 (Point Z in Figure 6F); thus, there is not a mechanism to replace nutrient-N lost to denitrification/anammox near the chemocline.

The advection model investigates the spatial variability of nitrogen isotopic distributions in the WCSB, assuming NH₄ ⁺ and NO₃ ⁻ as potential sources of bioavailable nitrogen in upwelling and stratified systems to account for a range of mixing scenarios. A component of N₂ fixation is present in all scenarios, but its extent varies with δ¹⁵N _initial and uptake ε values. Scenarios 1 and 2 represent the two end members of NH₄ ⁺ utilization by phytoplankton. The former depicts a scenario of mixing of surface and deep waters due to a “chemocline collapse” and therefore upwelled NH₄ ⁺ is isotopically light, whereas the latter demonstrates a situation of vertical stratification that results in a¹⁵N-enriched NH₄ ⁺pool being consumed by phytoplankton (Figures 7A,B). During stable periods, NH₄ ⁺ accumulates under the chemocline. When this stratification breaks down, the NH₄ ⁺gets released rapidly into the upper water layers. The third scenario holds true for modern OMZs, where the dissolved nitrogen pool is predominantly NO₃ formed by nitrification of NH₄ ⁺ as it reaches the surface (Figure 7C). Models predicting nutrient utilization during OAEs record effects of both NH₄ ⁺ and NO₃ ⁻ assimilation at different stages of redox stratification, contributing towards the final δ¹⁵N of exported biomass (Naafs et al., 2019; Uveges and Pearson, 2023). In the WCSB, an upwelling zone within the Peace River Embayment supported assimilation of a¹⁵N depleted nutrient pool (NH₄ ⁺). Portions of the Prophet Trough region and locations adjacent to the Sweetgrass Arch also have evidence for upwelled NH₄ ⁺. Across the interior parts of the epicontinental basin, on the Rundle and Madison Shelf, higher sediment δ¹⁵N values with low to moderate primary productivity suggest that N₂ fixation provided a greater portion of bioavailable nitrogen to sustain productivity. Without the influence of N fixation, the δ¹⁵N values would be even higher in the low-productivity continental shelf regions.

At the six locations with stratigraphic profiles representing shallow shelf and deeper Prophet Trough and Peace River Embayment provinces, there was a general relationship between δ¹⁵N values and C_org content (Figure 2). Low δ¹⁵N values are recorded in association with C_org maxima at the base of the Exshaw and gradually increase across the lower shale and upper siltstone deposits. An overall positive shift of 4–6‰ in δ¹⁵N is observed from the base of the Lower Shale to the contact between the upper siltstone and Banff Fm., possibly due to a redox switch from basin-wide anoxia during the onset of Exshaw deposition, to more oxygenated environments towards the upper Exshaw and Banff interval (Savoy, 1992). The End-Devonian OAE primarily sustained suboxic to anoxic conditions, based on TE data, and was probably more oxygenated relative to other major OAEs in the Phanerozoic. The decreasing TE enrichment up-core above the Exshaw Shale (Supplementary Figure S1) confirms more frequent oxygenation of bottom waters. Nitrogen isotope values, however, are primarily influenced by primary production, rather than by redox conditions as suggested by the strong correlation between δ¹⁵N and C_org and weak correlation with redox-sensitive trace elements from the base of the Exshaw Shale (Figure 5). This suggests that preservation or post-depositional bottom-water redox conditions have relatively little impact on nitrogen isotope signals.

The range of δ¹³C_org values at the base of the Exshaw Fm. (−26.8 to −29.5‰; Figure 4). indicates that the isotopic composition of OM was influenced by local factors in addition to the global atmospheric signature. A similar range of values was distributed stratigraphically through the Exshaw Fm., and δ¹³C_org was not correlated strongly with C_org either spatially or stratigraphically (Figures 2, 5). This range of values for the D-C boundary interval has been reported globally, including in the English River Formation in southeastern Iowa (Heath et al., 2021); the Kronhofgraben section in the Carnic Alps and Rhenish Massif in Germany (Kaiser et al., 2006; Agnieszka et al., 2020; Hartenfels et al., 2022); the Holy Cross mountains in Poland (Marynowski and Filipiak, 2007; Malec, 2013; Matyja et al., 2021), and the Cat Ba Island in Vietnam (Komatsu et al., 2014; Paschall et al., 2019; Shizuya et al., 2020). At the initiation of the black shale deposition, a shift towards lower δ¹³C_org values was observed in many D-C successions, including in the WCSB. Heath et al. (2021) reported an abrupt negative C_org isotope excursion at the onset of shale deposition, and a coincident increase in C_org content from the Burlington Shelf in southeastern Iowa. This is consistent with our sections at Loc. 9 and 11 on the Madison Shelf but not the deeper water locations in Prophet Trough, where δ¹³C_org decreased gradually in the shale, or in the Peace River Embayment, where there was no C isotope excursion in the shale unit (Figure 2).

Broadly, δ¹³C_org signals are a function of the isotopic composition of dissolved inorganic carbon (DIC), biological productivity, and growth temperatures. Local variations in the concentration and δ¹³C of the source DIC are especially influential, as δ¹³C_org values are sensitive to changes in the partitioning of ¹³C and ¹²C between organic carbon and carbonate (Velinsky and Fogel, 1999; Uveges et al., 2019). Microbially respired DIC that is ¹³C-depeleted when upwelled in epicontinental margins and assimilated by primary producers can influence the isotopic composition of the exported OM, causing lower δ¹³C_org values of marine sediment (Paschall et al., 2019; Pisarzowska et al., 2020). A broad trend of increasing carbonate δ¹³C values across the D-C boundary has also been recognized in association with the Hangenberg event (Brand et al., 2004; Kaiser et al., 2006; Myrow et al., 2013; Kumpan et al., 2015; Martinez et al., 2019; Heath et al., 2021). It has been interpreted as a global signal resulting from increased burial of ¹³C-depleted C_org and concomitant ¹³C enrichment of residual CO₂ in the atmosphere and DIC in the ocean (Buggisch and Joachimski, 2006; Bojar et al., 2013; Qie et al., 2016; Heath et al., 2021).

Lithologic variability within the Exshaw-Banff succession records variability in sea level during the D-C transition. An initial episode of transgression coincided with deposition of the C_org-rich Lower Shale of the Exshaw and a subsequent sea level fall and deposition of the upper siltstone (Macqueen and Sandberg, 1970). The siltstone unit was in turn followed by a second episode of transgression, and deposition of the calcareous Banff Shale. (Savoy, 1992; Robison, 1995; Figure 3, Supplementary Figure S4). This lithological sequence can potentially affect the occurrence of C isotope excursions associated with the Hangenberg event and/or the D-C boundary. In the WCSB, the D-C boundary occurs within the Exshaw Fm. (Caplan and Bustin, 1998). However, given the complex bathymetry of the WCSB, eustatic sea level change, and active tectonism at the time of deposition, it might not be definitively in the same lithologic unit at all locations. For example, at Loc. 21 and 22, the positive C isotope excursion appears to coincide with the Lower Shale/Upper Shale transition (Figure 2). At Loc. 11, there is a similar excursion at the Exshaw shale/siltstone transition. At Loc. 4, 9, and 11 there is a δ¹³C_org increase at the Exshaw/Banff boundary (Figure 2). The C_org content in all six locations with stratigraphic profiles record a decreasing trend within the Exshaw to the Banff succession (Figure 2). Carbon isotope values are generally stable in the shale or Lower Shale intervals; however, a shift to higher values and a correlative decrease in C_org is observed across the Exshaw-Banff boundary (Figure 2). A paleoclimatic shift from a greenhouse to icehouse climate may have been responsible for the observed trend in δ¹³C_org across the Exshaw-Banff boundary, since perturbations in the global carbon cycle are primarily related to atmospheric CO₂ levels (Becker et al., 2016a; Tuite et al., 2019; Zhang et al., 2021).

Nutrient cycling in past oceans might be analogous to the future, offering clues about how ocean ecosystems respond to shifts in nutrient availability. Environmental forcings that are currently causing coastal eutrophication and formation of oxygen minimum zones in modern oceans are mechanistically similar to forcings associated with some of the major marine biotic disruptions in the past (Paulmier and Ruiz-Pino, 2009; Zhou et al., 2016). The Devonian OAE was connected with the evolution of land plants and their extended root systems penetrating into the soil and amplifying the rates of continental weathering (Algeo and Scheckler, 2010). Consequently, there was an increase in the nutrient flux to the epicontinental seas, and globally widespread anoxia ensued. In modern oceans, warmer temperatures have intensified ocean stratification, which coupled with anthropogenic nutrient discharge in shallow seas has resulted in local stimulation of primary productivity and impinging coastal anoxia in many major continental margins and some restricted basins as well (De Pol-Holz et al., 2009; Domingues et al., 2011; Elfi Mollier-Vogel et al., 2012; Fulton et al., 2012; Brandt et al., 2023). Since continental margins are hotspots for marine biodiversity, prediction of nutrient circulation and transport pathways can help identify areas of the modern ocean that are particularly sensitive to nutrient imbalances. From this study we infer that high primary productivity is spatially tied to upwelling zones, with minimal local requirement for N₂ fixation. As with modern upwelling zones, the primary producers were likely eukaryotic algae. The parts of the basin distant from upwelling were more oligotrophic and required a greater contribution of N₂ fixation by cyanobacteria to balance the “nitrostat”.

Nitrogen isotope gradients in modern low oxygen environments like the western coasts of Peru and Ecuador suggest that advection patterns in nutrient utilization are characteristic of upwelling systems. The subsurface waters of the Eastern South Pacific coast host intensified OMZs due to high primary productivity and slow ventilation. As upwelling intensity varies spatially across the continental margin, it is reflected in sedimentary N isotope gradients that can be traced to different degrees of oxygen limitation (Elfi Mollier-Vogel et al., 2012). The Chilean margin also documents a strong latitudinal variability in sedimentary δ¹⁵N values, due to the lateral transport of ¹⁵N-enriched waters from an upwelling zone in the south of the OMZ towards higher latitudes (De Pol-Holz et al., 2009). The N-dynamics from modern upwelling zones confirm our interpretations of nutrient circulation across the D-C boundary in the WCSB. This phenomenon is not confined to the Devonian period alone, as anoxic oceans were prevalent throughout Earth’s history, and similar marine transgressive events have played a role consistently in facilitating the development of stratification and anoxia (Arthur and Sageman, 2004). Our interpretations of the influences of advection and upwelling on nitrogen isotopes in various depositional settings are not only applicable to past anoxic events but also enhance our understanding of nutrient circulation patterns in modern oxygen minimum zones. This knowledge is crucial for predicting how current and future changes in oceanographic conditions might impact marine ecosystems, particularly in regions susceptible to low oxygen levels. By integrating spatial variations in nitrogen fixation rates with nutrient transport and utilization models, we could potentially predict niche partitioning of eukaryotic algae versus cyanobacterial communities, allowing for a more comprehensive understanding of the interplay between nutrient availability and marine productivity patterns.