Dynamics of the Largest Carbon Isotope Excursion During the Early Triassic Biotic Recovery

Introduction

Following the Permian–Triassic boundary mass extinction (PTBME), and as a consequence of volcanic activity of the Siberian Large Igneous Province (SLIP), the global carbon cycle entered a protracted disequilibrium state (Payne et al., 2004; Galfetti et al., 2007a) spanning the entire Early Triassic. Prior to the SLIP, the carbon cycle was stable and controlled by a combination of local mechanisms distributed in time and space around Earth at the time (Bagherpour et al., 2017) excluding any consistent overriding and common control as was the case for the Early Triassic. This abrupt change of carbon cycle dynamics points to the crucial role of volatiles emitted from the SLIP in the disruption of the Late Permian equilibrium state. Based on the age of the Permian-Triassic boundary (PTB; Burgess and Bowring, 2015; Baresel et al., 2017) and of the early-middle Anisian boundary (Ovtcharova et al., 2015) this unstable state persisted for at least 5.4 Myr before waning during the middle Anisian. Until a new equilibrium state was restored during the Anisian, the recovery of the marine benthos lingered, thus markedly differing from the diversification-extinction crises of the nekton (Orchard, 2007; Brayard and Bucher, 2015) that were alternating with the ecological crises of terrestrial plants (Hochuli et al., 2016). The low competition within the shelly benthos most likely contributed to this delay (Hautmann et al., 2015). Diversified bivalve communities returned in the middle Anisian (Friesenbichler et al., 2019) while coral reefs diversity was not restored to levels equivalent to that of the late Permian diversity until the Ladinian (Stanley, 2003), i.e., about 6 Myr after resumption of a stable carbon cycle. Of paramount importance is to reconstruct the relative timing of the global excursion of the carbon cycle, of successive diversification-extinction cycles of the nekton, intercalated ecological swaps of land plants punctuated by fern spikes, of eustatic sea-level change, carbonate crises and deposition of black shales dominated by terrestrial organic matter (OM) on shelves, and of climate change.

Here we reconstruct the timing of the largest positive carbon isotope excursion (CIE) of the Early Triassic during the late Smithian and basal Spathian (Payne et al., 2004; Galfetti et al., 2007b). The magnitude of the associated Smithian-Spathian Boundary (SSB) extinction even surpassed that of the PTBME for the nekton and was the largest throughout the entire Triassic. High-precision U-Pb dating techniques (chemical abrasion, isotope dilution, thermal ionization mass spectrometry, and CA-ID-TIMS) are applied to single zircon crystals of volcanic ash beds intercalated with fossil-rich marine sediments in four sections of the Luolou Fm. (Nanpanjiang Basin, South China; Figures 1A, 2), herewith achieving a temporal resolution of the stratigraphic record at the 100 kyr level. Here we use a Bayesian age depth model to calibrate the duration of the CIE, of the black shale deposition, of the global SSB unconformity, and the age for the onset of cooling in low latitudes. This model permits (i) to detect hiatuses and changing sedimentation rate that both influence the shape and amplitude of the CIE, (ii) to exclude a volcanic winter as a cooling mechanism, and (iii) to evaluate other alternative cooling mechanisms such as CO₂ drawdown by silicate weathering or by the biological pump in the hypothetical frame of the cessation of a middle Smithian SLIP volcanic episode.

Figure 1. (A) Early Triassic paleogeographic map of the Nanpanjiang Basin (modified from Bagherpour et al., 2017). A- alluvial facies, B-shallow-water siliciclastics, C-carbonate platform, D- slope, and E-basin. Studied sections: 1 – Qiakong, 2 – Laren, 3 – Shanggang, and 4 – Lilong. (B) Sections with position and sample numbers of dated ash layers and carbonate carbon isotope records (δ¹³C data in Supplementary Table S3). Lithological units modified from Galfetti et al. (2008). Biostratigraphic age controls include ammonoid zones (white rectangles) and conodont Unitary Association zones (UAZ, colored; see on line material for further information).

Figure 2. Field photographs of the Smithian-Spathian boundary (SSB); white lines indicate the position of the boundary between units Vb and Va; (A) Qiakong section; topmost part of Daye Fm visible in lower left corner (person for scale); (B) Laren section; (C) Lilong section; (D) close-up of the Laren section, typical facies of unit IVb, with laminated black shale and rare, thin intervening gray limestone (Anasibirites multiformis/Xenoceltites pauciradiatus Zone); (E) Laren section, close-up view of the boundary between units IVb and Va; and (F) Laren section, ash layer CHIN10, note positive grading within the ash layer; for (E) and (F) hammer for scale.

Materials and Methods

Mineral Separation and CA-ID-TIMS U-Pb Dating

