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

The dynamics of the carbon cycle across different timescales is crucial for understanding past and present global climate changes. Following the Permian–Triassic boundary mass extinction (PTBME), the carbon cycle changed profoundly during the following 5.4 Myr, with magnitudes of changes comparable to those of the Precambrian. In pace with the successive cycles of the carbon budget, the recovery of the marine nekton underwent several evolutionary diversification and extinction cycles accompanied by eustatic sea-level changes and profound ecological reorganization of land plants, all indicative of climatic changes. Additional eruptive bursts of the Siberian Large Igneous Province (SLIP) are traditionally called upon as a plausible trigger for these climatic oscillations but firm evidence for coeval SLIP volcanism is still lacking. Based on new precise and accurate U-Pb zircon ages, we establish a high-resolution temporal calibration of the biggest positive carbon isotope excursion (CIE) spanning about 600 kyr in the late Smithian. The age of the Smithian-Spathian boundary (SSB) is established between 249.29 ± 0.06 and 249.11 ± 0.09 Ma. Our oldest U-Pb zircon ages indicate no overlap in time between the middle Smithian onset of the thermal maximum and the youngest available U-Pb zircon ages from SLIP volcanism. The constructed time line also indicates a duration of the global unconformity at the SSB that is compatible with glacio-eustatism. Potential cooling mechanisms such as a volcanic winter, the biological pump and the cessation of volcanism are discussed in the light of this new time line. In the low latitudes, the onset of the positive CIE remarkably predates the temperature drop by some 100 to 125 kyr. However, as long as the magnitude of such offset – if any – is unknown for the high latitudes, relations between the CIE and the cooling will remain an open question associated with largest Triassic extinction of the nekton.


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 2 drawdown by silicate weathering or by the biological pump in the hypothetical frame of the cessation of a middle Smithian SLIP volcanic episode.

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) 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 3 in 3 ml Savillex beakers before the partial dissolution for 18 h at180 • C in 40% HF and trace HNO 3 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 3 before fluxing overnight at 80 • C on a hotplate in 6 N HCl. Zircon were washed again several times in 3 N HNO 3 . Single zircon crystals were loaded into clean 200 µl microcapsules, spiked with ∼4-6 mg of the EARTHTIME 202 Pb-205 Pb-233 U-235 U tracer solution McLean et al., 2015), and dissolved in ∼70 µl 40% HF and trace HNO 3 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 3 PO 4 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 202 Pb/ 205 Pb ratio. All common Pb in the analysis was assumed to be procedural blank ( 206 Pb/ 204 Pb 17.6186 ± 0.3677, 207 Pb/ 204 Pb 14.7298 ± 0.4504, and 208 Pb/ 204 Pb 35.77 ± 1.07; 2σ). Uranium was measured as U-oxide in static mode on Faraday cups equipped with 10 12 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 (δ 13 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).
resistors or in dynamic mode on a MasCom secondary electron multiplier. Isotopic ratios were corrected for interferences of 233 (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 206 Pb/ 238 U dates were corrected for initial 230 Th/ 238 U disequilibrium using a Th/U of 3.5 ± 1 (2σ) for magma composition, augmenting their 206 Pb/ 238 U age by ∼80-100 kyr. All reported ages are weighted mean 206 Pb/ 238 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 202 Pb-205 Pb-233 U-235 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 206 Pb/ 238 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 agedepth 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 (δ 13 C carb , δ 18 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 (δ 13 C carb and δ 18 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. Inhouse 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 δ 13 C and δ 18 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 δ 13 C and ± 0.1 for δ 18 O values. The results are expressed in the δ-notation [δ = (R1/R2-1) × 1000] where R1 is the 13 C/ 12 C or 18 O/ 16 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 Frontiers in Earth Science | www.frontiersin.org 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 2 /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.

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, 1991Guex, , 2011Hammer, 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 buildin 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 . 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 δ 13 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 δ 13 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 nonobjective 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).

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 206 Pb/ 238 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 206 Pb/ 238 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 206 Pb/ 238 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 206 Pb/ 238 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 206 Pb/ 238 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 206 Pb/ 238 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 206 Pb/ 238 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 206 Pb/ 238 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 206 Pb/ 238 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 206 Pb/ 238 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 longprismatic zircon grains was determined, seven of them yielding a best estimate 206 Pb/ 238 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 206 Pb/ 238 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 206 Pb/ 238 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 206 Pb/ 238 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 FIGURE 3 | 206 Pb/ 238 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 238 U decay constant uncertainty).
Frontiers in Earth Science | www.frontiersin.org (ca. 6.3 m above the nodular limestone) were analyzed zircon are exceptionally low in radiogenic lead. The best estimate 206 Pb/ 238 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 ( 206 Pb/ 238 U = 249.110 ± 0.088/0.11/0.29 Ma 2σ; MSWD = 1.8, n = 7) and QIA07T ( 206 Pb/ 238 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 ( 206 Pb/ 238 U = 246.40 ± 0.018 Ma) and QIA07T ( 206 Pb/ 238 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 206 Pb/ 238 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 206 Pb/ 238 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.
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 δ 13 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  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(Black et al., , 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(Miller et al., , 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 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 206 Pb/ 238 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.
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.
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 2 originated from a volumetrically important middle Smithian eruptive burst of the SLIP.

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 2 sink during times of warm and humid climate. We adopt present-day CO 2 consumption rates of basaltic provinces from Dessert et al. (2003) and a minimal surface area (2.5 × 10 6 km 2 ) 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 2 must then range from 6660 Gt C to 115,380 Gt C (Table 1C). Because the amount of volcanogenic CO 2 emitted during a putative middle Smithian eruption is unknown, we use a value of 4.5 × 10 19 g CO 2 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 2 during a time interval of 600 kyr (blue cell in Table 1C). We obtain a consumption rate of 6.3 × 10 12 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 2 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 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. 206 Pb/ 238 U single-grain zircon ages of dated ash layers are indicated, for sample numbers see Figure 1. 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 .

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 δ 13 C in marine carbonates. The total shift in δ 13 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  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.
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 δ 18 O phosN  and associated consequences and feedbacks. Payne and Kump (2007) explained the unstable global C cycle during the Early Triassic with recurrent injections of volcanogenic (at δ 13 C = -5 ) and thermogenic C from OM sources (at δ 13 C = -25 ), augmented by associated feedback mechanisms (Svensen et al., 2009(Svensen et al., , 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 sealevel 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 . 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 glacioeustatic hiatus (e.g. Embry and Beauchamp, 2008;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  Dessert et al. (2003), *estimated numbers to match the Me 47,000 Gt of (1).
Frontiers in Earth Science | www.frontiersin.org 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(Burgess et al., , 2017Burgess 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). U/Pb analytics and Thierry Adatte for RockEval analyses. Ji Cheng and Kuang Guodun are deeply thanked for their longterm 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.