Latest Permian carbonate carbon isotope variability traces heterogeneous organic carbon accumulation and authigenic carbonate formation

. Bulk-carbonate carbon isotope ratios are a widely applied proxy for investigating the ancient biogeochemical carbon cycle. Temporal carbon isotope trends serve as a prime stratigraphic tool, with the inherent assumption that bulk micritic carbonate rock is a faithful geochemical recorder of the isotopic composition of seawater dissolved inorganic carbon. However, bulk-carbonate rock is also prone to incorporate diagenetic signals. The aim of the present study is to disentangle primary trends from diage-netic signals in carbon isotope records which traverse the Permian–Triassic boundary in the marine carbonate-bearing sequences of Iran and South China. By pooling newly produced and published carbon isotope data, we conﬁrm that a global ﬁrst-order trend towards depleted values exists. However, a large amount of scatter is superimposed on this geo-chemical record. In addition, we observe a temporal trend in the amplitude of this residual δ 13 C variability, which is reproducible for the two studied regions. We suggest that (sub-)sea-ﬂoor microbial communities and their control on calcite nucleation and ambient porewater dissolved inorganic carbon δ 13 C pose a viable mechanism to induce bulk-rock δ 13 C variability. Numerical model calculations highlight that early diagenetic carbonate rock stabilization and linked carbon isotope alteration can be controlled by organic matter supply and subsequent microbial remineralization. A major biotic decline among Late Permian bottom-dwelling organisms facilitated a spatial increase in heterogeneous organic carbon accumulation. Combined with low marine sulfate, this resulted in varying degrees of carbon isotope overprint-ing. A simulated time series suggests that a 50 % increase in the spatial scatter of organic carbon relative to the average, in addition to an imposed increase in the likelihood of sampling cements formed by microbial calcite nucleation to 1 out of 10 samples, is sufﬁcient to induce the observed signal of carbon isotope variability. These ﬁndings put constraints on the application of Permian–Triassic carbon isotope chemostratigra-phy based on whole-rock samples, which appears less reﬁned than classical biozonation dating schemes. On the other


Introduction
Carbon isotopes in carbonate rock are a pivotal tool for understanding the ancient biogeochemical carbon cycle and a stratigraphic aid in determining the age of sedimentary deposits.Individual fossil carbonate shells are the preferred recorder of the isotope composition of marine dissolved inorganic carbon (DIC).However, some deposits, such as those of Precambrian age or those which were formed during biotic crises or where shelly fossils are absent, are not suitable for such a single-component approach (Veizer et al., 1999;Prokoph et al., 2008).Under these circumstances, bulk-rock samples are a widely used alternative recorder (Saltzman, 2001;Brand et al., 2012a, b).This approach has attracted criticism due to the complex multicomponent nature and high potential for diagenetic alteration, which may result in a mixture of primary and secondary diagenetic signals.Major sources of contamination that are known to affect the primary marine carbon isotope signal of bulk-carbonate rock are degradation of organic matter and methane oxidation (Marshall, 1992;Reitner et al., 2005;Wacey et al., 2007;Birgel et al., 2015).
Carbonate sediments are especially prone to chemical alteration during early compaction and lithification when the highly porous sediment and unstable carbonate polymorphs (e.g., aragonite and high-Mg calcite) are subjected to cementation and dissolution (Bathurst, 1975;Irwin et al., 1977;Marshall, 1992).The potential for diagenetic alteration of biogenic high-Mg calcite and aragonite results from elevated substitution of foreign ions in the lattice and a higher degree of dislocations (Busenberg and Plummer, 1989).
On carbonate platforms subjected to rapid sea level fluctuations, lithification might be controlled by interaction with meteoric fluids and oxidized terrestrial organic matter (Brand andVeizer, 1980, 1981;Banner and Hanson, 1990;Knauth and Kennedy, 2009).The latter mode of carbonate rock cementation with allochthonous dissolved carbonate sources (meteoric water) in a relatively open system has long been viewed as the predominant process of carbonate rock stabilization (Bathurst, 1975).These notions were fuelled by an over-representative focus on Pleistocene carbonate platforms (e.g., the Bahamas) and their associated diagenetic stabilization under the influence of glacial-interglacial-induced sea level fluctuations (James and Choquette, 1983).However, this mechanism is less relevant for rock lithified under greenhouse conditions when sea level changes are dampened (Bathurst, 1993;Munnecke and Samtleben, 1996;Melim et al., 2002).
Alternatively, carbonate rock alteration is driven by marine-derived fluids that evolved through interaction with precursor sediments and ambient chemical conditions (Munnecke and Samtleben, 1996;Melim et al., 2002).Foremost, the reactive and biologically active upper section of the sediment column can be envisioned as leaving an imprint on carbonate chemistry, including the carbon isotope composition.This carbon isotope signal will be retained in postdepositional carbonate cements and is controlled by organic carbon (OC) fluxes, pH, alkalinity and the nature of the (sub-)sea-floor microbial communities (Melim et al., 2002;Reitner et al., 2005;Wacey et al., 2007;Birgel et al., 2015;Schobben et al., 2016;Zhao et al., 2016).
A more nuanced view on the nature of bulk-rock carbon isotope composition would be to envision diagenetic overprinting on carbonate rock chemistry as a spectrum where pure primary and strictly diagenetic end-member states are considered as a continuum, with bulk-rock chemistry falling somewhere in between these two extremes (Marshall, 1992).In this case, the aim should be to discern the degree to which an original imprint could have been retained in the bulk-rock signal.Embracing both ends of this spectrum could illuminate new aspects of marine sedimentary systems and ancient ocean chemistry, as well as providing a more refined understanding of the mechanisms that promote the diagenetic alteration of carbonate sediments.Clear examples of the importance of interpreting diagenetic signals have been given by the assertion that authigenic carbonate production might play a primary role in the global biogeochemical carbon cycle (Schrag et al., 2013) and in the isotopic imprint of terrestrial biomass on Precambrian carbonates (Knauth and Kennedy, 2009).