Detailed U-Pb age determinations were carried out on single zircon grains from volcanic ash beds of four sedimentary sections (Figure 1B); Laren (GPS coordinates 24°36′25.30″N 106°52′40.70″E); Shanggang (24°48′44.40″N 106°32′31.90″E); Lilong (24°54′25.69″N 106°32′30.37″E); and Qiakong (25°51′26.38″N 107°18′32.09″E). Zircon was recovered from each ash layer by crushing, milling and sieving to <250 μm prior gravity separation (Wifley Table), Frantz magnetic separation, and heavy-liquid separation using methylene iodide. From some ash layers zircon crystals were selected and mounted in epoxy resin and polished to reveal internal growth textures, followed by cathodo-luminescence (CL) imaging using a JEOL JSM7001F thermal field emission scanning electron microscope (SEM) at the University of Geneva (Supplementary Figure S4; CL images provided on request). CL images were used to select zircon crystals without inherited cores, major cracks, inclusions or disturbed zoning prior chemical abrasion. Selected zircon were annealed for 48 h at 900°C and washed several times with 3 N HNO₃ in 3 ml Savillex beakers before the partial dissolution for 18 h at180°C in 40% HF and trace HNO₃ in a pressurized Parr vessel containing fifteen 200 μl microcapsules, in order to minimize the effects of post-crystallization loss of radiogenic lead (Mattinson, 2005). Zircon grains of the Shanggang section were partially dissolved for only 12 h at 210°C. The chemically abraded zircon crystals were transferred in 3 ml Savillex vials and cleaned with 3 N HNO₃ before fluxing overnight at 80°C on a hotplate in 6 N HCl. Zircon were washed again several times in 3 N HNO₃. Single zircon crystals were loaded into clean 200 μl microcapsules, spiked with ∼4–6 mg of the EARTHTIME ²⁰²Pb–²⁰⁵Pb–²³³U–²³⁵U tracer solution (Condon et al., 2015; McLean et al., 2015), and dissolved in ∼70 μl 40% HF and trace HNO₃ at 210°C for 48 h. After dissolution, samples were dried down at 140°C and re-dissolved in 6 N HCl at 180°C for 12 h in 200 μl microcapsules, dried down again at 140°C, and re-dissolved in 3 N HCl. Uranium and Pb were separated in a single micro-column anion exchange chemistry and collected in 7 ml Savillex beakers with a drop of 0.035 M H₃PO₄ and dried down. Samples were loaded on single outgassed Re filament with a modified Si-gel emitter (Gerstenberger and Haase, 1997). Lead isotopes were measured on a Thermo Scienitfic TRITON thermal ionization mass spectrometer in dynamic mode using a MasCom discrete-dynode secondary electron multiplier with a deadtime of 22.5 ns. Lead isotope results were fractionation corrected using the measured ²⁰²Pb/²⁰⁵Pb ratio. All common Pb in the analysis was assumed to be procedural blank (²⁰⁶Pb/²⁰⁴Pb 17.6186 ± 0.3677, ²⁰⁷Pb/²⁰⁴Pb 14.7298 ± 0.4504, and ²⁰⁸Pb/²⁰⁴Pb 35.77 ± 1.07; 2σ). Uranium was measured as U-oxide in static mode on Faraday cups equipped with 10¹² Ω resistors or in dynamic mode on a MasCom secondary electron multiplier. Isotopic ratios were corrected for interferences of ²³³U¹⁸O¹⁶O on ²³⁵U¹⁶O₂ using an ¹⁸O/¹⁶O of 0.00205 and a sample ²³⁸U/²³⁵U ratio of 137.818 ± 0.045 (Hiess et al., 2012). Data were statistically treated with the Tripoli and U-Pb_Redux software (Bowring et al., 2011) ET2535 tracer calibration version 3.0 was used. All zircon ²⁰⁶Pb/²³⁸U dates were corrected for initial ²³⁰Th/²³⁸U disequilibrium using a Th/U of 3.5 ± 1 (2σ) for magma composition, augmenting their ²⁰⁶Pb/²³⁸U age by ∼80–100 kyr. All reported ages are weighted mean ²⁰⁶Pb/²³⁸U ages with associated 2σ uncertainties but without the systematic uncertainties associated with the tracer calibration and decay constants unless otherwise stated. We decided to use the U-Pb- and Th-disequilibrium corrected age of the youngest zircon population with a mean square of weighted deviates (MSWD) of around 1 (Wendt and Carl, 1991) as the best estimate for the deposition age of an ash layer.

In this study we re-dated two ash layers from Laren previously dated by Galfetti et al. (2007b) and Ovtcharova et al. (2006): CHIN10 (250.55 ± 0.51 Ma), and CHIN40 (251.22 ± 0.2 Ma) were re-dated applying state of the art analytical techniques and EARTHTIME ²⁰²Pb-²⁰⁵Pb-²³³U-²³⁵U tracer solution, that allows a direct comparison to ages published for the PTME (Burgess et al., 2014), and EMTB (Lehrmann et al., 2015; Ovtcharova et al., 2015). Besides the significantly improved uncertainty on the weighted mean ages, we notice a million-year scale age offset that is outside analytical uncertainty. Ash layer CHIN40 immediately underlies the Flemingites beds (Figure 1) and is now dated ∼0.57 ± 0.23 Myr younger compared to Galfetti et al. (2007b) CHIN10 and the associated Tirolites/Columbites beds are now dated ∼1.67 ± 0.51 Myr younger. The offset in age may partly arise from the larger uncertainty on sample CHIN10 dated by Ovtcharova et al. (2006) which masks the natural age variability within the sample that is resolvable with our higher precision, due to revised isotope tracer calibration, as well as the improved chemical abrasion procedures that allow to discard lead loss as a potential cause for “too young” individual dates.

Age Depth Model