Background: Permian-Triassic carbonate carbon isotope records
Temporal trends in carbonate carbon isotope composition have been recorded at varying timescales from million-year secular trends (Veizer et al., 1999;Zachos et al., 2001) to dramatic sub-million-year events, such as at the Paleocene-Eocene Thermal Maximum (Dickens et al., 1995).The sedimentary record of the Permian-Triassic (P-Tr) boundary interval, marked by the largest mass extinction of the Phanerozoic (Erwin, 1993;Alroy et al., 2008), is also characterized by pronounced carbon isotope excursions, which are almost exclusively recorded in bulk rock.However, chemical signatures from this time period include records with a wide variety of amplitudes and shapes and span centimetres to metres of rock sequence (e.g., Xie et al., 2007;Korte et al., 2010).Equally diverse are the proposed drivers of these geochemical signals, ranging from prolonged episodes of volcanism and the eustasy-controlled erosion of organic-rich shelf sediments to abrupt blooms of methane-producing microbes (Renne et al., 1995;Berner, 2002;Rothman et al., 2014).
Variable-amplitude carbon isotope excursions in shelf to basin transects have been linked to an intensified biological pump (Meyer et al., 2011;Song et al., 2013).This notion contrasts with the viewpoint of collapsed primary productivity (or a "Strangelove" ocean) as a driver for the demise of bottom-dwelling and infaunal marine biota (Rampino and Caldeira, 2005).The loss of these geobiological agents resulted in a reduced thickness of the sedimentary mixed layer and a global increase in the occurrence of laminated sediments (Wignall and Twitchett, 1996;Hofmann et al., 2015).Contrary to disrupted primary productivity, these geological features would also agree with a scenario of enhanced OC remineralization and resulting widespread marine deoxygenation (Meyer et al., 2008(Meyer et al., , 2011)).These contrasting scenarios of global-scale environmental deterioration are difficult to reconcile.Nevertheless, the reduction of sediment mixing is an unequivocal feature that warrants consideration when interpreting P-Tr chemical records.
An elevated oceanic carbonate inventory is a predicted aspect of an ocean without pelagic calcifiers, which only started proliferating during the Mesozoic and would effectively buffer short-term perturbations to P-Tr marine carbon isotopes (Kump and Arthur, 1999;Zeebe and Westbroek, 2003;Rampino and Caldeira, 2005;Ridgwell, 2005;Payne et al., 2010).This does not, however, exclude local departures from this dynamic equilibrium or depth-related isotope differences, such as those forced by the biological pump (Kump and Arthur, 1999;Rampino and Caldeira, 2005).High carbonate ion concentrations are invoked to explain the advent of (microbial) carbonate sea-floor structures (e.g., thrombolites, stromatolites and fan-shaped structures) in the extinction aftermath and the ensuing Early Triassic recovery phase (Baud et al., 1997;Rampino and Caldeira, 2005;Riding and Liang, 2005;Pruss et al., 2006;Kershaw et al., 2007Kershaw et al., , 2009;;Leda et al., 2014).While conditions favoured the formation of these structures, poorly buffered calcified metazoans (e.g., brachiopods and corals) were proportionally more affected by the end-Permian mass extinction (Knoll et al., 2007).This biotic shift has been interpreted as a prevalent carbonate factory turnover from skeletal to microbial (e.g., Kershaw et al., 2009).
A carbonate factory turnover, reduced sediment mixing and increased OC sinking fluxes invite consideration of a scenario in which these individual parameters act as synergistic effects on carbonate formation and diagenetic stabilization.Moreover, the imprint of 12 C-depleted carbon on latest Permian carbonate rock produced by the introduction of authigenic carbonate has been connected to a systematic carbon isotope offset between bulk rock and brachiopod shells (Schobben et al., 2016).Diagenetic and authigenic carbonate sources can, however, result in a range of carbon isotope values relating to the specific microbial community and sedimentary environment (Irwin et al., 1977).Hence, we hypothesize that spatial carbon isotope variability at the P-Tr boundary interval relates to an increased importance of mi-crobially controlled calcite nucleation.Moreover, we postulate that organic carbon sinking fluxes and subsequent in situ remineralization by microbes determine the trajectory of carbonate rock stabilization.Combined geographic differences in the OC sinking fluxes (km scale) and sediment mixing (cm scale) might have generated spatially heterogeneous dispersion of organic carbon.This spatial pattern of OC distribution links to observed carbonate carbon isotope variability by modulating the chance of sampling variable isotope signals along lateral lithological transects.We adopt this conceptual model as a working hypothesis for interpreting stratigraphic carbon isotope patterns in P-Tr carbonate rock.In addition, however, this scenario likely has implications for the interpretation of bulk-rock carbon isotope patterns during other periods of the Precambrian and Phanerozoic.
To carry out this investigation, we studied carbonatebearing sequences of Permian to Triassic strata located in Iran and China using a compilation of published and new data.The effect of microbial metabolism on sediments and porewater can be numerically approximated by reactive transport models (Soetaert et al., 1996;Boudreau, 1997;Meysman et al., 2003;van de Velde and Meysman, 2016).Similarly, diagenetic models have proven useful in delineating trajectories of bulk-rock Sr and Ca isotope stabilization (Fantle andDePaolo, 2006, 2007).For the purpose of our study, we will combine aspects of these models in an effort to estimate the potential for the microbially mediated spatial modification of bulk-rock carbon isotope signals.

Material selection and stable carbon isotope measurements
The P-Tr limestone sequences of Iran are particularly suited to study spatial δ 13 C carb variations, as lateral lithological homogeneity excludes a strong control from palaeoceanographic conditions or selective preservation potential (see the Supplement).In addition to the classical P-Tr sites of Shahreza, Zal and Ali Bashi in Iran, we sampled several other sites for which we present the first δ 13 C carb results.These sites include the Aras Valley profile (39.015A total of 463 stable carbon isotope measurements were made at the University of Copenhagen (UC) and the Museum für Naturkunde (MfN), Berlin.Glass reaction vessels (Labco ® ) containing the sample powders were flushed with helium, and carbonate was left to react with 50 (UC) or 30 (MfN) µL of anhydrous phosphoric acid (∼ 102 %) for at least 1.5 h.Carbon and oxygen isotope values were measured from resulting CO 2 using an Isoprime triple collector isotope ratio mass spectrometer in continuous-flow setup (UC) or a Thermo Finnigan Gasbench (GS) II linked to a Thermo Finnigan MAT V isotope ratio mass spectrometer (IRMS; MfN).Pure CO 2 (99.995 %) calibrated against international IAEA standards (NSB-18 and NSB-19) was used as a reference gas.At the UC, results were corrected for weightdependent isotope ratio bias using multiple measurements of an in-house standard LEO (Carrara marble; δ 13 C = 1.96 ‰) covering the entire range of signal intensities encountered in the samples.External reproducibility was monitored through the replicate analysis of the in-house standards LEO and Pfeil STD (Solnhofen limestone) at the UC and MfN, respectively.Long-term accuracy was better than 0.1 ‰ (2 SD; better than 0.2 ‰ 2 SD for oxygen isotopes).All carbon isotope values are reported in ‰ relative to VPDB and in standard δ notation.

Complementary data collection
Our new data are complemented by published bulk-carbonate δ 13 C data from multiple P-Tr localities in Central Iran (Abadeh (Heydari et al., 2000;Korte et al., 2004a;Horacek et al., 2007;Korte et al., 2010;Liu et al., 2013), Shahreza (Korte et al., 2004b;Heydari et al., 2008;Richoz et al., 2010) and north-western Iran (Ali Bashi mountains and Zal; Baud et al., 1989;Korte et al., 2004c;Korte and Kozur, 2005;Kakuwa and Matsumoto, 2006;Horacek et al., 2007;Richoz et al., 2010;Schobben et al., 2016).For a more global perspective and a comparative approach, we additionally extracted published data from the P-Tr Global Boundary Stratotype Section and Point (GSSP) at Meishan, China (Baud et al., 1989;Chen et al., 1991;Xu and Yan, 1993;Hoffman et al., 1998;Jin et al., 2000;Cao et al., 2002;Gruszczyński et al., 2003;Zuo et al., 2006;Kaiho et al., 2006;Huang et al., 2007;Riccardi et al., 2007;Xie et al., 2007).Data were extracted from tables and supplementary files or provided by the authors and where necessary read from figures with the open-source software xyscan (Ullrich, 2016).This compilation adds 2077 data points to our new dataset.The analytical uncertainty of the collected data, when given, ranges between 0.02 and 0.2 ‰.We included replicate studies on the same site, as these data are potentially valuable assets to test the reproducibility of stable isotope investigations and to shed light on the effect of the bulk-rock multicomponent nature on δ 13 C carb composition.

Data projection
To facilitate a direct comparison of δ 13 C carb between different sites with differences in total sediment thickness, we converted the stratigraphic height to a dimensionless timeline.This approach preserves a more direct connection to the sampled rock sequence rather than converting to absolute ages or maintaining original stratigraphic heights.The conversion to a dimensionless time grid was carried out by using the lower and upper bounds of individual conodont assemblage biozones as tie points to which a relative distance for an individual δ 13 C carb data point is assigned (see the Supplement).The chronological scheme used here divides the P-Tr interval into biostratigraphic units (A-K), which enables the correlation of individual sequences in both geographic regions (Table 1).
To evaluate trends in the collected data, a sliding window with a bandwidth equivalent to the dimensionless grid unit has been applied to the δ 13 C carb data.This method ensures the extraction of temporal trends at equivalent time resolution to the biozonation units (Fig. 1).
To cancel out the potential for an uneven spread of data points, we also applied a subsampling routine on the Iranian and Chinese datasets.Subsampling was performed with the sliding window procedure, as described above, before applying summary statistics.Subsequently, the median values were calculated for each of the subsampled δ 13 C carb populations.For simplification, the resulting median trend line and 95 % confidence intervals (CIs) are calculated and visually weighted (Fig. 2) (visually weighted plot or watercolour regression; Hsiang, 2012;Schönbrodt, 2012).The latter is accomplished by calculating a kernel density, and both the drawn median trend line and CI interval are weighted by colour saturation based on this kernel density value.In addition, the median trend line and the CI interval have been assigned contrasting colours to further enforce the visual signal.Hence, saturated and contrasting parts of the depicted regression curve (i.e., sections with high visual weight) depict areas with the highest fidelity of the temporal geochemical pattern based on the subsampling routine.In addition to the median trend lines, the same subsampling approach and visual weighting has been applied to calculate and plot the interquartile range (IQR; the value range containing 50 % of the sample population) and 95 % interpercentile range (IPR; the value range containing 95 % of the sample population).Graphs are plotted with the open-source programming platform R (R Core Team, 2016) and with the aid of the R packages; ggplot2 (Wickham, 2009), reshape2 (Wickham, 2007), plyr (Wickham, 2011) and gridExtra (Auguie, 2016).3.2 Reactive transport modelling

Geochemical model formulation and biogeochemical reactions
Organic matter availability fuels in situ metabolic pathways linked to calcite nucleation.Upon entering the sediment, organic matter is microbially mineralized, which is coupled to the reduction of electron acceptors.These electron acceptors are consumed in a well-defined sequence based on their thermodynamic energy yield: O 2 , SO 2− 4 and methanogenesis (Berner, 1964;Froelich et al., 1979).Each of these microbial metabolisms will imprint a specific carbon isotope signature on the porewater DIC and thus create a potential source for diagenetic alteration of carbonate C isotope signatures (Irwin et al., 1977).
These mineralization processes can be described by mass balance equations (Eqs. 1 and 2), which can subsequently be solved numerically via the method of lines (Boudreau, 1996;Soetaert and Meysman, 2012).The numerical solutions of these equations solve the age-depth (t and z) relationship of deposited sediments in terms of solids (S i ), solutes (C i ) in pore fluids and their respective reactions between the solid and liquid phase.At the top of the sediment pile the porewaters communicate with ocean water so that dissolved elements can diffuse according to their concentration gradients.In addition to diffusion transport processes, continued sedimentation supplies the sediment with organic carbon and calcium carbonate.
In Eqs.
(1) and (2), ϕ is porosity, D i is the effective diffusion coefficient and D b and α(z) are parameters associated with bioturbation (see below); ν i,k is the stoichiometric coefficient of species C i in reaction R k .Note that we express reactions as mol per unit time per volume sediment and concentrations as mol per volume porewater or volume solid phase.Therefore ϕ and 1-ϕ are introduced as unit conversions.The model includes three different modes of transport: sedimentation (represented by downward advection of solutes v and solids w), molecular diffusion and biological transport (bioturbation).Since we only consider cohesive sediments, the only advective transport is burial; i.e., new sediment is added on top of the sediment column, and sediment at the bottom of the column is buried (transported out of the column).Molecular diffusion for porewater solutes is expressed by Fick's first law, where D i (the effective diffusion coefficient) is calculated following the definition as given in Boudreau and Meysman (2006).

Biozones
Meishan (Chen et al., 1991) Meishan (Hoffman et al., 1998) Meishan (Xu & Yan 1993) Meishan (Zuo et al., 2006) Meishan A (Kaiho et al., 2006) Meishan B (Huang et al., 2007) Meishan B (Jin et al., 2000) Meishan B (Xie et al., 2007) Meishan D (Baud et al., 1989) Meishan D (Cao et al., 2002) Meishan The individual carbonate carbon isotope values are placed on a dimensionless timeline that marries both geographic areas in the most acceptable biochronological scheme (Table 1).The solid black line represents the subsampled median trend line.The dashed blue line depicts the seawater δ 13 C DIC curve as obtained from the time series simulation (Sect.5).The stratigraphic placement of the sea level changes and the extinction horizon as well as the biostratigraphic framework can be found in the Supplement.
where D 0 is a function of temperature and salinity and has been calculated with the R package marelac (Soetaert et al., 2016).Bioturbation is implemented as two separate processes: bio-mixing and bio-irrigation (Kristensen et al., 2012) following conventional descriptions (Boudreau, 1984;Meysman et al., 2010).The bio-mixing is modulated over the depth interval to account for the effects of sediment reworking by metazoans in the top layer of the sediment pile.In this formulation it is assumed that D b remains constant in a layer with thickness z b , after which it attenuates with depth by following an exponential relation with coefficient λ D b and background bio-diffusivity (D b0 ) at the seawater-sediment interface (Soetaert et al., 1996).
Bio-irrigation exchanges porewater with the overlying water via burrow flushing.This is implemented via a non-local exchange process (Boudreau, 1984).
The quantity α(z) represents the depth-dependent irrigation intensity, and the solute concentrations of the bottom water and at depth are given by C ow and C(z), respectively.The attenuation of bio-irrigation is expressed by an exponential relation.
In this formulation, α 0 is the irrigation coefficient at the sediment-water interface, and x irr is the attenuation depth coefficient.Most of the macrofaunal activity takes place in the upper few centimetres of the sediment (as animals are dependent on food resources that rain down via the water column).Therefore, both bio-mixing and bio-irrigation processes are most intense near the sediment-water interface.
The aim of this model is to give a parsimonious description of the potential effect of early diagenetic reactions on the isotope signature of carbonates.Therefore, we do not consider nitrification or metal cycling, as this would increase the complexity of the model and calculations; in addition, this would not fundamentally alter our conclusions.Organic mat- Secondary redox reactions Aerobic oxidation of methane Carbonate precipitation and dissolution

R carb
Metabolism-produced carbon isotope values according to (Irwin et al., 1977), except organoclastic methanogenesis-sourced CH 4 subsequently employed by AOM and ArOM, which is a conservative estimate between −65 ‰ and near-quantitative sedimentary OC conversion to CH 4 , thereby approaching the parent OC-C isotope composition of −25 ‰.
ter can be mineralized via different pathways: aerobic respiration, sulfate reduction and organoclastic methanogenesis (in which organic matter acts both as electron donor and acceptor; Table 2).The reduction of sulfate produces sulfide, which can be oxidized by oxygen via canonical sulfide oxidation (CSO) (Meysman et al., 2003;van de Velde and Meysman, 2016).Methanogenesis produces methane, which can be oxidized by sulfate in a process called anaerobic oxidation of methane (AOM) or by oxygen (ArOM) (Boetius et al., 2000).The mineralization reactions are expressed via standard limitation-inhibition formulations (Soetaert et al., 1996).The Monod constants (K O 2 and K SO 2− 4 ) in these formulations determine the inhibition or limitation of a certain electron acceptor.For example, at high concentrations of oxygen, oxygen reduction will be efficient and other mineralization pathways will be inhibited.At lower concentrations, oxygen reduction will be less efficient, and other pathways will gain importance (Soetaert et al., 1996;Van Cappellen and Wang, 1996;Berg et al., 2003;Meysman et al., 2003).This sequential alternation of limitation and subsequent inhibition of metabolic pathways results in chemical zonation of the sediment profile (Froelich et al., 1979).The reoxidation reactions are given as second-order rate laws in which the kinetic constants (k ArOM , k AOM and k CSO ) determine the reaction rate of the respective reactions (Table 2; Boudreau, 1997;Meysman et al., 2003).
To mimic calcium carbonate recrystallization following an age-depth relationship, dissolution and carbonate production reaction rates, expressed as Eq. ( 7), have been incorporated in the model with no net carbonate dissolution or precipitation (Fantle andDePaolo, 2006, 2007).
In contrast to the models of Fantle andDePaolo (2006, 2007) we envision that only a fraction of the whole carbonate rock determines the reactivity towards dissolution and recrystallization.This fraction of 0.2 is considered to be the diagenetic carbonate fraction (f dia ) of the total rock volume (e.g., Schobben et al., 2016).Specific metabolically produced DIC δ 13 C values are assigned for organoclastic microbial sulfate reduction, AOM and organoclastic methanogenesis (see Table 2).These metabolic pathways have been considered to induce calcite nucleation by changing the surrounding carbonate chemistry, forming an organic matrix or destroying calcification-inhibiting organic mucus (Berner et al., 1970;Visscher et al., 2000;Arp et al., 2001;Ries et al., 2008;Heindel et al., 2013;Birgel et al., 2015).The peak activity of these metabolic pathways coincides with the most reactive upper layers of the sediment column, i.e., highly porous and high metastable calcite polymorph content (Irwin et al., 1977;Marshall, 1992); this is mathematically expressed by Eq. ( 7).Solute diffusion tracts and the sedimentation of the carbon isotopes of DIC and solid-phase carbonate are defined as 13 C and 12 C, respectively, and are calculated separately.Metabolically introduced carbon (v metabolism ) by reaction (R k ) and its respective effect on porewater DIC carbon isotope composition is formulated by Eqs. ( 8) and ( 9) for the heavy and light isotope, respectively.
We assign a 13 C / 12 C ratio (r) to the metabolism-specific produced DIC from the organic carbon remineralization reactions according to the data in Table 2.The mathematical formulation of metabolism-introduced DIC ensures that the isotope composition of ambient porewater DIC is altered and ensures mass balance in the model carbon input and output.

Numerical solutions and parameters
The model includes nine state variables: organic matter (CH 2 O), solid carbonate ( 13 C carb and 12 C carb ), DIC ( 13 C DIC and 12 C DIC ), dissolved oxygen (O 2 ), dissolved sulfate (SO2− 4 ), dissolved methane (CH 4 ) and dissolved sulfide (HS − ).We numerically solved Eqs. ( 1) and ( 2) on the opensource programming platform R (R Core Team, 2016) with a finite difference approach by expanding the spatial derivatives of partial differential equations over a sediment grid (Boudreau, 1997).Through the application of the R package ReacTran (Soetaert and Meysman, 2012), a sediment grid 10 m thick with 20 000 layers of identical thicknesses was generated.The finite difference approach transforms the partial differentials into ordinary differentials which are subsequently integrated by a stiff solver routine vode (Brown and Hindmarsh, 1989) within the R package DeSolve (Soetaert et al., 2010).Boundary conditions were chosen to reflect Permian oceanic conditions based on the current state of research (Sect.2) and complemented by current knowledge of shallow marine environments.The upper boundary conditions were set to concentrations for normal modern bottom water conditions: 0.28 µmol cm −3 O 2 and zero for the reduced species (CH 4 and HS − ).A lower seawater SO 2− 4 value of 4 µmol cm −3 (compared to the modern value of 28 µmol cm −3 ) was chosen based on evidence of a reduced global marine S reservoir in the latest Permian to Early Triassic (Luo et al., 2010;Song et al., 2014;Schobben et al., 2017).On the other hand, DIC was set at 4.5 µmol cm −3 (higher than the modern value of 2.2 µmol cm −3 ) based on model calculations of a Late Permian ocean without pelagic calcifiers (Ridgwell, 2005).The solid phases entering the model at the top of the sediment stack are set at 730 µmol cm −2 yr −1 for both the OC flux (F OC ) and the carbonate flux (F carb ) in the baseline model, which are typical average shelf sedimentation values (Müller and Suess, 1979;Reimers and Suess, 1983;Sarmiento and Gruber, 2004).For the baseline condition, sedimentation rates (v and w) are fixed at 0.2 cm yr −1 and the δ 13 C composition is set at +5 ‰ (VPDB) based on primary carbon isotope values from pristine preserved brachiopod calcite from the Ali Bashi and Meishan sections, as well as sites in northern Italy (Korte et al., 2005;Brand et al., 2012b;Schobben et al., 2014).For the bio-mixing and bio-irrigation parameters, the values for the Palaeozoic proposed by Dale et al. (2016) were chosen: D b0 = 5 cm 2 yr −1 , z b = 2 cm, α 0 = 50 yr −1 , x irr = 1 cm.This will approximate the sedimentary conditions of a sea floor inhabited by a Permian benthic community.The relation between the isotope composition of DIC and solid-phase carbonate is given by the following thermodynamic relation (Eq.10) established by Emrich et al. (1970) in which we assume temperature (T ) to be constant at 25 The lower boundary conditions for both solid and solute species were set at a "no gradient" boundary to ensure that materials are only transported by burial to deeper parts of the sediment column.The kinetic parameters are summarized in Table 3 and are based on published models.

Bulk-carbonate carbon isotope signal imprinting and its sensitivity towards environmental parameters
The previously defined parameters provide a baseline model to test the sensitivity towards certain boundary conditions in bringing forward carbon isotope compositional changes in the porewaters.Two pathways of carbonate rock formation and stabilization are considered to be of importance in capturing these porewater carbon isotope compositional changes.
1.The first is diagenetic carbonate alteration, which is envisioned as ongoing recrystallization and dissolution with depth after burial, as formulated in Eq. ( 7), and with both processes operating in equilibrium (Fantle andDePaolo, 2006, 2007).Under this mechanism carbonate is constantly exchanged between the solid (carbonate precipitate) and aqueous phase (porewater DIC) with ongoing burial, thereby to some degree buffering the isotopic perturbations caused by microbial mineralization.This pathway of carbonate sediment progression can be viewed as the classical interpretation of limestone stabilization: dissolution of less stable components and subsequent occlusion of the produced pore spaces by cements (Bathurst, 1993;Munnecke and Samtleben, 1996;Melim et al., 2002).a time of unusual carbonate sedimentation; e.g., stromatolites, thrombolites and microbially induced cementation (Baud et al., 1997;Woods et al., 1999;Leda et al., 2014).However, the mechanism directly responsible for this type of biologically induced calcification is ambiguous and might be related to extracellular polymeric substances and the orientation of functional groups on this organic matrix (creating nucleation sites), destruction of organic compounds that inhibit calcification or a direct consequence of metabolism and associated changes in ambient carbonate chemistry (Visscher et al., 2000;Arp et al., 2001;Bundeleva et al., 2014).We further assume that cryptic forms of carbonate precipitates, e.g., rock-binding cements (e.g., Richoz et al., 2010), might have escaped detection when selecting carbonate rock for carbon isotope studies.
These two pathways are henceforward referred to as diagenetic carbonate alteration and authigenic carbonate addition, respectively.Note that porewater chemistry and its control on carbonate formation are not explicitly formulated in the model and are instead expressed by a reaction rate that can vary with age (Eq.7) (Fantle andDePaolo, 2006, 2007).This formulation is justifiable due to the uncertainties in the ultimate control of biological activity in steering calcification, as cited above.Diagenetic carbonate alteration is determined by taking the evolved solid-phase δ 13 C at 10 m of sediment depth where carbon isotope exchange between porewater and solid carbonate reaches equilibrium (see also Fantle andDePaolo, 2006, 2007, for a comparative Sr and Ca isotope equilibrium with depth).Although transient seawater chemistry changes can occur over the time needed to complete carbonate stabilization (i.e., before reaching equilibrium at 10 m of depth), we consider the adopted isotope value to be a good approximation of a partially recrystallized carbonate rock, as diffusion and the majority of metabolic processes are confined to only a small interval of the sediment column (upper ∼ 1 m of the sediment pile).Authigenic carbonate addi-tion is regarded to take place in the upper 0.1 m of the sediment column through precipitation from porewater DIC that is most severely perturbed by metabolic activity and largely unbuffered by exchange with carbonate sediments through dissolution and recrystallization.The evolved bulk-rock carbon isotope value is then calculated by using the mass balance of the primary and authigenic carbonate (f auth ) components.
An important parameter inducing changes in the vertical structure of a diagenetic redox profile, and ultimately carbonate carbon isotope composition during carbonate rock formation and stabilization, is the F OC .By systematically changing this parameter we test the sensitivity of the diagenetic baseline model to the OC sinking flux and the consequential organic matter remineralization trajectories, i.e., the dominant microbial communities involved in local in situ C-cycling.The effect of the OC sedimentary regime on the carbonate isotope system is probed by looking at the C isotope offset, the primary carbonate precipitate (settling through the water column) and the bulk-rock end-member ( 13 C primary−bulk ).In addition to the OC flux, we performed sensitivity experiments for the sedimentation rate (w and v), the concentrations of bottom water DIC and sulfate and the fraction of diagenetic and/or authigenic carbonate (f dia and/or f auth ) incorporated in the bulk rock.For each separate sensitivity experiment we return the changing parameter to the initial value, as defined for the baseline model (Sect.3.2.2),except for the parameter of interest and the OC flux.In addition, we tested the influence of bio-diffusion (D b0 ) and bioirrigation (α 0 ) under varying depth profiles of sediment mixing on 13 C primary−bulk .

First-order temporal trends
By compiling δ 13 C carb data and unifying them in a biochronological framework, we can distinguish first-order temporal features of the combined Iranian and the Meishan P-Tr (sub-)localities in China (Fig. 1).The subsampling routine prevents modifications through unequal sampling intensity; hence this temporal pattern is unbiased by the sampling strategy applied in individual studies.A gradual decline towards 4 ‰ lower δ 13 C carb can be discerned starting from the middle Changhsingian, with minimum values (0 ‰) reached in the earliest Triassic and a return to 2 ‰ higher δ 13 C carb commencing above the I. isarcica zone (unit I).Nonetheless, disparity in regional trends is visible as a double-peaked negative excursion marking the P-Tr transitional beds of Meishan, whereas such a signature is absent in time-equivalent boundary beds deposited at the Iranian sites.
By comparing the combined subsampled median trend lines of individual studies (constructed following the same statistical routine as outlined in Sect.3.1.3),it is possible to determine whether the first-order isotope trend is a consistent feature of the analysed stratigraphic sequences.A large portion of the compared individual datasets is marked with a high coefficient of determination (r 2 ; Fig. 3).We can conclude that the geochemical signal obtained is reproducible and, as such, observable in rock collected during separate sampling campaigns over several decades.However, the comparative study also suggests that the correlation coefficients are consistently weaker in the Meishan data relative to the Abadeh data.This lack of reproducibility might stem from limited stratigraphic coverage; see Gruszczyński et al. (2003), who report on a composite dataset of two sublocalities, Meishan D and Z, only encompassing Changhsingian strata.Sampling campaigns of reduced stratigraphic coverage might not be able to capture the first-order trend.Nevertheless, region-specific differences in signal reproducibility mark an overall greater disparity between individual Meishan δ 13 C carb records and might point to a higherorder of δ 13 C variability.

Residual carbon isotope variability
Subsampled and time-sliced median interpolations are depicted in Fig. 2. Trends with the highest confidence are highlighted by a well-defined white line, and the width of the CI is represented by a blue colour that becomes less intense with increasing data spread.Combined, these graphing features result in a more blurred image with a larger spread between the median values of individual subsampling routines (Sect.3.1.3).From these figures we can discern three features in the second-order δ 13 C excursions: (1) the residual δ 13 C variability is seemingly random or stochastic and defined δ 13 C excursions are unreproducible across lithological successions from different geographic locations (Iran) or studies targeting the same site (e.g., Meishan section; Fig. 1).
(2) There is a returning temporal pattern towards decreased fidelity of the median trend line (blurred white trend line) connected with increased CI dispersion (i.e., residual δ 13 C variability seen as less saturated blue colour tones) across the extinction horizon and ranging into the Early Triassic at both geographic locations.Increased maximum δ 13 C value ranges from less than 1 and 2 ‰ (IQR) or 2 and 3 ‰ (IPR) in the pre-extinction beds to more than 2 and 3 ‰ (IQR) or 5 and 8 ‰ (IPR) around the extinction horizon for Iran and China, respectively, further enforces the observation of a globally significant peak in residual δ 13 C variability.This pattern seems to recover after the P-Tr boundary, with generally higher confidence median trend lines and a trend to smaller maximum δ 13 C value ranges of less than 2 and 3 ‰ (IQR) or 5 and 8 ‰ (IPR) for Iran and China, respectively.(3) In addition, a significant regional offset can be discerned in δ 13 C value ranges; the Chinese profile displays systematically higher residual δ 13 C variability.
Studies focussing on boundary events (e.g., the P-Tr transition) tend to channel sampling effort at a focal point around the presumed faunal turnover and horizon of palaeoenvironmental change (Wang et al., 2014).Sampling effort is known to have a profound impact on studies that use fossil data to track animal diversity through time (Alroy et al., 2008;Mayhew et al., 2012).By correcting for sampling effects through subsampling (Sect.3.1.3),we cancelled out the potential effect of over-representative data accumulation on the temporal trend of stochastic residual δ 13 C carb variability.Sampling effort might also have affected fossil collecting, thereby biasing conodont biozonation schemes, which heavily relies on the first-appearance datum of certain fossil species (Wang et al., 2014;Brosse et al., 2016).As such, it would be conceivable that longer-duration biozones could capture a more temporally variable marine DIC δ 13 C compared to shorter-duration biozones.When evaluating the average thickness and duration of our chosen biochronological units, there is an apparent decrease in unit thickness and a shorter duration when approaching the extinction interval (Fig. 2).A comparison of biozone stratigraphic thickness and duration points to a relationship between higher δ 13 C carb dispersion and smaller and shorter-duration units (Fig. 2).This inverse relationship suggests that δ 13 C carb variability is not controlled by the increased potential sample size and the inherent risk of sampling a more temporally variable isotope signature.Hence, we exclude the applied sampling strategy and the biochronological framework as causal factors behind the observed temporal trend of stochastic δ 13 C carb variability.
By systematically changing F OC we see a defined relationship between bulk-rock C isotope alteration and the predominant mode of organic matter remineralization, e.g., microbial sulfate reducers, AOM or methanogenesis (Fig. 5).The sensitivity experiments further reveal that the fraction of diagenetic precipitate incorporation and bottom DIC have only a limited effect on the difference between the primary carbonate precipitate and the bulk-rock end-member C isotope value ( 13 C primary−bulk ).In contrast, marine sulfate levels modulate the total range of observed δ 13 C (∼ 2-3.5 ‰) and introduce a switch in the system, which narrows the range of OC accumulation over which the maximum deviation in 13 C primary−bulk occurs (Fig. 5).In a similar fashion, elevated sedimentation (v and w) causes a shift in peak 13 C primary−bulk values.However, this shift in peak values causes a modest increase in the total attained C isotope alteration with varying OC accumulation under higher sedimentation rates (v and w).
A systematic study of OC accumulation results in comparable trajectories of diagenetic C isotope modification on bulk rock through authigenic carbonate addition (Fig. 6).However, authigenic sea-floor cementation results in a 10 ‰ range of attainable bulk-rock end-member C isotope values, which is a wider range than obtained for diagenetic carbonate alteration.This style of carbonate cementation is also strongly controlled by seawater sulfate concentration.Lower than modern marine dissolved sulfate values (28 µmol cm −3 ) allow for a large range of δ 13 C end-member values to be generated under a smaller range of F OC (Fig. 6).As opposed to diagenetic carbonate alteration, heightened sedimentation (v and w) only causes a shift in the switch from positive to negative 13 C primary−bulk in the end-member rock but does not change the total range of C isotope alteration.However, both elevated bottom water DIC and a decreased size of the authigenic fraction diminishes the maximum range of 13 C primary−bulk .
Both bio-mixing and bio-irrigation by pre-extinction Permian benthic fauna (D b0 = 5 cm 2 yr −1 , z b = 2 cm, α 0 = 50 yr −1 , x irr = 1 cm) cause minimal modulation of 13 C primary−bulk .The influence of modern benthic fauna (D b0 > 5 cm 2 yr −1 , z b = 3 cm, α 0 > 50 yr −1 , x irr = 2 cm) would suppress δ 13 C modifications caused by OC-steered diagenetic carbonate alteration and authigenic carbonate addition.On the other hand, the absence of benthic fauna (postextinction situation) allows for an unconstrained impact of the previously cited OC-controlled trajectories of C isotope modification during bulk-rock formation and stabilization under P-Tr environmental conditions.

Simulation of a virtual carbon isotope time series
In order to understand the obtained temporal patterns in residual δ 13 C carb variability, a series of reactive transport models have been solved to steady state under varying OC sedimentation regimes.The virtual time series approach is an amalgamation of multiple sets of individual reactive transport model runs (50 model iterations = 1 set = 1 time unit) sliding across a timeline of length 100 and with time increments of length one (Fig. 7).
By performing sensitivity tests (Sect.3.2.3),we have established that there is a systematic relation between OC arriving at the seabed and the magnitude of bulk-rock C isotope alteration, which is expressed as 13 C primary−bulk (Fig. 5).Hence, the initial F OC can be estimated because the 13 C primary−bulk can be approximated by the C isotope offset between calcite from well-preserved Permian brachiopods and time-equivalent bulk-rock samples (∼ 1.0 ‰) (Schobben et al., 2016) and equates to 802 µmol cm −3 yr −1 .The F OC is subsequently modulated for each set (i.e., time unit) based on the observation that OC accumulation increases by a factor of 4 across the P-Tr transition (Algeo et al., 2013) and by linearly scaling this parameter to the observed residual δ 13 C carb variability (IQR) as obtained from Fig. 2, yielding a continuous range for F OC (802 ≤ F OC ≤ 3206 µmol cm −3 yr −1 ). .Sensitivity experiment for diagenetic carbonate alteration through equilibrium recrystallization with depth designed to investigate the forcing effect of (a) the fraction of diagenetic carbonate incorporated (f dia ), (b) sulfate content of overlying water, (c) DIC content of overlying water and (d) sedimentation rate (v and w) on the carbon isotope offset between the diagenetic end-member rock and the primary calcite over a range of F OC .Panels (e) and (f) depict sensitivity experiments for bio-irrigation (α 0 ) and bio-diffusion (D b0 ), respectively, on diagenetic altered carbonate for changing modes of sediment reworking and irrigation by biota under a normal OC flux (730 µmol cm −3 yr −1 ).An increasingly suppressed alteration of the end-member carbon isotope signal is observed for the pre-extinction (Dale et al., 2016) and modern (van de Velde and Meysman, 2016) depth profiles compared with sediment that is not inhabited by metazoans (post-extinction).
With regard to examining the effects of spatially distinct OC sinking fluxes and/or benthic fauna on the lateral distribution of OC, we forced individual model runs (within one time unit) with slightly deviating F OC values.The maximum attained F OC variability is constrained by adopting a lognormal skewed density distribution (R package emdbook; Bolker, 2008Bolker, , 2016) ) around the former established F OC , which represents the population's median value (Fig. 7), and by randomly selecting a value from the created distribution.A log-normal skewed density distribution has been chosen to represent natural variation in the F OC -population, which would be signified by variation that is skewed towards heavy sea-floor OC loading in rare instances, whereas most variation is smaller in magnitude and can modulate the flux towards both smaller and higher contributions.Temporal variations in this parameter are attained by modulating the width of the sample population ranging up to a factor of 1.5.This lateral OC dispersion parameter is also linearly scaled to the observed residual δ 13 C carb variability (IQR) as obtained from Fig. 2.This results in sample population variations .Sensitivity experiment for authigenic carbonate addition through sea-floor cementation designed to investigate the forcing effect of (a) the fraction of authigenic carbonate (f auth ) incorporated, (b) sulfate content of overlying water, (c) DIC content of overlying water and (d) sedimentation rate (v and w) on the carbon isotope offset between the diagenetic end-member rock and the primary calcite over a range of F OC .Panels (e) and (f) depict sensitivity experiments for bio-irrigation (α 0 ) and bio-diffusion (D b0 ), respectively, on authigenic carbonate addition for changing modes of sediment reworking and irrigation by biota under a normal OC flux (730 µmol cm −3 yr −1 ).An increasingly suppressed alteration of the end-member carbon isotope signal is observed for the pre-extinction (Dale et al., 2016) and modern (van de Velde and Meysman, 2016) depth profiles compared with sediment that is not inhabited by metazoans (post-extinction).
of 524 ≤ median absolute deviation ≤ 3889 µmol cm −3 yr −1 when combined with the initial F OC range for homogeneously and more heterogeneously dispersed sedimentary OC, respectively.The initial carbon isotope composition of F carb is based on the median δ 13 C carb values of the first time unit of the pooled Iranian dataset (Fig. 1).This value is subsequently corrected by −1.0 ‰ based on the systematic isotope offset between bulk carbonate and pristine carbonate from brachiopod shells (Schobben et al., 2016) in order to obtain the C isotope value of F carb .The systematic relationship of OC sedimentation with 13 C primary−bulk , as obtained from Fig. 5, yielded the primary δ 13 C of calcium carbonate particles arriving at the seabed for successive time units.
The processes governing increased authigenic sea-floor cementation in post-extinction strata are largely speculative.However, undisturbed substrates are a feature required to enable the growth of calcifying microbial mats, regardless of the mechanism inducing the global-scale proliferation of these sea-floor precipitates (e.g., Garcia-Pichel et al., 2004).Considering this universal control, it is justified to assume a direct link between bioturbation intensity and authigenic seafloor cementation.This assumption links both forms of carbonate rock formation and stabilization: equilibrium recrystallization and authigenic sea-floor cementation.The probability (p) of producing authigenic sea-floor cements has been assumed to vary linearly in the range 0.01 ≤ p ≤ 0.1 following the upper two terciles of the range attained by the previously defined OC dispersion parameter.This definition determines the response of introducing authigenic sea-floor cement δ 13 C to the sample pool.
In order to further simulate the mass extinction of marine metazoans across the end-Permian extinction on sedimentary conditions, the parameters that define bio-mixing and bioirrigation (Sect.3.2.2) are set to zero from time unit 60 onward (equivalent to biozone G; Fig. 1 and Table 1).In addition, it has been postulated that sedimentation rates increased over the P-Tr transition (Algeo and Twitchett, 2010), so the sedimentation rate (of v and w) was modulated at time unit 70 (equivalent to biozone H; Fig. 1 and Table 1).However, a conservative 2-times sedimentation increase is adopted, as the previously cited estimates largely hinge on calculations based on the duration of the respective stages, which are subject to continuous modification (Li et al., 2016).The f dia and/or f auth for both trajectories of carbonate rock stabilization is kept constant at 0.2.The final bulk-rock carbon isotope value is obtained by using the mass balance from the calculated primary and diagenetic and/or authigenic carbonate components.
A total of 5050 reactive transport models have been solved in order to match the length and resolution of the cumulative "real" δ 13 C time sequence from sites in Iran (Fig. 1).The resulting virtual time series of model runs forced by previously established value ranges for the OC sinking flux, sediment homogeneity and authigenic mineral addition is depicted in Fig. 8.Although these steady-state solutions will not always be truly representative of a dynamic sedimentary environment (e.g., sedimentation rate changes and sediment reworking), we conclude that the conceptualized model (Sect.2) and selected parameter values can explain trends in residual carbon isotope variability.

Stochastic carbon isotope signatures of the Permian-Triassic boundary beds
Global comparisons of P-Tr carbon isotope records have uncovered disparate δ 13 C trends and excursions (Fig. 1).These deviations are often encountered on refined stratigraphic and lateral scales, complicating interregional high-resolution correlations (Heydari et al., 2001;Hermann et al., 2010).Subaerial exposure and projected trajectories of bulk-carbonate stabilization under the influence of meteoric fluids and high water-rock ratios might be invoked to explain this disparity (Heydari et al., 2001;Li and Jones, 2017).Other studies have pointed to stratigraphic variations in predominant mineral rock composition as a source of δ 13 C carb modulations, citing the mineral-specific isotope offset between aragonite, dolomite and low-Mg calcite (Brand et al., 2012b;Heydari et al., 2013;Schobben et al., 2016;Li and Jones, 2017).
In this study, we provide an alternative explanation for disparate second-order δ 13 C carb excursions involving the observed stochastic residual δ 13 C variability on a refined spatial scale (Fig. 1).We regard this test to be of importance, as a transgression marks the boundary beds and Lower Triassic deposits (Fig. 1).Furthermore, subaerial exposure and consequential early rock stabilization with the high waterto-rock ratios (open system) of the studied localities is unlikely to have left an isotopic overprint.Under such circumstances, regionally extensive carbonate sediment dissolution through contact with an undersaturated solution could have sourced the underlying sediment pile with saturated fluids for cementing (Bathurst, 1975).However, evidence of such exposure in the form of karstification is absent in the studied locations (Leda et al., 2014;Yin et al., 2014).A recent study on the Meishan section suggests that recrystallization (i.e., zoned dolomite crystals) and the modification of the C isotope composition of certain stratigraphic levels was forced by a short-lived regression and consequential exposure to meteoric water (Li and Jones, 2017).However, the assignment of such petrological indices to specific diagenetic environments is often fraught with uncertainties, and it is not uncommon for alleged diagnostic features to occur in a range of diagenetic environments (Melim et al., 2002).Although a complete exclusion of isotope resetting by meteoric water is not possible, it is possible that a large portion of the global carbonate archive lithified without being subjected to meteoric fluids, therefore justifying the exploration of alternative modes of isotope resetting.
Systematically 12 C-depleted bulk rock relative to brachiopod calcite (Schobben et al., 2016) further counters the notion of a variable predominance of an aragonite precursor for Permian carbonate rock (Brand et al., 2012b;Heydari et al., 2013).Nonetheless, secondary calcite and dolomite additions might be a causal factor behind observed δ 13 C variability.However, these additions can also be the result of in situ microbially steered mineral formation, which incorporates metabolically produced DIC isotope signatures (Bontognali et al., 2014;Li and Jones, 2017).As such, our conceptualized model is not mutually exclusive with regards to carbonate polymorph compositional changes, but it does not require rapid global secular marine ion inventory changes.Moreover, our model outcome predicts δ 13 C variation driven solely by small biologically steered mineral contributions (f dia and/or f auth = 0.2; Fig. 6), whereas the recorded δ 13 C fluctuations are up to 10 times higher in amplitude than otherwise predicted for thermodynamically driven mineralspecific isotope compositional offsets (Brand et al., 2012b;Heydari et al., 2013).Note also that the mineral-specific isotope offset would require bulk-rock samples to display effectively complete (100 %) mineral compositional changes to a pure carbonate polymorph composition.In contrast, our model predicts a mechanism that can introduce temporal variation in spatially heterogeneous carbonate polymorph assemblies and associated isotope compositional differences through post-depositional processes driven by spatially and temporally variable OC accumulation comparable to modern ranges (Müller and Suess, 1979).Depleted marine dissolved sulfate concentrations would have exacerbated the sensitivity of carbonate rock towards only small deviations in this allochthonous C source (Figs. 5 and 6).Combined with a less well-mixed sediment and microbially steered cementation, this model can account for the entirety of the observed residual δ 13 C carb variability.

The carbon isotopic composition of the Permian-Triassic DIC reservoir
The first-order P-Tr negative carbon isotope trend can be explained by changes in the sources and sinks of the long-term (> 100 ky) carbon cycle, e.g., a reorganization of organic carbon burial and volcanic CO 2 outgassing (Renne et al., 1995;Kump and Arthur, 1999;Berner, 2002).Carbon-cycle partitioning between the deep and shallow ocean through the marine biological pump have been proposed to explain the depth-gradient isotope differences (Meyer et al., 2011;Song et al., 2013) or a rapid (< 100 ky) C isotope excursion (Rampino and Caldeira, 2005).On the one hand, the stochastic residual δ 13 C carb variability on a confined centimetremetre stratigraphic scale is inherently difficult to reconcile with transient perturbations in the dynamic C-cycle equilibrium enforced by repeated biological reordering of DIC between the shallow and deep ocean (i.e., consecutive events of global productivity increase and wholesale collapse).On the other hand, if lateral variations in marine DIC-C isotope composition account for δ 13 C carb variability, water column mixing would likely erase spatial compositional differences on smaller scales (< km scale).For example, residual δ 13 C carb variability recorded at Meishan (Fig. 1) on a centimetre to metre scale cannot be reconciled with the heterogeneity of the overlying water mass, and less well-mixed sed-iments become a more acceptable alternative.Although the geographically disparate C isotope signals of different sites in Iran (> km scale) could instead be accounted for by the heterogeneity of water masses, the temporal trend in the amplitude of the residual δ 13 C carb variability is similar for both China and Iran (Fig. 2).An overarching control of the hereinvoked OC-steered diagenetic mechanism is the most parsimonious model to explain this geographical coherence.Sea-floor methane clathrate dissociation, coal and organicrich sediment combustion through igneous intrusions and ocean Ni fertilization stimulating methane production by methanogens are viable sources of transient (< 100 ky) isotope-depleted carbon contributions to the atmosphereocean system (Erwin, 1993;Knoll et al., 1996;Krull and Retallack, 2000;Retallack and Jahren, 2008;Svensen et al., 2009;Rothman et al., 2014).Although all of the cited mechanisms are conceivable triggers for defined second-order δ 13 C carb excursions, they are less likely to create spatially divergent isotope trends and repeated temporally distinct fluctuations; instead, they essentially reflect bed-to-bed variation.The existence of transient carbon isotope excursions is also unlikely considering the predicted elevated DIC levels of an ocean without pelagic calcifiers, resulting in conservative behaviour of the DIC pool towards transient C isotope perturbations (i.e., increasing the system response time; Kump and Arthur, 1999;Zeebe and Westbroek, 2003;Rampino and Caldeira, 2005).Nonetheless, second-order (transient and depth-gradient) C isotope signals of a primary origin might still be imprinted in the isotope records.However, rock formed under the previously outlined marine depositional conditions (with intensified spatially heterogeneous trajectories of C isotope alteration) increases the chance of sampling spatially variable and diagenetically modulated δ 13 C signals, and these would obscure primary signals.

Perspectives on the Permian-Triassic carbonate archive: stratigraphy, biological evolution and the global biogeochemical carbon cycle
These findings put renewed constraints on the application of whole-rock δ 13 C as a high-resolution stratigraphic tool.Diagenesis forced by a recrystallization process under the influence of meteoric fluids, high fluid-to-rock ratios and oxidized terrestrially derived organics is generally considered to be the main driver of bulk-rock C isotope alteration (Brand and Veizer, 1981;Marshall, 1992;Brand et al., 2012a, b).Rocks without petrological evidence of meteoric fluids percolating through voids and interacting with sediment grains, e.g., palaeokarsts, blocky calcite and pendant cements, are usually regarded as pristine and minimally isotopically altered (Veizer et al., 1999;Flügel, 2010;Brand et al., 2012a, b).Regarding the observed P-Tr residual whole-rock δ 13 C variability and model predictions, we have shown that under certain marine and sedimentary situations marine carbonate rock diagenesis can introduce significant δ 13 C fluctuations to the whole-rock isotope composition.The insight gained here undermines the idea of using second-order δ 13 C fluctuations as stratigraphic markers, which are often confined to limited stratigraphic intervals (e.g., individual limestone beds) without a proper understanding of the potential effects of diagenetic carbon isotope modification.At least for the data presented here, we can demonstrate that shifts in bulk-rock δ 13 C smaller than classical biozonation schemes have to be considered with caution when they coincide with the sedimentological features of reduced bioturbation (e.g., lamination) and likely do not serve as a globally universal marker.Instead, long-term first-order trends cross-cutting formational boundaries (see Supplement) have a higher likelihood of representing secular variation in seawater DIC isotope composition and are largely unaffected by marine diagenetic processes (Fig. 1).Extrapolation of our model results to Mesozoic and Cenozoic bulk-carbonate records signified by comparatively reduced δ 13 C variability might reflect buffering of highly variable diagenetic C isotope contributions under elevated marine sulfate levels and physical sediment mixing by benthic fauna.
The carbonate chemical signal towards increased δ 13 C variability supports the global significance of three parameters during the deposition of P-Tr carbonate rock: reduced sediment mixing, authigenic precipitates and heterogeneous OC accumulation.Although we have postulated a connection between sediment lamination and C isotope variability, the relation of these observations to the end-Permian mass extinction is unconstrained.The cataclysmic mechanism that drove the extinction of benthic fauna is still a matter of debate, but bottom water anoxia is one of the main contenders.This could have been caused by decreased seawater oxygenation instigated by stagnating ocean circulation (Wignall and Twitchett, 1996;Winguth et al., 2015).Alternatively, a higher OC sinking flux might have led to increased O 2 drawdown through aerobic respiration and H 2 S production by sulfate-reducing microbes, leading to widespread marine euxinia (Meyer et al., 2008;Algeo et al., 2013;Schobben et al., 2015).However, δ 13 C variability complies most compellingly with at least continuing OC accumulation over the studied interval and thereby largely undermines scenarios of reduced primary productivity (e.g., Rampino and Caldeira, 2005).On the other hand, the proliferation of post-extinction authigenic sea-floor precipitates (e.g., thrombolites, stromatolites and fan-shaped structures) likely represents a symptom of reduced sediment disturbance and elevated carbonate saturation (Kershaw et al., 2009).Since some of the more conspicuous post-extinction sea-floor structures have been connected with microbial communities and distinct metabolism-mediated isotope signatures (e.g., Heindel et al., 2013), they might be equally important constituents of the observed spatial δ 13 C variability.Geographic differences observed as systematically higher residual δ 13 C variability at the Chinese localities are suggestive of a regionally distinct hydrographic and depositional setting.This regional feature could be linked with evidence of episodic Late Permian euxinic conditions in the South China basin (Grice et al., 2005;Cao et al., 2009).Comparatively higher and more spatially variable OC loading of the sea floor in the South China basin under these circumstances might be triggered by enhanced primary productivity or OC sinking fluxes (Algeo et al., 2011(Algeo et al., , 2013)).On the other hand, the depositional sites found in Iran lack evidence for bottom water anoxia and favour other drivers behind the marine faunal extinction, e.g., thermal stress or ocean acidification (Leda et al., 2014;Schobben et al., 2014;Clarkson et al., 2015).
In addition to capturing regional biological expressions of the global faunal disruption, our approach illuminates potentially important feedbacks of the Earth system that control the global-scale carbon cycle.Recently, attention has increased on the potential significance of authigenic carbonate production for the global carbonate reservoir (Schrag et al., 2013) and the addition of authigenic carbonate to Permian and Triassic carbonate-bearing sequences (Schobben et al., 2016;Zhao et al., 2016).The numerical exercise highlights the link between local sedimentary OC-cycling and bulkrock diagenetic stabilization, hinting at the potentially significant portion of remineralized OC sequestered by authigenic carbonates.Additionally, bulk-rock residual δ 13 C variability might serve as an indicator of infaunal biological activity, effectively preserving an ecological signal, which could in turn be an indicator of local sedimentary C-cycling.This sub-cycle is strongly controlled by bioturbation and controls the amount of buried organic matter and the supply of electron acceptors (Canfield, 1989;Canfield et al., 1993).These fluxes drive carbon removal from the exogenic carbon reservoir and potentially require a revision of generally accepted constraints on these budgets, including during events with postulated major C-cycle perturbations, such as the latest Permian extinction.

Conclusions
The δ 13 C carb record straddling the P-Tr boundary of Iran and South China can be decomposed in a first-order temporal trend with a negative trend marking the transition beds and superimposed second-order residual δ 13 C variability.The primary goal of the current study was to delineate the nature of the second-order δ 13 C fluctuations, which harbour a reproducible temporal pattern towards increased variability across the extinction beds.By investigating the origin of these second-order geochemical signals, it is possible to better constrain the limits of carbon isotope chemostratigraphy, most notably the time resolution on which it can be applied.Having established that second-order C isotope variability traces local effects rather than the global carbon reservoir, we discerned that mineralogical heterogeneity and meteorically steered cementation are unsatisfactory candidates to explain this spread in δ 13 C. Hence, the possibility of metabolism-steered carbonate addition is more likely, as the anaerobic microbial pathway introduces isotopically depleted (or enriched) carbon to the porewater and the same microbes are often connected with the production of diagenetic and authigenic carbonates.Moreover, organic carbon degradation occurs in a part of the sediment column signified by high reactivity of the primary carbonate sediments.By testing these assumptions with reactive transport models and with the construction of a virtual time series, we can recreate the observed variability in the end-member bulk-rock δ 13 C composition by inducing the heterogeneous spatial distribution of sedimentary OC, i.e., a reduction of physical sediment mixing or spatially divergent OC sinking fluxes.Authigenic sea-floor cementation, possibly mediated by microbial mat communities, is another source that might have introduced diagenetic overprints in the P-Tr boundary beds.On the other hand, these mechanisms can still preserve long-term trends in oceanic DIC carbon isotope composition.These findings strengthen the notion that bulk carbonate δ 13 C is a fruitful source of information that records secular trends in carbon cycle evolution and potentially informs us about the biosphere, notably infaunal animal activity, sea-floor microbial communities and OC delivery to the sea floor.

Figure 1 .
Figure 1.Compiled published and new data for multiple P-Tr rock sequences in Iran (a) and China (b).The individual carbonate carbon isotope values are placed on a dimensionless timeline that marries both geographic areas in the most acceptable biochronological scheme (Table1).The solid black line represents the subsampled median trend line.The dashed blue line depicts the seawater δ 13 C DIC curve as obtained from the time series simulation (Sect.5).The stratigraphic placement of the sea level changes and the extinction horizon as well as the biostratigraphic framework can be found in the Supplement.

Figure 2 .
Figure 2. Visually weighted data plot for the Iranian (a) and South China (b) P-Tr (sub)localities (Sect.3.1.3)depicting the median trend line, the IQR (50 %) and the IPR (95 %) δ 13 C value ranges.Trends with lowest fidelity are marked by a blurring of colours and less contrasted colours (based on the CI of multiple subsample routines).The dashed and the solid line represent the extinction horizon and P-Tr boundary,respectively.The saturation level of the green tiles in the upper two panels equals a more extended biozone thickness (0.12-32.00 m) and longer duration (0.02-1.00 My).Grey tiles represent intervals with no available biostratigraphic data.See Fig.1for colour legend of the stratigraphic units.

Figure 3 .
Figure 3.A comparative analysis of first-order temporal trends in δ 13 C carb from individual studies for the Abadeh and Meishan (sub-)localities.The saturation level of the individual tiles measures the correlative strength of stratigraphic trends obtained during individual sampling campaigns based on the coefficient of determination (r 2 ).

Figure 4 .
Figure 4. Diagenetic depth profiles of Late Permian sea-floor sediments of (a) oxidized and reduced solutes under a low shelf OC flux 500 µmol cm −3 yr −1 , (b) DIC and carbonate δ 13 C forced with a low shelf OC flux 500 µmol cm −3 yr −1 , (c) oxidized and reduced solutes under a high OC flux 1200 µmol cm −3 yr −1 , (d) DIC and carbonate δ 13 C forced with a high OC flux 7000 µmol cm −3 yr −1 .(e) Situational sketch of pre-extinction bulk-carbonate accumulation with a normal OC flux and active benthic fauna and (f) situational sketch of the postextinction sedimentation with elevated OC accumulation, removal of metazoan benthic fauna and consequentially reduced sediment mixing (artwork by Mark Schobben; http://cyarco.com).
Figure5.Sensitivity experiment for diagenetic carbonate alteration through equilibrium recrystallization with depth designed to investigate the forcing effect of (a) the fraction of diagenetic carbonate incorporated (f dia ), (b) sulfate content of overlying water, (c) DIC content of overlying water and (d) sedimentation rate (v and w) on the carbon isotope offset between the diagenetic end-member rock and the primary calcite over a range of F OC .Panels (e) and (f) depict sensitivity experiments for bio-irrigation (α 0 ) and bio-diffusion (D b0 ), respectively, on diagenetic altered carbonate for changing modes of sediment reworking and irrigation by biota under a normal OC flux (730 µmol cm −3 yr −1 ).An increasingly suppressed alteration of the end-member carbon isotope signal is observed for the pre-extinction(Dale et al., 2016) and modern (van de Velde and Meysman, 2016) depth profiles compared with sediment that is not inhabited by metazoans (post-extinction).
Figure6.Sensitivity experiment for authigenic carbonate addition through sea-floor cementation designed to investigate the forcing effect of (a) the fraction of authigenic carbonate (f auth ) incorporated, (b) sulfate content of overlying water, (c) DIC content of overlying water and (d) sedimentation rate (v and w) on the carbon isotope offset between the diagenetic end-member rock and the primary calcite over a range of F OC .Panels (e) and (f) depict sensitivity experiments for bio-irrigation (α 0 ) and bio-diffusion (D b0 ), respectively, on authigenic carbonate addition for changing modes of sediment reworking and irrigation by biota under a normal OC flux (730 µmol cm −3 yr −1 ).An increasingly suppressed alteration of the end-member carbon isotope signal is observed for the pre-extinction(Dale et al., 2016) and modern (van de Velde and Meysman, 2016) depth profiles compared with sediment that is not inhabited by metazoans (post-extinction).

FFigure 7 .
Figure 7. Schematic depiction of the work flow behind the virtual carbon isotope time series.Model sets consist of 50 separate reactive transport model solutions driven by a randomly selected F OC .The randomly generated organic OC flux values translate into a density distribution with a right-skewed tail.Increased F OC and a widening of the value range is used to simulate increased spatial OC heterogeneity.

Figure 8 .
Figure 8. Subsampled virtual carbon isotope time series with model results plotted as "visually weighted" trends.Each time unit incorporates 50 model iterations with a randomly chosen OC flux from a log-normal skewed density distribution.The time series simulates an increase of 4 times the average OC flux, a 1.5 times increase in the spatial heterogeneity of OC accumulation and increased likelihood of producing sea-floor authigenic carbonate to a maximum of 1 in 10.Bioturbation (expressed as bio-mixing and bio-irrigation) is modulated at the extinction horizon to mimic the effect of a dieoff among benthic organisms.In addition, sedimentation (v and w) is increased twofold over the P-Tr transition (see main text for sources).

Table 2 .
Biogeochemical reaction equations and kinetic rate expressions.

Table 3 .
Parameter values for the kinetic constants of the reactive transport model.