The age depth-model was obtained by the “rbacon” code (Blaauw and Christen, 2011; complete rbacon code in Supplementary material). The best estimate U-Pb date for the deposition age of an ash layer serves as input date for the age depth model. Zircon single grain analyses which are significantly younger or older than the 2σ uncertainty of the best estimate of ²⁰⁶Pb/²³⁸U age and/or further contradict the stratigraphic order of the volcanic ash beds were excluded from the final age-depth model and assigned to either (i) insufficiently mitigated loss of radiogenic lead, or (ii) significantly longer residence time of zircon in the magmatic reservoir and/or sedimentary recycling. The “rbacon” code uses a semiparametric approach to model the age depth relationship and a non-Gaussian autoregressive process to estimate sediment accumulation (Blaauw and Christen, 2011). Additionally, the code allows inclusion of a priori known observations, e.g., facies changes, or discontinuities. We assumed a hiatus of a maximum duration of 500 kyr between Unit VIb and Va for each section, see Supplementary Table S2 for detailed model settings and outputs. The default t-distribution of dates (optimized for radiocarbon dates) was changed to a normal distribution, otherwise the U-Pb dates were considered as outliers by the algorithm. We tested the effect of changing the settings of the prior for the calculation of the accumulation rate (e.g., changing the shape of the accumulation rate to 1.5; #acc.shape; see rbacon code in the Supplementary Material) as well as changing the segment length. The two parameters of the prior distribution of the accumulation autocorrelation distribution were changed from the default 0.7 and 4 value to 0.2 (#mem.mean) and 1.5 (#mem.strength), assuming a low memory effect. The segment length was adjusted to the smallest possible to allow the model to run. The accumulation rate was inferred from the U-Pb dates and stratigraphic position (#acc.mean; Supplementary Table S2). The two parameters with the biggest impact on the final age-depth model are the stratigraphic position and the memory effect. We inferred a hiatus between Unit VIb and Va in each section. If no such hiatus is included, the model showed an unrealistically high age discrepancy between the sections and some ash layers turn out as outliers.

Carbon Isotopes and Oxygen Isotope Measurements of Carbonates

High-resolution sampling was carried out in Qiakong, Shanggan, and Laren (Figure 1 and Supplementary Table S3) sections for stable isotope measurements (C and O) of bulk micrite (δ¹³C_carb, δ¹⁸O_carb). Samples were carefully cleaned, cut and drilled with a diamond-tipped drill to produce a fine powder. Carbonate carbon and oxygen isotopes compositions (δ¹³C_carb and δ¹⁸O_carb) were performed at the Institute of Earth Surface Dynamics, University of Lausanne with a GasBench II connected to a Thermoquest Finnigan DeltaPlus XL mass spectrometer, using a He-carrier gas system according to a method adapted after Spötl and Vennemann (2003). For calcite, a reaction temperature of 70°C was used and samples were reacted for 1 h. In-house standards of calcite were treated in the same way run interspersed with the samples in the same sequence. Samples were normalized using the in-house Carrara Marble standard calibrated against δ¹³C and δ¹⁸O values of NBS-19 (+1.95 and -2.20‰, respectively, relative to VPDB). External reproducibility for the analyses estimated from replicate analyses of the in-house standard (n = 8 per run of 32 samples) was better than ±0.08‰ for δ¹³C and ± 0.1‰ for δ¹⁸O values. The results are expressed in the δ-notation [δ = (R1/R2-1) × 1000] where R1 is the ¹³C/¹²C or ¹⁸O/¹⁶O ratio in the sample and R2 is the corresponding ratio of the standard V-PDB, in parts per thousand (‰).

Organic Geochemistry and Palynofacies

The characterization and quantification of preserved OM were performed on powdered whole rock samples at the Institute of Earth Sciences of the University of Lausanne, Switzerland, using a Rock-Eval 6 and following the method described by Behar et al. (2001). The samples were placed in an oven and first heated at 300°C under an inert atmosphere and then gradually pyrolysed up to 650°C. After the pyrolysis was completed, the samples were transferred into another oven and gradually heated up to 850°C in the presence of air. The Rock-Eval 6 results include: total organic carbon content (TOC, wt.%), hydrogen index (HI, mg HC/g TOC, and HC = hydrocarbons), oxygen index (OI, mg CO₂/g TOC), and T_max values (°C, OM thermal maturity indicator). The HI and OI parameters are used to distinguish between algal/bacterial (type I and II kerogen) and terrestrial and/or reworked (type III and IV) OM. The T_max, OI, and HI parameters can only be interpreted for TOC values ≥ 0.25 wt.% (Espitalie et al., 1985). Palynofacies analysis has been performed on nine samples to cover the entire positive CIE at the University of Zürich. A total of nine samples have been cleaned, crushed and weighed (∼13–18 g) and subsequently treated with hydro- chloric and hydrofluoric acid according to standard palynological preparation techniques (Traverse, 2007). Oxidation was performed using Schulze’s reagent (Traverse, 2007). The remains were sieved over an 11 μm mesh screen. For detailed method description see Hermann et al. (2011) and references therein.

Results

Biochronology

Accuracy of biochronological ages is of uttermost importance in constructing a reliable and robust frame in time and space for the biotic and abiotic events associated with the end-Smithian extinction. This goal has been facing two major obstacles. Firstly, the traditional use of conodont interval zones whose bases are defined by first occurrences (FO) of index species and whose tops are defined by the base of the next zone. These continuous zones assume that the FO of index species are synchronous, which is more the exception than the rule. Secondly, the incompleteness of the record, which has been generally overlooked but which opens a window in terms of global eustatic sea-level changes. Triassic ammonoid biochronology hinges on the principle of discontinuous and mutually exclusive maximal associations zones, which is an approach also shared by the more elaborate Unitary Association zones (UAZs; Guex, 1991, 2011; Hammer, 2013; Monnet et al., 2015). Correlation by means of such discontinuous biochronozones also directly bears directly on the reconstruction of the C-isotope signal, which can be distorted by gaps and whose absolute values may vary laterally (see review by Zhang et al., 2019 for the SSB). Moreover, any time line (e.g., the SSB) defined by an interval zone comes with no uncertainty, but with frequent contradictions (crossing over of FOs) between different sections. In contrast, the time interval containing a time line comes with a build-in uncertainty (i.e., the interval of separation between two zones Chen et al., 2019 for a SSB example) that can be subsequently narrowed down (by acquisition of additional data) but with no contradiction (crossover) when defined by maximal association zones and UAZs (see Brosse et al., 2016 for a PTB example).

Ammonoid zones based on maximal associations are here updated and new conodont Unitary Association Zones for the Nanpanjiang Basin are presented in Figure 1B. The most complete succession of ammonoid zones of the Luolou Fm. is documented in Laren (Galfetti et al., 2007a). However, the Luolou succession consistently misses a single zone between the late Smithian Glyptophiceras-Xenoceltites (GXZ) and the earliest Spathian T. n. gen A. (TAZ) zones. This new zone (NZ) is of late Smithian affinity has so far only been documented in NE Nevada where it intercalates between GXZ and TAZ (Bucher, ongoing work). This zone is recorded in a peculiar tectonic setting where high subsidence and sedimentation rates compensated for a global eustatic sea level fall and the associated hiatus (see detailed discussion of the biochronology in the Supplementary Material). Late Smithian ammonoid faunas have a cosmopolitan distribution, thus providing a superb tool for global correlation (Brayard et al., 2006). The new conodont UAZs bracketing the SSB provides an independent and additional time control in the Luolou Fm. (Supplementary Figures S1, S2). With the exception of NE Nevada, the present ammonoid/conodont biochronological frame for the SSB in the Luolou Fm. is the most highly resolved one, and is diagnostic of the low latitude record (Hammer et al., 2019). In Qiakong and Lilong, the SSB is bracketed by conodont UAZ5 and UAZ6 and falls in the separation interval between GXZ and TAZ. TAZ unequivocally records the base of the evolutionary radiation of tropical Tirolitidae in addition to the cosmopolitan Bajarunia euomphala. All relevant biochronologic data and processing are given in the Supplementary Material.

Carbon Isotope Compositions and Organic Matter Contents

Carbon isotope variations in middle Smithian to early Spathian carbonates obtained from the four studied sections (Figure 1) are consistent and synchronous with other SSB sections worldwide as established by ammonoid biochronological control. The four analyzed sections yield a spike-shaped positive CIE peaking either slightly below or at the boundary between the black shales (Unit IVb) and the nodular limestone (Unit Va). The particular shape of his CIE is known from almost all low latitude SSB sections and its variable amplitude suggests truncation of the CIE peak, thus pointing to a global hiatus (Hammer et al., 2019). In Qiakong, where sedimentation rates are higher and the record more complete, the CIE spans the upper half of the black shales and peaks at a δ¹³C value of +4.5‰. The onset of the CIE is taken where values depart from the preceding -2‰ plateau. This onset is gradual in Qiakong, whereas it starts abruptly at the base of the black shales in Laren and Shanggan. Similar to Qiakong, the CIE peaks at the top of the black shales in Laren and Lilong, but with values of only +3‰ and +3.5‰, respectively. The end of the positive CIE is taken where values first reach 0‰, a value defining the next negative plateau of middle Spathian age and that extends into the Hellenites Zone (Galfetti et al., 2007a). A value of 0‰ is first reached at the boundary between TCZ and PZ in Laren (Figure 1).

Carbon isotope compositions of organic carbon, Rock-Eval, TOC, and palynofacies analyses from Qiakong are presented in Supplementary Material, Supplementary Figure S3 and Supplementary Table S4, and provide evidence for the nature of the OM. Absence of major change in the composition of the OM supports the primary origin of the carbon isotope compositional changes. The parallel trend between organic and carbonate δ¹³C values suggests that variations within the carbon cycle were equilibrated over the atmosphere, linking the organic records to the inorganic records. Rock-Eval analyses plot within the fields of kerogen type II and III and point to a dominantly terrestrial origin of OM. The mature late Smithian OM is dominated by woody particles (both translucent and opaque) in this shelf setting, with a weak upward increasing trend, indicating a preeminent terrestrial origin. This late Smithian OM may conceivably be reworked from older strata or result from bacterial decomposition. The black shale unit is also accompanied by a reduction of the biogenic carbonate fraction in the sediment, suggesting a carbonate crisis on shelves.

U/Pb Geochronology of Zircon From Volcanic Ash Beds

Zircon U/Pb dates were analyzed from volcanic ash beds of all four sedimentary sections shown in Figure 1B. The complete data set can be found in Supplementary Table S1, and is shown in Figure 3. The choice of the relevant zircon dates for the calculation of the mean age is an essential and somewhat non-objective step. Below, we are giving our best estimates for the depositional age of a give ash layer, on the basis of the youngest cluster of zircon dates, after removal of some clear outliers from radiogenic lead loss. This “best estimate” age is compared to the age of the youngest zircon (lead loss outliers removed) and the mean age estimated from the maximum number of samples that still satisfy the statistical criteria of Wendt and Carl (1991) in Supplementary Table S5. Computing age-depth models with these different selection approaches will give an indication of the error introduced by this subjective choice (see below).

Figure 3. ²⁰⁶Pb/²³⁸U single-grain zircon analyses of volcanic ash layers from sections Qiakong, Laren Shanggang, and Lilong. MSWD = mean square of weighted deviates. Horizontal bars: single-grain zircon analysis (2σ error). Weighted mean ages are calculated from grains marked in black, analyses marked in gray are excluded from the weighted mean age calculation. Vertical bar: weighted mean age and 2σ error (x = internal, y = external uncertainty including tracer calibration, z = external uncertainty including tracer calibration, and ²³⁸U decay constant uncertainty).

Laren

Five ash layers from Laren were sampled from the latest Smithian to the early Spathian. Twenty-two zircon grains of ash layer CHIN10 were analyzed. Ash layer CHIN10 is a prominent, dm thick, coarse-grained volcanic ash layer thinning upward, which is interbedded within the marly limestone ca. 8 m above the SSB. This ash layer within the Luolou Fm. has a very large lateral distribution throughout Guangxi and southern Guizhou. Long prismatic zircon without any prominent feature or inheritance yielded the youngest U-Pb ages. A subset of zircon grains has significantly older ²⁰⁶Pb/²³⁸U ages (up to 1013 Ma), and exhibits sector zoning and/or xenocrystic cores (Supplementary Figure S4). Finally, three grains were considered to calculate the best estimate ²⁰⁶Pb/²³⁸U date of 248.853 ± 0.086/0.11/0.29 Ma 2σ (MSWD = 0.51; X/Y/Z error notation according to Schoene et al., 2010). Sample 230214AT is located at the lithological boundary between the nodular and wavy limestone 15 cm below ash layer CHIN10. Zircons are short to long prismatic. Sixteen grains were dated, the youngest three yield a final ²⁰⁶Pb/²³⁸U age of 248.86 ± 0.23/0.33/0.35 Ma 2σ (MSWD = 0.06, n = 3). Nine grains of ash layer 230214BT (2.75 m below CHIN10) were analyzed, the best estimate is a ²⁰⁶Pb/²³⁸U date of 248.97 ± 0.16/0.16/0.31 Ma 2σ (MSWD = 0.89, n = 3). Six grains of ash layer LAR206 (0.5 m below the base of the nodular limestone) were dated, three grains yielding the best estimate ²⁰⁶Pb/²³⁸U date of 249.35 ± 0.13/0.15/0.30 Ma 2σ (MSWD = 0.079, n = 3). Eight grains of CHIN40 were analyzed, four of which yielded a best estimate ²⁰⁶Pb/²³⁸U date of 250.647 ± 0.064/0.09/0.28 Ma 2σ (MSWD = 0.84, n = 4).

Shanggang

Five ash layers were sampled from the middle-late Smithian to the early Spathian. Ash layer SHA301T is ca. 8.4 m above the base of the nodular limestone, at the top of Unit Vb. Six grains were analyzed, four of them yielding a best estimate ²⁰⁶Pb/²³⁸U date of 248.86 ± 0.11/0.13/0.30 Ma 2σ (MSWD = 0.21, n = 4) for the ash layer deposition. Ash layer SHA306T is located ca. 2.7 m above the base of the nodular limestone, yielding a variety short to long prismatic zircon grains. The age of six grains was determined, of which five yielded a best estimate ²⁰⁶Pb/²³⁸U date of 249.071 ± 0.092/0.11/0.29 Ma 2σ (MSWD = 0.66, n = 5). Because of low angle faulting the topmost 0.6 m of the black shale below the nodular limestone are crushed. Ash layer SHA359T is located ca. 0.8 m below the base of the nodular limestone. Six grains were dated yielding a best estimate ²⁰⁶Pb/²³⁸U date of 249.33 ± 0.10/0.12/0.29 Ma 2σ (MSWD = 0.25, n = 6). Ash layer SHA339T is ca. 5.8 m below the base of the nodular limestone. Six long prismatic zircon grains were analyzed yielding a best estimate ²⁰⁶Pb/²³⁸U date of 249.74 ± 0.12/0.13/0.30 Ma 2σ (MSWD = 0.62, n = 6). Ash layer SHA303T is ca. 9.3 m below the base of the nodular limestone. The age of eight long-prismatic zircon grains was determined, seven of them yielding a best estimate ²⁰⁶Pb/²³⁸U date of 250.116 ± 0.085/0.11/0.29 Ma 2σ (MSWD = 1.1, n = 8).

Qiakong

The sample collection from the Qiakong section is composed of six thin ash layers that do not show signs of post-depositional lateral transport. Twelve long prismatic (∼200 μm) zircon of ash layer QIA13T are analyzed (ca. 14 m above the base of the nodular limestone). Zircon displays no inheritance at the resolution we are able to achieve. The youngest grain yielding a ²⁰⁶Pb/²³⁸U date of 247.83 ± 0.17 Ma was excluded from the final age calculation due to suspected loss of radiogenic lead. The best estimate ²⁰⁶Pb/²³⁸U date of QIA13T is 248.250 ± 0.065/0.091/0.28 Ma 2σ (MSWD = 1.2, n = 3). Sample QIA12T is located at the lithological boundary between the nodular limestone (unit Vb) and the overlying fine laminated marly limestone (unit Vc) and could possibly correlate with ash layer 230214AT in the Laren section as well as with SHA301T in the Shanggang section (Figure 1). Nine long prismatic zircon grains from QIA12T (∼200 μm) were analyzed. The youngest cluster of three grains yield a ²⁰⁶Pb/²³⁸U age of 248.893 ± 0.069/0.094/0.28 Ma 2σ (MSWD = 0.56, n = 3). Eleven zircon grains of ash layer QIA11T (ca. 6.3 m above the nodular limestone) were analyzed zircon are exceptionally low in radiogenic lead. The best estimate ²⁰⁶Pb/²³⁸U was calculated on the basis of the two youngest grains at 248.95 ± 0.49/0.50/0.57 Ma 2σ (MSWD = 0.0056; n = 2). The lithological boundary between the black shale and the nodular limestone is bracketed by ash layers QIA9T (²⁰⁶Pb/²³⁸U = 249.110 ± 0.088/0.11/0.29 Ma 2σ; MSWD = 1.8, n = 7) and QIA07T (²⁰⁶Pb/²³⁸U = 249.292 ± 0.063/0.091/0.28 Ma 2σ; MSWD = 0.36, n = 4), which are considered to be representative of zircon crystallization and ash bed deposition. Both younger zircon of ash layer QIA9T (²⁰⁶Pb/²³⁸U = 246.40 ± 0.018 Ma) and QIA07T (²⁰⁶Pb/²³⁸U = 248.58 ± 0.22 Ma), respectively, most likely suffered from a small amount of lead loss and were excluded from the age calculation. QIA3T is a reddish thin layered ash bed occurring 4 m below the base of the nodular limestone (Figure 1B). Twelve grains have ²⁰⁶Pb/²³⁸U dates from 253.4 ± 1.8 Ma to 371.1 ± 4.2 Ma, which are all considered to be inherited and are excluded from the age model. These grains of Triassic to Permian age are evidently reworked into the late Smithian black shale.

Lilong

Ash layer LIL508 is a thin ash bed occurring within the black shale ca. 0.5 m below the base of the nodular limestone. Zircons are short to long prismatic (ca. <100 μm). Twelve grains were measured, yielding dispersed ²⁰⁶Pb/²³⁸U dates over an age range from 249.56 ± 0.26 Ma to 263.38 ± 0.36 Ma. Therefore, all these grains are of detrital origin and are excluded from the age model. Interestingly, QIA3T, and LIL508 are the only two samples containing reworked zircon grains and are from unit Vb, providing another line of evidence that the black shale was deposited during a lowered based level accompanied by erosion in the hinterland.

Age-Depth Model

The stratigraphic position of each ash bed and the selected mean ages of the youngest cluster of zircon dates were used to calculate separate age-depth models for Qiakong, Laren, and Shanggang (Figure 4). The model of each section reveals an increase of sedimentation rate across the SSB and the presence of a hiatus at or near the top of the late Smithian black shale. Based on the ammonoid biochronological definition, the SSB is within the interval of separation between GXZ and TAZ (since NZ is missing in South China) in the topmost part of unit IVb, which is bracketed by ashes QIA07T (249.292 ± 0.063 Ma) and QIA09T (249.110 ± 0.088 Ma). The shortest estimate for the duration of the SSB gap in the most expanded Qiakong section derived from the age depth model is of 64 ± 104 kyr, a figure well within the range of Cretaceous “cold snaps,” and ephemeral ice sheets (Miller et al., 2005). The base of the TAZ is considered equivalent to the base of the nodular limestone (Figure 1B) and therefore yields a minimum estimate for the age of the SSB (249.167 ± 0.087 Ma; Figure 4). Calculating age-depth models on the base of (i) youngest zircon date of each sample, and (ii) minimum cluster at maximum acceptable MSWD value (see above) yields average age estimates of 249.133 ± 0.132 Ma and 249.263 ± 0.101 Ma, respectively, for the SSB. The three approaches thus yield coincident SSB model age estimates within their respective uncertainties.

Figure 4. Bayesian “rbacon” age models for Qiakong Laren and Shanggang. The red line represents the model mean age with its error envelope at 95% confidence interval in gray. Data are reported as Th-corrected ²⁰⁶Pb/²³⁸U ages with 2σ internal uncertainty. Two groups of ash layers yielded ages that overlap within error, thus allowing the positioning of two time-lines: a first line including QIA07T, LAR206T, and SHA359T at the top of the late Smithian Unit IVb, and a second line linking QIA12T, CHIN10, and SHA301T of late early Spathian age at the lower boundary of Unit Vc. The position and age of the “Black Band Marker,” an isochronous 10–15 cm thick shale layer 1.0 to 3.5 meters above the boundary between units Va and Vb in Qiakong (Figure 1) is indicated as well. The base of the nodular limestone rests upon an unconformity and is interpolated from all three sections at an age of 249.167 ± 0.087 Ma.

The bracketed age of the SBB is at marked variance with the ca. 1 Myr younger age derived from astronomical tuning (Li et al., 2016a, b). Omission of stratigraphic gaps and the utilization of interval zones defined by frequently diachronous first appearance of index conodont species provide potential explanations for this younger age. Interestingly, this divergence in the age of the SSB decreases when astronomical tuning is based on carbon isotope composition (Fu et al., 2016). Based on Qiakong, a duration of 622 ± 137 kyr is obtained for the deposition of the black shale. With an offset of 135 ± 163 kyr between the initiation of the black shale deposition and that of the CIE (Figures 5, 6), any quantification of the delayed response time of the δ¹³C record with respect to the onset of black shale deposition remains elusive. The preserved part of the positive limb of the CIE in Shanggang and Laren amounts to 296 ± 161 kyr and 251 ± 143 kyr, respectively, (Figure 5). Here again, the extremely short duration (19.21 kyr) of the positive limb of the CIE derived from astronomical tuning of carbon isotope composition (Fu et al., 2016) is perplexing. The age model also allows us to quantify the duration of each UAZ and their interval of separation in each section (see Supplementary Table S2). Effective cooling in the tropics only started with the GXZ (Goudemand et al., 2019) but it is conceivably expected to have started somewhat earlier at high latitude. This late Smithian cooling was also accompanied by the onset of the ecological recovery of terrestrial plants such as gymnosperms (Hochuli et al., 2016). The durations of the CIE and of the cooling phase presented in this study exclude a volcanic winter generated by sulfur-rich aerosol shielding because of the short residence time (months to years) of sulfur-compounds in the atmosphere and stratosphere (Barnes and Holmann, 1997; Self et al., 2006; Jones et al., 2016). However, sulfur emissions associated with protracted pyroclastic activity of flood basalt may have triggered longer periods of cooling (Black et al., 2012, 2018). As argued by these authors, a rather short punctuated cooling event embedded in a prolonged period of volcanic activity is unlikely to be resolved within the currently achieved precision of U-Pb dating. The apparent duration of the SSB unconformity in Qiakong, Laren and Shanggang is of 64 ± 104, 149 ± 143, and 127 ± 137 kyr, respectively, and excludes large scale plate tectonic controls that operate at a much longer time scale. The consistent biochronological age of this worldwide gap (Galfetti et al., 2007a, b; Hammer et al., 2019) is also at variance with the timing of regional block faulting tectonics within the Nanpanjiang Basin (e.g. Sun et al., 2015). Only glacio-eustatic sea level changes (Miller et al., 2005, 2011) are in agreement with the timing of this global unconformity. In Qiakong, the model also yields 95 ± 122 kyr for the part of the ascending limb of CIE comprised between +2 and +4.5‰ (pink interval, Figure 5). This interval is partly truncated by the longer SSB unconformity in Laren, Shanggang, and Lilong, thus explaining the smaller amplitude of the CIE. Interestingly, in Jiarong III, another deep water section located at the tectonically active northern rim of the basin (Sun et al., 2015), the amplitude of the CIE compares with that of Qiakong. Hence, the precise paleogeographic location within the basin has obvious consequences for the completeness of the CIE. Ash samples within the black shales (LIL508T, QIA3T; Figure 3) typically contain detrital zircons of Early to Late Permian age, indicating a lowered base level accompanied by erosion.

Figure 5. Durations of the hiatus, of the carbon isotope excursion (CIE) and age of the early Spathian “Black Band marker.” Color code: pink represents the part of the CIE in the Qiakong section that is missing in Laren and Shanggang; blue: duration of preserved carbon excursion in Laren and Shanggang; white: duration between onset of nodular limestone and black shale marker; and purple: duration of pre-CIE black shale deposition in Qiakong. The Black Band Marker provides an independent isochronous time line. ²⁰⁶Pb/²³⁸U single-grain zircon ages of dated ash layers are indicated, for sample numbers see Figure 1.

Figure 6. Conceptual model for the timeframe of the subsequent five-phase evolution of global climate, sea-level and ammonoid diversity from the mid-Smithian to basal Spathian, relating the carbon isotope record of this study to δ¹⁸O_phosN (Goudemand et al., 2019) and associated consequences and feedbacks.

Equipped with this calibration, the time compatibility of two cooling feedback mechanisms, i.e., silicate weathering and marine organic carbon cycling, can now be quantitatively evaluated. The results of these calculations are summarized in Table 1. Feedback mechanisms could conceivably account for some or all of the late Smithian cooling. Such assessments rest on the unwarranted assumption that the excess of atmospheric CO₂ originated from a volumetrically important middle Smithian eruptive burst of the SLIP.

Table 1. Mass balance calculations.

Discussion

Mechanisms of CO₂ Removal From the Atmosphere

Silicate Weathering

The weathering of large subaerial exposures of basaltic flows emplaced during the activity of Large Igneous Provinces can act as an efficient CO₂ sink during times of warm and humid climate. We adopt present-day CO₂ consumption rates of basaltic provinces from Dessert et al. (2003) and a minimal surface area (2.5 × 10⁶ km²) of the SLIP from Fedorenko et al. (1996) for a C cycle box model (Table 1C). Considering a 600 kyr duration for the positive limb of the late Smithian CIE, the amount of consumed CO₂ must then range from 6660 Gt C to 115,380 Gt C (Table 1C). Because the amount of volcanogenic CO₂ emitted during a putative middle Smithian eruption is unknown, we use a value of 4.5 × 10¹⁹ g CO₂ for the degassing of SLIP eruption (derived from Sobolev et al., 2011) as a substitute to derive the consumption rate needed to pump this volume of CO₂ during a time interval of 600 kyr (blue cell in Table 1C). We obtain a consumption rate of 6.3 × 10¹² mol C/a, a value that is well within the modern range listed by Dessert et al. (2001). Therefore, the timing of resorption of any putative middle Smithian volcanogenic CO₂ via silicate weathering appears as a plausible process for generating the late Smithian cooling.

Burial of Organic Carbon, the “Biological Pump”

In the following, we assess how much carbon has to be removed from the atmosphere by the biological pump in order to generate a positive CIE with an amplitude of 6‰. Mass balance calculations (Table 1B) yield a drawdown of 25,500 Gt C (Ma) from the atmosphere to be stored in the sedimentary sink. This buried amount is approximately half of the ∼47,000 Gt C that was injected into the atmosphere during the main eruptive event of the SLIP (Sobolev et al., 2011; Me in Table 1A), a value used here for our Smithian case. Hence, it leaves the final atmospheric C reservoir with an excess of 103,500 Gt C (Mf in Table 1B). The 64 ± 104 kyr long hiatus at the top of black shales in Qiakong is equivalent to about 3264 Gt C missing from the budget, which must be added to the 25,500 Gt C estimate for the preserved part of the drawdown (Mb in Table 1B), thus yielding a grand total of ca. 29,000 Gt C for the entire late Smithian. This calculation suggests that the biological pump alone was not able to remove more than ca. 30% of the total atmospheric C and thus was not sufficient for eliciting a late Smithian cooling of 7–8°C on the northern Indian margin (Goudemand et al., 2019).

Implications for Triassic Climate Models

A Timeline for the CIE and Temperature Around the SSB

The calibrated succession of climatic upheavals around the SSB is subdivided into five phases, summarized in Figure 6. The 622 ± 137 kyr long late Smithian positive CIE is best defined in the most complete Qiakong section (Figure 1B). It encompasses the time interval between the end of the middle Smithian negative peak marking the exit from the global thermal maximum (Phase 1), and the positive peak at the base of the Spathian (Phase 5). In Qiakong, the precise timing of the SSB and of Phase 4 are blurred by a gap including the NZ ammonite zone, which is included within the global glacio-eustatic hiatus straddling the boundary. The massive and protracted burial of OM during the late Smithian cooling is considered the cause of a global increase of δ¹³C in marine carbonates. The total shift in δ¹³C between Phase 1 and 5 is of +6‰, thus yielding an average rate of change of approximately +1‰ per 100 kyr. The comparably smaller +4‰ CIE that followed the negative end-Permian peak in deep water and continuous sections (Baresel et al., 2017) had an estimated duration of 30 kyr (Shen and Bowring, 2014) thus resulting in a much higher rate of +13‰ per 100 kyr. As there is a general consensus that the PTB events were linked to the eruptive phase of the SLIP, the ten-fold longer order of magnitude of the late Smithian positive CIE may likely discard a causal relation with a short-term episode of volcanism that would have injected sulfur-rich volatiles into the atmosphere. In the case of the Triassic-Jurassic Boundary (TJB) as a further example, a positive CIE preceded the extinction and was also associated to cooling and a brief phase of glacio-eustatism (Schoene et al., 2010). However, the unknown duration of this positive CIE precludes any calculation of rates.

Reasons for C Cycle Instability

Payne and Kump (2007) explained the unstable global C cycle during the Early Triassic with recurrent injections of volcanogenic (at δ¹³C = -5‰) and thermogenic C from OM sources (at δ¹³C = -25‰), augmented by associated feedback mechanisms (Svensen et al., 2009, 2015). Among the latter, we may identify variations of global temperature and latitudinal temperature gradient, marine anoxia, burial of organic carbon, weathering and riverine C, P input into the oceans, and sea-level oscillations. Exact timing established for the SLIP activity (Burgess et al., 2017; Augland et al., 2019), however, precludes invoking volcanic drivers younger than 250.5 Ma, incompatible with our timeline for an end-Smithian CIE starting at about 249.6 Ma (Figures 4–6 and Supplementary Table S2), and with our middle Smithian ages (samples SHA303 and 309). Weak mercury anomalies of Smithian age in the Boreal realm were also tentatively proposed as evidence for additional volcanism, however, Hammer et al. (2019) showed that these are not of late Smithian but of middle Smithian age and that their volcanic origin is largely uncertain.

Arguments for a Cold Latest Smithian

A late Smithian cooling (starting with the GXZ) in the low latitudes is supported by several independent lines of evidence. Biogeographic distribution of conodonts indicates that cold-water forms (segminiplanate taxa) became unusually abundant within shallow tropical waters (Leu et al., 2019). Following the middle Smithian spore spike of global scope, the ecological recovery of terrestrial plants already started during the late Smithian (Hochuli et al., 2016). Sea surface temperatures measured from pristine conodont apatite from biostratigraphically well-constrained sections of the northern Gondwana shelf (Pakistan) indicate a 7 to 8°C decrease during late Smithian to early Spathian times (Romano et al., 2013; Goudemand et al., 2019). A decrease in the chemical index of alteration at the SSB (Shitouzai section, Nanpanjiang Basin) was documented by Zhang et al. (2015) thus providing another independent argument supporting a cooling episode. Finally, as documented in the present work, reworked Permian zircons occur exclusively within the late Smithian black shales, indicating a lowered base level accompanied by erosion in the hinterland, in agreement with a regressive phase ending by a global glacio-eustatic hiatus (e.g. Embry and Beauchamp, 2008; Davies et al., 2018; Davies and Simmons, 2018).

A Volcanic Origin for the Middle-Early Late Smithian Thermal Maximum?

A recent string of studies investigating geochemical proxies across the end-Smithian extinction (Lyu et al., 2019; Shen et al., 2019; Song et al., 2019; Stebbins et al., 2019a, b; Zhang et al., 2019) all revolve around the scenario of a volcanogenic origin for the middle to early late Smithian thermal maximum. All these works also proposed cessation of volcanism as the primary trigger for the latest Smithian and basal Spathian cooling. Below, we address exclusively key points pertaining to volcanism indicators and timing.

Smithian weak mercury anomalies from a diverse array of sections form a central argument in support of volcanism. These anomalies were causally related to the middle Smithian –early late Smithian thermal maximum by Shen et al. (2019). The modest magnitude of these Smithian anomalies stands in marked contrast with that associated to the PTBME, which is contemporaneous with the main eruptive phase of the Siberian LIP. Unlike the PTBME anomaly, Smithian weak Hg/TOC anomalies do not appear to be synchronous but are scattered throughout middle and late Smithian times. Therefore, clear evidence for substantial middle Smithian volcanism that would lead to a thermal maximum is still wanted (see section “discussion” in Hammer et al., 2019).

Although the upper boundary of the Siberian Traps is an erosional one, no SLIP ages (Burgess et al., 2014, 2017; Burgess and Bowring, 2015; Augland et al., 2019) appears to be younger than 250.60 ± 0.22 Ma an age which is older than those obtained here for the middle Smithian (conodont UAZ2) ashes SHA303 (250.116 ± 0.085 Ma) and SHA339 (249.74 ± 0.12 Ma). The thermal maximum encompasses Phases 1 and 2 and the onset of the cooling is with Phase 3 (Figure 6). Hence, the duration of the thermal maximum is estimated to be in the order of ∼650 kyr (250.1–249.45 Ma). This finding leaves us with two hypotheses: (i) the thermal maximum is younger than the youngest SLIP age, which thus excludes a volcanic origin, and alternatively, (ii) if one admits that the volcanic hypothesis is true, it then implies that the youngest part of the SLIP is missing or has not discovered and dated yet. An argument against the second alternative is provided by the lack of any substantial Smithian mercury anomaly: this suggests that, if there indeed was volcanism, it was volumetrically negligible in comparison to earlier eruptive phases and it is unlikely to explain a ∼650 kyr long thermal maximum.

Expanding on the volcanic hypothesis, Shen et al. (2019) then related the cooling documented across the SSB to the cessation of volcanism. However, younger Hg/TOC anomalies of comparable amplitude during the Spathian were also recorded from the same sections studied by these authors, thus undermining a direct control of Smithian and Spathian climatic changes by volcanism, even if assuming that Smithian, and Spathian weak Hg/TOC anomalies are a reliable proxy for volcanism. Additional volcanic bursts of the SLIP have long been invoked as the usual ad hoc suspect controlling Early Triassic climatic changes (e.g., Ovtcharova et al., 2006) but no clear mercury anomaly nor U-Pb ages support this hypothesis.

The offset between the inception of the positive CIE and the onset of the cooling amounts to 100 to 125 ky in the low latitude (see Figure 6). However, no such data are available for the high latitude. Therefore, it is not yet possible to establish if the cooling was synchronous or not across the entire range of latitude, as the CIE apparently was. As long as the latitudinal pattern of temperature change remains unknown, feedbacks between temperature and the carbon cycle during the Smithian and Spathian cannot be objectively addressed at a global scale.

Conclusion

The high-resolution timeline proposed in this study allows detailed understanding of the sequence of environmental and biological changes during a 2 Myr long period across the SSB. Our model summarized in Figure 6 proposes links between (i) global temperature, (ii) sea-level variations, (iii) variation in terrestrial and marine bio-productivity, and of burial of organic carbon, (iv) weathering intensity, continental runoff and precipitation of inorganic carbon, and eventually (v) evolutionary crises and diversification pulses of marine biota. Silicate weathering alone is shown to be a quantitatively plausible mechanism explaining the late Smithian global cooling. Although the burial of organic carbon (the “biological pump”) must have partly contributed to the cooling, it cannot explain the late Smithian cooling if considered as the sole mechanism. We also point out that a volcanic origin of the immediately preceding thermal maximum still lacks robust evidence, a consequence of which is that volcanic forcing of Early Triassic climate remains an ad hoc and weak explanation.

Data Availability Statement

All data are contained in Table 1 and Supplementary Tables S1–S5. Background data, metadata, CL images of dated zircon crystals, and further analytical protocols not contained in the methods section are available on request from the authors.

Author Contributions

HB and US designed the research. HB, ML, BB, and NG collected the field data and the samples. PW carried out the U-Pb geochronology work. TV, BB, ML, ES-H, and NG collected the other data. All authors contributed the data interpretation. PW, HB, and US wrote the manuscript with the help of all other authors.

Funding

This study was supported through the Swiss National Science Foundation grant numbers 156424 and 182007 (US), and 160055 (HB).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

All authors thank the Geneva and Zürich groups for technical support, especially Maria Ovtcharova for invaluable support in U/Pb analytics and Thierry Adatte for RockEval analyses. Ji Cheng and Kuang Guodun are deeply thanked for their long-term help in the field. The final version of this work benefited from constructive and helpful reviews by TA and SB, as well as by the Frontiers editor MC.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2020.00196/full#supplementary-material

References