Geochemistry of deep Tunguska Basin sills, Siberian Traps: correlations and potential implications for the end-Permian environmental crisis

Abstract

A vast portion of the plumbing system of the Siberian Traps Large Igneous Province (STLIP) is emplaced in the Tunguska Basin, where borehole data reveal ubiquitous and abundant sills with great lateral extension. These intrusions intersect Cambrian–Ordovician evaporite, carbonate and siliciclastic series, and locally coal-bearing Permian host rocks, with a high potential for thermogenic gas generation. Here we present new geochemical data from 71 magmatic and 4 sedimentary rock samples from the Tunguska Basin center and periphery, recovered from 15 deep sills intercepted by boreholes. The studied samples are all low-Ti basalt and basaltic andesites, confirming absence of high-Ti and alkaline STLIP magmatism in the Tunguska Basin. The sills derive from picritic parental melts produced by extensive melting of a mantle source with recycled crustal components below a thinned lithosphere (50–60 km), within the spinel stability field. The mantle source was dominantly peridotitic, with enriched pyroxenitic domains formed by recycled lower crust, in agreement with previous models for the main tholeiitic STLIP phase. Limited amounts (up to 5%) of highly radiogenic granitoids or moderately radiogenic metapelites were assimilated in upper crustal magma reservoirs. After emplacement, sills intruded in Cambrian evaporites assimilated marlstones and interacted with the evaporitic host rocks, probably via fluids and brines. This is the first time that such process is described in subvolcanic rocks from all across the volcanic basin. The sills are correlated geochemically with the established chemostratigraphy for the on-craton STLIP lava piles and intrusions (Norilsk region). Sills correlated with the Morongovsky–Mokulaevsky Fm. and the Norilsk-type intrusions are the most voluminous, present all across the central Tunguska Basin, and bear the strongest evidence of interaction with evaporites. Massive discharge of thermogenic volatiles is suggested by explosive pipes and hydrothermal vent structures throughout the Tunguska Basin. We propose that this voluminous pulse of magmatism is a good candidate for the hitherto unidentified early intrusive phase of the STLIP, and may link the deep Tunguska basin sills to the end-Permian environmental crisis.

Introduction

The Siberian Traps Large Igneous Province (STLIP; ca. 251–252 Ma; Fig. 1), synchronous with the end-Permian mass extinction (Burgess et al. 2017; Dal Corso et al. 2020) and home to world-class Cu–Ni–PGE ore deposits (Ryabov et al. 2014), is one of the most intensely studied LIPs. However, due to limited accessibility, we still lack geochemistry and age information on sills constituting the non-exposed and deep-seated portion of the STLIP in the Tunguska Basin (TB).

Fig. 1

a Borehole locations (blue diamonds) are plotted on a map of the Siberian Traps LIP (modified after Svensen et al. 2018), along with locations of previously studied sills (yellow symbols). The Siberian Traps crop out in several regions. Four of these are roughly on-craton, namely Norilsk, Putorana, Maymecha-Kotuy and Nizhnaya Tunguska (translated to Lower Tunguska). Other expressions of the Siberian Traps magmatism appear in regions off-craton, such as the West Siberian Basin, and the Taimyr Peninsula in the North (Reichow et al. 2005, 2016). b Schematic logs of the investigated boreholes organized from W to E from 219 to 277; borehole 260 is placed last for its detached southern position. Sampled sills are indicated by black stars, sampled host rocks by white stars. c Lava formations and chemostratigraphic classification of STLIP tholeiites in the Norilsk area—simplified stratigraphic log (modified after Ivanov et al. 2008). Names of the formations: iv – Ivakinsky; sv – Syverminsky; gd – Gudchikhinsky; hk – Khakanchansky; tk – Tuklonsky; nd – Nadezhdinsky; mr – Morongovsky; mk – Mokulaevsky; hr – Kharaelakhsky; km – Kumginsky; sm – Samoedsky. Nrl – Norilsk 1 intrusion, positioned after Pavlov et al. (2019)

The intrusive components of LIPs are attracting increasing attention from the scientific community since thermogenic production of volatiles from the country rocks by magma/host rock interaction was recognized as volumetrically significant in the appraisal of LIPs gas emissions (Svensen et al. 2004, 2009, 2018). The numeorous pipes/vents and volcaniclastic products within the STLIP suggest a link between its sub-volcanic plumbing system and surface products (Jerram et al. 2016; Polozov et al. 2016), adding further evidence for a significant host basin influence in the STLIP evolution. Volcanic basins and LIP plumbing systems are, therefore, key players for linking magmatic events, gas production, and environmental crises. For the atmospheric impact of a LIP, the tempo of volcanism is a crucial factor, perhaps more so than a cumulative appraisal of the volumes of volcanics and volatiles involved (Sprain et al. 2019). As much as they can be considered geologically brief events, LIPs usually last much longer than the related extinction intervals (< 100 kyr). Hence, the extinction-triggering event is probably a shorter subset of an entire LIP episode (Burgess et al. 2017). Recent U–Pb geochronology was focused on providing a high-precision reassessment of the STLIP, looking at timing of intrusive vs. effusive activity (Burgess and Bowring 2015; Burgess et al. 2017; Augland et al. 2019). Most ages obtained for the intrusive phase of the STLIP post-date the extinction interval, with the exception of the Norilsk 1 and Karaelakh intrusions (Burgess et al. 2014; Burgess and Bowring 2015). However, the start of intrusive activity is not constrained, hence it is possibile that the extinction-triggering part of the Siberian Traps is a yet undetected or unsampled early intrusive batch (Burgess et al. 2017). Increased penetration of the subsurface with boreholes (Svensen et al. 2018) recently gave a more complete picture of the intrusive STLIP in the TB. Sills comprise a significant volume of the STLIP (between 37 and 44%; Vasiliev et al. 2000; Ross et al. 2005), and they intercept volatile-rich lithologies (evaporite, marl, limestone and coal-bearing siliciclastics) with high potential for thermogenic gas discharge into the end-Permian atmosphere (Svensen et al. 2018). The TB sills are excellent candidates for representing the undetected intrusive pulse invoked by Burgess et al. (2017), but evidence for their interaction with host rocks is scarce (Svensen et al. 2009), and no geochemical data is yet available for deep sills from localities covering a vast part of the Tunguska Basin. We have analyzed sill samples (chips and cores) from boreholes drilled to 4 km depth in 5 localities distributed across the central part of the TB. This geochemical study sheds light on the petrogenesis of these sills and their interaction with crustal rocks. Hinging on proposed geochemical correlations with the chemostratigraphy etablished in the Norilsk region, we discuss the timing of relative emplacement of the sills network and its implications on the STLIP thermogenic gas production scenarios.

Geological context for the Tunguska Basin and chemostratigraphy of the Siberian Traps

The Siberian Traps LIP includes large volumes of basaltic lavas (Fedorenko et al. 1996; Reichow et al. 2005) and volcaniclastic products (Jerram et al. 2016), as well as a vast sub-volcanic plumbing system of dikes and sills, connected to explosive pipe and hydrothermal vent structures (Fig. 1; Polozov et al. 2016; Svensen et al. 2009, 2018). The STLIP crops out in four main regions, mostly within the Siberian Craton (Sharma 1997): Norilsk, Putorana, Maymecha-Kotuy and Nizhnaya (Lower) Tunguska (Fig. 1). Off-craton, the STLIP extends to the West Siberian Basin (buried basalts; Saunders et al. 2005; Reichow et al. 2005) and in the Taimyr peninsula (Reichow et al. 2016; Augland et al. 2019). Volume estimates for the STLIP products vary widely, with a minimum of 2 million km³ (Vasiliev et al. 2000). Preserved and original volumes are difficult to constrain due to the remoteness of most outcrops and to the fact that a significant part of the province resides below the surface (Reichow et al. 2002, 2005, 2009; Svensen et al. 2018). The short-lived main (tholeiitic) pulse of volcanic and subvolcanic magmatism was dated to ca. 252.3–251.3 Ma (Burgess and Bowring 2015; Augland et al. 2019), and straddles the extinction interval at the Permo-Triassic boundary and the concomitant negative δ¹³C excursion (dated in Meishan to 251.941 ± 0.28–251.88 ± 0.28 Ma; Burgess et al. 2014).

The Tunguska volcanic basin is an extensive system of interconnected sill complexes and associated hydrothermal vents. The sills intruded the vast and deep (3–12.5 km) Tunguska sedimentary Basin, composed of Riphean to Early Permian sedimentary series, including evaporitic, siliciclastic, carbonatic, and hydrocarbon-bearing lithologies, in addition to thick coal deposits. The basement below the TB does not crop out, but Proterozoic gneiss and leucogranites were sampled as xenoliths (Czamanske et al. 2000) by explosive pipes close to Norilsk. Sills are abundant throughout the TB, and locally constitute up to 65% of the stratigraphy (Kontorovich et al. 1997; Ulmishek 2001; Svensen et al. 2018). These sills mainly crop out along the margins of the province and along the main rivers, but are intersected by numerous stratigraphic, exploratory and commercial boreholes. Svensen et al. (2018) reported sills in 93.5% of 284 boreholes (including the 6 boreholes studied here), concluding that sills are present everywhere beneath the TB, although they may pinch out laterally.

The effusive STLIP started as a low-volume subalkaline high-Ti phase (TiO₂ > 2 wt%; Lower Sequence), continued as a low-Ti tholeiitic series (Middle and Upper Sequence) and closed with a very heterogeneous alkaline phase. Tholeiitic and alkaline magmas from the Siberian Traps originated from different mantle lithologies and sources and followed separate paths and processes en route to the surface (Arndt et al. 1998; Elkins-Tanton et al. 2007; Sobolev et al. 2011). This interpretation is mainly based on studies of the lava sequences in the Norilsk area, where chemostratigraphic units and magma types were defined based on major and trace element and isotope data for the lavas (Lightfoot et al. 1990, 1993; Wooden et al. 1993; Fedorenko et al. 1996), and for the local intrusive bodies (Naldrett et al. 1992; Hawkesworth et al. 1995; Arndt et al. 2003). Similar signatures can be tracked over long distances in the STLIP (Reichow et al. 2005; Augland et al. 2019), and it is not uncommon for LIP magmas to show uniform geochemical compositions over large areas, even over > 1000 km distance (Marzoli et al. 2018; Bédard et al. 2020). We note that as more data on the STLIP became available, independent names for different regions were assigned to the various magma units (cf. Fig. 2 in Sharma 1997 and Sibik et al. 2015). Here, we use the nomenclature defined from the Norilsk region studies, which can be taken as type-locality because it represents the most complete section of volcanic rocks in the STLIP.

Fig. 2

Mineral-whole rock Fe-Mg equilibrium diagrams for olivine (a) and clinopyroxene (b). Mineral compositions (Fo and Mg#) are compared to whole-rock Mg#, calculated with FeO set to 85% of the total iron in the whole-rock. The curves display the range of equilibrium compositions between mineral and melt for Fe-Mg partition coefficients of 0.30 ± 0.03 for olivine (Roeder and Emslie 1970) and 0.25 ± 0.05 for clinopyroxene (cf. Bédard 2014). The mineral chemistry full dataset is available in Supplementary information 2

Near Norilsk, eleven lava formations (see summary table in Fig. 1c; Fedorenko et al. 1996; Lightfoot et al. 1990, 1993; Wooden et al. 1993; Hawkesworth et al. 1995) and five major magma types (Naldrett et al. 1992) were defined. Lightfoot et al. (1990, 1993) grouped the lavas into Lower (high-Ti), Middle-Transitional (low-Ti), and Upper (low-Ti) sequences (Fig. 1c). The Lower, high-Ti sequence comprises the alkaline/subalkaline, high magnesian Ivakinsky, Syverminsky and Gudchikhinsky formations. Above these, the Middle-Transitional low-Ti sequence includes the Tuklonsky and Nadezhdinsky lava suites. These are followed by the Upper sequence, including five suites of mainly tholeiitic basalts (Morongovsky to Samoedsky fromations) with uniform geochemical signatures, and representing the main tholeiitic pulse of the STLIP. The intrusions found in the Norilsk region have been classified into five types and correlated with the lava piles (Naldrett et al. 1992; Hawkesworth et al. 1995; Arndt et al. 2003; Pavlov et al. 2019; Yao and Mungall 2021). The most prominent intrusive types are the dolerites and differentiated intrusions (type 3), spanning the entire Norilsk region, the ore-bearing Norilsk-type intrusions (5A), and the Lower Talnakh-type (5B; Naldrett et al. 1992). Type 3 intrusions are further subdivided and correlated with specific Upper series formations, as is the case for the Ambarninsky type II dolerites, correlated with the Kharaelakhsky Fm. (Hawkesworth et al. 1995). Lower Talnakh-type intrusions are geochemically similar to the Middle sequence Nadezhdinsky basalts (enriched trace elements ratios and Nd isotopes, low Cu), to the West Siberian Basin and Taimyr basalts (Reichow et al. 2005, 2016). Norilsk-type intrusions are geochemically akin to the Upper Sequence Morongovsky–Mokulevsky basaltic suites, with which they are coeval (Pavlov et al. 2019; Burgess and Bowring 2015). The last phase of the STLIP, intrusive and effusive, is strongly alkaline, and so far documented only in the Maymecha-Kotuy region (Arndt et al. 1998; Sobolev et al. 2009b; Fig. 1) and in the Taimyr peninsula (Vernikovsky et al. 2003; Augland et al. 2019).

Despite some still unexplained features such as the lack of regional uplift, there is general consensus on a plume-derived origin for this LIP (Sobolev et al. 2011 and references therein). However, different petrogenetic histories are reflected by the geochemistry of the STLIP rocks. The Lower sequence high-Ti lavas were interpreted as derived from small degree deep melting of pure pyroxenite in the presence of garnet (Sobolev et al. 2011). Lavas of the Upper and Middle-Transitional low-Ti sequences reflect shallower and more extensive melting of a dominantly peridotitic mantle source. These magmas show arc-like features (e.g. low Nb, high Pb, Sr, Ba, fairly radiogenic Sr and Nd isotopes) that have been explained as either inherited directly from a subduction-modified lithospheric mantle (Fedorenko et al. 1996; Lightfoot et al. 1993; Ivanov et al. 2008), or as an effect of crustal contamination in deep crustal magma chambers (Wooden et al. 1993; Arndt et al. 2003; Reichow et al. 2005). High temperature deep melting of a volatile-rich source produced the last alkaline phase of the STLIP (Elkins-Tanton et al. 2007).

Materials and methods

Sampled borehole

The samples analyzed in this study were collected from five boreholes (Fig. 1b; coordinates in Supplementary information 1) housed in the core storage of the Siberian Scientific Research Institute of Geology, Geophysics and Mineral Resources (Novosibirsk). We logged, sampled and studied selected intervals in the Nizhneimbakskaya 219 borehole (Total depth TD 4000 m), Zheldonskaya 260 (TD 4509 m), Suringdinskaya 274 (TD 2679 m), Majgunnskaya 275 (TD 3500 m) and Verkhneilimpejskaya 277 (TD 3030 m). The sills intercepted by each borehole are identified by letters, starting from the shallowest sill as “A” (Fig. 1b). A total of 71 samples were studied from 15 individual sills. Borehole 219 was drilled at the western limit of the TB, where only intrusive products of the STLIP are preserved. Five sills were intercepted, for a cumulative thickness of 710 m. All but the shallowest one (not available for sampling) were emplaced in Cambrian evaporites, including halite, carbonates and anhydrite. Boreholes 274 and 275 were drilled SE of 219, where intrusive and volcanoclastic products are preserved. Borehole 274 intercepted 4 sills. The thickest sill (minimum 310 m; not sampled) crops out and intrudes Permian clastic, organic-rich sediments. The three sills below are thinner (total thickness 145 m), and intrude Cambrian evaporites. Seven sills are found along borehole 275 (total thickness of 478 m). The shallowest sill (A) intrudes Ordovician marls. The remaining 6 sills intrude marl- and dolostone-dominated Cambrian series. Only the deepest sill was available for sampling. Borehole 277 was drilled ca. 500 km to the east of 275, in an area with abundant pyroclastic deposits and breccia pipes, and it intercepted 5 sills (Svensen et al. 2018 reported 707 m as cumulative sill thickness, which is a minimum value because the borehole stops within sill 277E). The shallowest sill is emplaced in Silurian marls and is recorded in the logs as a single sill. However, sedimentary cuttings recovered between 325 and 440 m, suggest that sill 277A is in fact a composite intrusion of two sills separated by sedimentary rocks (hereinafter referred as sills 277A1 and 277A2). The other 4 sills are emplaced in Cambrian evaporites (halite and dolostone). Borehole 260 was drilled in the south, ca. 500 km south of 277 and 200 km north of Bratsk, in an area heavily pierced by pipe structures (Svensen et al. 2018). Two thin sills (total thickness 107 m) were intercepted within Cambrian evaporites (halite and dolostone).

Analytical techniques

Samples selected for geochemical analyses were either core fragments or cuttings selected from the freshest available material (Fig. S1a, b;e; note that S refers hereinafter to figures included in Supplementary information 1). Any visible alteration patches and fracture sides were eliminated by saw from the cores. Cuttings of 1–5 mm were carefully hand-picked under binocular microscope after sieving to avoid fragments of other origin. Chips coarser than 2 mm were reduced in size under a fly press, using plastic bags and plexiglass plates to avoid contact with the press metal. Samples were repeatedly washed with deionized water in ultrasonic bath and milled using an agate ball-mill.

Whole-rock major and trace elements were measured by Inductively Coupled Plasma Optical Emission (ICP-OES) or ICP Mass Spectrometry (ICP-MS) using the 4Lithores method at Activation Laboratories Ltd (Canada). Sample powders were mixed with a flux of lithium metaborate and tetraborate and fused in an induction furnace and mixed with 5% nitric acid until completely dissolved. Major elements and selected trace elements (Ba, Be, Sc, Sr, V, and Zr) were analyzed by ICP-OES (Thermo Jarell-Ash ENVIRO II ICP or a Varian Vista 735 ICP), whereas the remaining trace elements were analyzed by ICP-MS (Perkin Elmer Sciex ELAN 6000, 6100, or 9000 ICP/MS). Along with the samples, certified reference materials, a method reagent blank, and sample replicates were also analyzed (Supplementary information 2). Calibration was performed using six prepared USGS and CANMET certified reference materials for the ICP-OES analysis, and 12 for the ICP-MS analysis. For the ICP-OES analysis, one of the standards was analyzed for every group of ten samples, while the ICP-MS instrument was recalibrated every 40 samples. Analytical uncertainties for major elements are estimated at < 5%. For trace elements, uncertainties are 100% at the detection limit, 20% at 10 times the detection limit, and 10% at 100 times the detection limit (Supplementary information 2). Loss on Ignition (LOI) was determined by igniting sample powders in a chamber furnace at 1000 °C for 2 h, and weighing the samples before and after heating.

Whole-rock Sr–Nd isotopic compositions were determined for a subset of 35 samples on the same sample aliquots analyzed for major and trace elements, using a Thermal-Ionization Mass Spectrometer (TIMS; Finnigan MAT 262) at the Bergen Geoanalytical Facility, University of Bergen, following the procedures described in Pin et al. (1994). Pulverized samples were heated in a chamber furnace to 1000 °C for 1 h and decomposed in HF–HNO₃ in Savillex beakers on a heating plate. The dissolved and dried samples were re-dissolved in 2 N HNO₃ prior to chromatographic separation. Sr and Nd were separated from matrix elements by ion selective extraction chromatography on three different chromatographic materials supplied by Eichrom. Rare-Earth Elements (REE) were first retained on TRU Spec columns. In the next step, the retained REE elements were purified on LN Spec resin to selectively recover Nd. After removal and recovery of REE, Sr was extracted from matrix effluent by the ion specific resin, Sr Spec. Nd and Sr isotopic rations were corrected for mass fractionation using ¹⁴⁶Nd/¹⁴⁴Nd ratio of 0.7219 and ⁸⁸Sr/⁸⁶Sr ratio of 8.375209. Repeated measurements of the Strontium Carbonate Standard (SRM987) at the time of analysis (May 2018 and September 2018) yielded an average ⁸⁷Sr/⁸⁶Sr ratio of 0.710229 ± 8 (2SD; n = 3) and 0.710237 ± 13 (2SD; n = 7), respectively. Hence, ⁸⁷Sr/⁸⁶Sr measured values for the samples were adjusted for daily variations between the measured values of SRM987 and the reference ⁸⁷Sr/⁸⁶Sr value of this standard (0.710240). Repeated isotopic measurements of the La Jolla Reference Standard at the time of analysis yielded an average ¹⁴³Nd/¹⁴⁴Nd ratio of 0.511850 ± 11 (2SD; n = 2) and 0.511850 ± 9 (2SD; n = 5), equivalent to the accepted reference value of 0.511850 for La Jolla, negating the need for further correction. Analytical errors (2σ) for the measured Sr isotopic ratios are reported in Table 1 and vary between 0.000006 and 0.000009. The uncertainties for the values re-calculated to 251 Ma taking into account radiogenic in-growth are propagated including uncertainties on trace element determinations for Rb and Sr and on the decay constant. They vary between 0.000011 and 0.000089. Analytical errors for ¹⁴³Nd/¹⁴⁴Nd are also reported in Table 1 and vary between 0.000005 and 0.000013. For two sedimentary samples ¹⁴³Nd/¹⁴⁴Nd values were not measured due to extremely low Nd concentration.

Table 1 Major and trace elements and Sr–Nd isotopic compositions for selected TB dolerites

Detailed chemical traverses and spot analyses on selected crystals were obtained using a Cameca SX100 electron microprobe at the University of Oslo, operated in wavelength-dispersive (WDS) mode at 15 kV, 15 nA for most mineral phases and at 10 nA for plagioclase, with a focused electron beam and a counting time of 10 s on peak. Matrix corrections were performed according to the PAP procedure implemented in the Cameca software (Merlet 1994). Calibration standards were as follows: wollastonite (Si, Ca), orthoclase (K), albite (Na), synthetic Al₂O₃ (Al), pure iron metal (Fe), pyrophanite (Mn, Ti) synthetic periclase (Mg), eskolaite (Cr), nickel oxide (Ni), synthetic alforsite (Cl), fluorite (F), and Durango apatite (Ca, P, F). Where possible, F, Cl, Na, and K were analyzed first.

Petrography and mineral chemistry

Sill samples vary in texture from microporphyritic (Fig. S1f, rare), with microphyric groundmass, to fine- or medium-grained ophitic dolerites (Fig. S1c, d; common), with a tendency for clinopyroxene to form oikocrysts. Coarser grain size is observed for the deepest sills (e.g. 219E; 277E; Fig.S1a), where some medium-grained gabbros are found (Fig. S1a). For the sills with numerous available samples (e.g. 219C, 277B, 277C), a shift in grain size within a single sill is recognized, from fine-grained or porphyritic at the margins to medium-grained, gabbroic, at the center, although never exceeding 5 mm grain size for the largest crystals.

The sill mineralogy is typically doleritic, with ubiquitous clinopyroxene and plagioclase, frequent olivine and rare orthopyroxene. Oxides (ilmenite and Ti-magnetite) are present, in very large (occasionally poikilitic) or small crystals, sometimes showing dendritic or resorbed textures. Magmatic sulfides are rare and small, but pyrite and chalcopyrite are regularly observed occupying large secondary veins, or pockets between larger crystals in coarse-grained rocks (arrow in Fig. S1a), or at the contact with the host sedimentary rocks, where they can reach mm size (arrow in Fig. S1b). The large clinopyroxene oikocrysts are sometimes visible macroscopically as rounded domains (Figs. S1d, S2). Oikocrystic augite encloses randomly oriented plagioclase laths and very rarely olivine crystals, but never opaque minerals (Fig. S1c, d). In some samples, plagioclase of different generations appears, as small crystals included in the pyroxenes, but also as larger, zoned crystals. In various samples, evolved micropegmatites occupy late-stage interstitial pockets isolated among larger crystals. The mineralogy of these domains is typical of evolved (granitic) melts, with quartz, (Cl-rich) biotite, chlorite, K-feldspar, apatite and small occurrences of incompatible-element-bearing minerals (e.g. zircon, baddeleyite, sphalerite). The crystallization sequence of the primary dolerite assemblage as suggested by the textures is: olivine, plagioclase, pyroxenes, and opaque phases, followed by crystallization of the evolved interstitial domains, for which euhedral quartz, high-Cl in biotite and the presence of chlorite suggests that crystallization straddled orthomagmatic and hydrothermal stages. Detailed mineral chemistry of the contacts with host rocks and of the late-stage pockets will be addressed in a separate contribution.

In several dolerites, minor alteration has caused substitution of plagioclase by clay minerals or sericite and local growth of iddingsite and serpentine at the expense of olivine. Chlorite and amphibole sometimes replace pyroxene. Chlorite appears also as a substitute for biotite. Late-stage hydrothermal fluids precipitated calcite, scapolite, amphiboles and phyllosilicates in millimetric to centimetric veins. Two samples from a magmatic enclave, possibly an autolith, were found within sill 219D (Fig. S2). This enclave is relatively coarse-grained (up to 5 mm crystals), with a melanocratic portion (sample 219–2007.0A), dominated by dark blades of amphibole and few interstitial pink and white feldspars, and a leucocratic portion (sample 219–2007.0B), very altered, volumetrically dominated by feldspar, and showing some porosity. Thin sections were made through the dolerite–enclave contacts, for both portions (Fig. S3), highlighting a magmatic (possibly cumulitic) texture for the enclave, which is also strongly altered with sericitized plagioclase, and pyroxene substituted by chlorite and amphibole. The surrounding dolerite is fine-grained and well preserved, except in areas close to fractures, where chlorite and secondary biotite have developed. Four sedimentary host rocks from borehole 260 were analyzed. Sample 260/2200 is an impure dolomitic limestone, 260/2998.45 and 260/3000.3 are dolostones, and 260/3001.2 is a mixed anhydrite/halite evaporite. These samples do not show strong compositional or textural evidence of contact metamorphism, but many contact zones in hand sample (from borehole 260 as well as others studied here) show features of contact phenomena including re-crystallized reaction zones, fracturing at contact margins, and concentration of sulfides (Fig. S1b).

Olivine crystals were analyzed in 20 samples (Supplementary information 2). Olivine shows large compositional variations, ranging in forsterite between Fo₈₆ and Fo₁₉. Most olivine crystals have compositions around Fo_50–70 in our dolerites, and most crystals are unzoned or show mild zonations. Most olivine crystals are in disequilibrium with the whole rock, for a Fe/Mg olivine/melt partition coefficient of 0.30 ± 0.03 (Roeder and Emslie 1970; Fig. 2a), except those from sills 260A-B and 275G. Normal zoning is observed for sample 260/2207.7 (Fo_75–29), 219/1973.9 (Fo_34–27) and 219/1976.2 (Fo_62–49). Olivine crystals from sample 260/3025-3030 are inversely zoned (Fo_54–74), with the rims in equilibrium with the whole rocks and the cores less Mg-rich. Olivine crystals from rapidly cooled (chilled margin) samples from sill 260A-B contain melt inclusions.

All analyzed dolerites contain pyroxene. Augite (Wo_17–50En_49–46Fs_42–9) is the most common pyroxene phase, with rare occurrence of subcalcic augite. Considering a Fe–Mg partition coefficient of 0.25(± 0.05; Bédard 2014), most augites are not in equilibrium with the whole rock, having lower Mg# than predicted by equilibrium conditions, except those from 260/2210.3 and 260/3008.2 (Fig. 2b). Most clinopyroxene crystals are unzoned, some present a slight normal zoning (Fe enriched at the rims, e.g. 260/2207.7 from Wo₄₄En₃₉Fs₁₆ to Wo₃₇En₂₆Fs₃₆), and augites from 260/3008.2 show a complex zoning (decreasing Ca–Mg and increasing Fe from core to intermediate portions of the crystal, followed by a strong inverse zoning). Similar inverse zoning was observed for olivine from the same rock. Augites from 260/3008.2 are rich in Na₂O (0.37–1.46 wt.%, as opposed to < 0.45 for the other crystals, some groundmass crystals excepted), TiO₂ (up to 2 wt.%), and CaO (20.1–25.4 wt.%). Orthopyroxene (Wo_1-5Fs_42-71En_56-27) is present in samples 219/3891.7, 219/3930.6 and 277/495 and was observed as microcrystalline matrix phase in sample 260/3025–3030. We calculated crystallization pressures of 8 ± 3kbar for clinopyroxene from sill 260B and 2 ± 1 kbar for clinopyroxene from sill 260A (Neave and Putirka 2017).

Plagioclase is ubiquitous and was analyzed in 22 samples, with compositions ranging from An₈₂ (bytownite) to nearly pure albite (An₂). Strong normal zoning is common. Compositions approaching pure albite are exclusively found at the rims, especially in crystals bordering late-stage micropegmatitic pockets. In contrast to what is observed for olivine and clinopyroxene, crystals from sample 260/3008.2 are normally zoned. K-feldspar was analyzed within the evolved interstitial pockets.

Whole-rock geochemistry

Major elements

Whole-rock major elements were anlayzed on 69 dolerites, 2 gabbroic enclaves and 4 sedimentary rocks (representative analyses in Table 1; complete dataset in Supplementary information 2). With SiO₂ ranging between 47.56 and 53.09 wt% and total alkali (Na₂O + K₂O) between 2.37 and 4.56 wt%, most dolerites classify as basalts on a total alkali-silica (TAS) classification diagram (Fig. 3c; Le Maitre 2002), with four samples classifying as basaltic andesites. All but 5 dolerites are subalkaline, plotting at or below the alkaline/subalkaline separation curve (Irvine and Baragar 1971). Rocks from sills 219B, 219C, 219E, 277A and 277D depict trends of increasing alkali with increasing silica, whereas sills 260A, 260B and 275G show a negative correlation. Loss on Ignition (L.O.I.) for most rocks is low (< 2 wt%), with the exception of samples 260B-3004.9 (2.57 wt%) and 275G/2634.2 (3.48 wt%). L.O.I. increases with increasing Na₂O and K₂O concentration, but is not correlated with CaO content. The dolerites are mildly to strongly evolved basalts, with MgO from 3.68 to 9.60 wt%, (most samples have MgO > 5wt%), resulting in Mg# from 64 to 34 (Mg# = molar 100 × [(MgO/40.3)/(MgO/40.3 + (FeO_tot × 0.85/71.85)], assuming Fe³⁺/Fe_tot = 0.15). All the dolerites show low TiO₂ (< 2 wt%; Fig. 3a), except two evolved samples from the upper portion of sill 277D and the melanocratic enclave. Dolerites from 260B have higher TiO₂ than the other rocks for a given MgO. MgO correlates negatively with TiO₂, Fe₂O_3tot, P₂O₅ and SiO₂, and positively with Al₂O₃. CaO concentrations are too scattered to depict any trend (Fig. 3b; Fig.S4). The major element compositions of our dolerites fall within the field of previously published data for Siberian Traps sills (Fig. 3), and with the lavas of the Middle and Upper sequence (Fedorenko et al. 1996).

Fig. 3

a) Whole-rock variation diagrams of TiO₂ and CaO against MgO. Plotted are also the 2 specimens from the coarse-grained enclave from sill 219D. Previously studied sills compositions are also shown, either as a compositional field (Hawkesworth et al. 1995; Egorova and Latypov 2013; Kiselev et al. 2014; Sibik et al. 2015; Krivolutskaya et al. 2020) or with data points (Usol’skii and Tolstomysovskii sills from Ivanov et al. 2008; sills of the Bratsk area from Black et al. 2012). Data by Egorova and Latypov (2013) extend to higher MgO than our dolerites (up to 12 wt% MgO), due to the cumulitic character of many of their samples. Symbols grow darker with increasing depth of the sills. A full set of major element variation diagrams is available in Supplementary information 1 (Fig. S4). Liquid lines of descent are calculated for 60% fractional crystallization after emplacement of the sills, starting from the most primitive dolerites with the Magma Chamber Simulator (Bohrson et al. 2020; Heinonen et al. 2020) at 1kbar at anhydrous and hydrous conditions. Fractionating phases for anhydrous conditions from 219/1410–15: 5% ol, 30% cpx, 25% plg; from 260/2225–30, 30% plg, 30% cpx. Fractionating phases for hydrous onditions from 219/1410–15: 13% ol, 22% cpx, 23% plg, 2% spinel; from 260/2225–30: 8% ol, 24% cpx, 24% plg, 4% spinel. c) The analyzed dolerite samples (69) are plotted on a portion of the TAS classification diagram (Le Maitre 2002). Analyses are recalculated on a volatile-free basis. The dashed discriminating line for alkaline character is from Irvine and Baragar (1971). A mixing line between 260A dolerite and 10% impure salt (Pang et al. 2013) is shown. Tick marks every 1% mixing

Trace elements

The TB dolerites show variable REE abundances (La_CH = 18–73, Yb_CH = 12–31 times chondritic concentrations; McDonough and Sun 1995), mostly with smooth, gently sloped patterns broadly parallel with each other (Fig. 4; Fig. S5a–c). Most samples show (La/Yb)_CH from 1.46 to 2.54, (La/Sm)_CH from 0.90 to 1.47, and variable Eu anomalies (0.6–1.2; calculated as Eu* = Eu_CH/(Sm_CH × Gd_CH)^0.5; Supplementary information 2). Gabbroic, evolved rocks sampled at the top of sills 219B, 219E, 277D and 277E are the most enriched in REE (Fig. 4a–b), with mildly negative Eu anomalies (0.80–0.90). In sill 219D, the central samples and the gabbroic enclaves have the highest REE concentrations. Sills 260A and 260B show very uniform within-sill compositions, but 260B rocks are more enriched. Sill 277A1 marks a clear exception (Fig. 4b), with much higher LREE concentrations and higher LREE/HREE ratios, unlike all the other rocks analyzed (e.g. La/Sm_CH = 1.82–2.09; La/Yb_CH 3.84–4.12). Both portions of the enclave from sill 219D show the most enriched REE concentrations of the dataset and Eu negative anomalies (not shown).

Fig. 4

Representative chondrite- and primitive mantle-normalized Rare-Earth Elements (a, b) and Incompatible Trace Elements (c, d) multi-element diagrams for the analyzed dolerites. The diagrams showing the full dataset are in Supplementary information 1 (Fig.S5a-c). The ITE compositional field of the sediments analyzed by Sibik et al. (2015) and the average composition of continental crust (Rudnick and Gao 2014) are shown in panel d) for comparison with the host rocks analyzed from borehole 260

Incompatible trace elements (ITE) normalized to the Primitive Mantle (Sun and McDonough 1989) show uniform patterns, with strong negative Nb, and mild negative Ti, Th and P anomalies, and with positive spikes in Pb (Pb data available only for few samples) and Sr (Figs. 4c–d; 7c–d; S5a–c), i.e. typical arc-like signatures shared by many continental flood basalts (Wang et al. 2016; Pearce et al. 2021). The REE-enriched gabbroic evolved samples show a negative Sr spike and lack a negative P anomaly (Fig. 4d). This is accompanied by petrographic evidence for late crystallization of apatite in these rocks, mostly within the evolved late-stage pockets, with phosphorous (incompatible) being enriched in the most evolved melts. The dolerites from the 260 sills have very uniform ITE patterns, but samples of sill 260B never show a positive Sr anomaly (Fig. 7d). Sill 277/A1 shows higher abundance of the more incompatible elements and lower concentration of the least incompatible ones, with a pronounced negative Ti amonaly (Fig. 4d). The coarse-grained enclaves from sill 219D show very enriched ITE patterns (particularly in Zr, Hf and U; Fig. 4e). The leucocratic enclave is depleted in the mobile large-ion lithophile elements Rb, Ba and K with respect to the melanocratic one, possibly as an effect of the observable hydrothermal alteration. The host rocks are variably enriched in Sr and U, and depleted in Nb and Ti. Sample 260/2200 is enriched in K, whereas the halite-dominated sediment is depleted in K. ITE are strongly depleted in our evaporitic host rocks (Fig. 4e), as opposed to terrigenous lithologies (Pang et al. 2013; Sibik et al. 2015).

Isotopes

Most dolerites form a cluster with ⁸⁷Sr/⁸⁶Sr_i between 0.7048 and 0.7058 and ¹⁴³Nd/¹⁴⁴Nd_i between 0.51242 and 0.51235 (Fig. 5), overlapping the compositional fields for Morongovsky-to-Samoedsky Formations (Upper sequence) and for dolerites and differentiated intrusions (Naldrett et al. 1992; Hawkesworth et al. 1995; Fedorenko et al. 1996). Samples from 260A show the most uncontaminated compositions (0.70455; 0.51244). Dolerites from sill 219B, the shallowest from borehole 219, are shifted to the lower ¹⁴³Nd/¹⁴⁴Nd_i (0.51234) and ⁸⁷Sr/⁸⁶Sr_i (0.7046) end of the cluster, overlapping with the Dolerites and differentiated intrusions field (type 3 from Naldrett et al. 1992). Sill 277A1 is isotopically very different, showing high radiogenic Sr coupled with low radiogenic Nd (277/300 = 0.70797–0.512024 and 277/325 = 0.708214–0.512048), similar in composition to the intrusions of Lower Talnakh-type, and approaching compositions of the Nadezhdinsky lava Fm. (Fig. 5a; Naldrett et al. 1992; Hawkesworth et al. 1995; Yao and Mungall 2021).

Fig. 5

a) Sr–Nd isotope values and b) ⁸⁷Sr/⁸⁶Sr vs. Sr (ppm) of a subset of dolerites (n = 32) plus the leucocratic enclave 219/20070B. Isotopic values are corrected for radiogenic in-growth to 251 Ma, using trace element concentrations obtained by ICP-MS on the same aliquot sample. Compositional fields for the Nadezhdinsky (Nd) and Upper Sequence lava formations (Morongovsky to Samoedsky, Mr–Sm) and for the Norilsk and Lower Talnakh-type intrusions are from Naldrett et al. (1992), Hawkesworth et al. (1995), Fedorenko et al. (1996), and Sharma (1997). Previously published data from the cratonic STLIP are plotted for comparison (intrusive rocks from: Hawkesworth et al. 1995; Arndt et al. 2003; Ivanov et al. 2008; Kiselev et al. 2014; Krivolutskaya et al. 2020; Yao and Mungall 2021, and references therein; effusive rocks from: Lightfoot et al. 1990; Sharma et al. 1991; Sharma 1997; Wooden et al. 1993; Fedorenko et al. 1996; Krivolutskaya et al. 2020). Composition of intrusive and effusive rocks from Taimyr are also plotted (Reichow et al. 2016). Fields for OIB and MORB compositions are compiled from the GEOROC geochemical database (georoc.mpch-mainz.gwdg.de/georoc). ECAFC trajectories are calculated for i) upper crustal basement granitoids (gneiss and leucogranites; ⁸⁷Sr/⁸⁶Sr_i = 0.734; Sr = 280 ppm; ¹⁴³Nd/¹⁴⁴Nd_i = 0.5115; Nd = 16 ppm; Czamanske et al. 2000). Tick marks every 1%; ii) upper crustal metapelites (⁸⁷Sr/⁸⁶Sr_i = 0.7091; Sr = 233 ppm; ¹⁴³Nd/¹⁴⁴Nd_i = 0.5123; Nd = 26 ppm; Czamanske et al. 2000; ECAFC path reaches 2% assimilation with 60% crystallization). Mixing line with marlstone (Pang et al. 2013) is calculated until 10% mixing, tick marks every 1%. The parameters used for the ECAFC model are listed in Supplementary information 2. The compositions of these and other considered assimilants are shown in Supplementary Fig. S10

We observe two potential trends in the Sr–Nd isotopic space. A steep trend towards low radiogenic ¹⁴³Nd/¹⁴⁴Nd_i and high radiogenic ⁸⁷Sr/⁸⁶Sr_i is indicated by 277A1 samples. Conversely, four dolerites (from sills 219C, 219D and 277E) and the leucocratic enclave are shifted from the main cluster towards more radiogenic Sr isotope values at similar ¹⁴³Nd/¹⁴⁴Nd_i, forming a subhorizontal trend in Sr–Nd space. Rocks of this latter trend overlap the compositions of Norilsk-type intrusive rocks (Fig. 5a; Hawkesworth et al. 1995; Arndt et al. 2003; Yao and Mungall 2021). Strontium concentration and ⁸⁷Sr/⁸⁶Sr signatures correlate positively (Fig. 5b). The data from single sills do not describe an isochron on a ⁸⁷Rb/⁸⁶Sr vs. ⁸⁷Sr/⁸⁶Sr plot, suggesting that the system was perturbed, potentially by late-stage hydrothermalism. No correlation is observed between ⁸⁷Sr/⁸⁶Sr and whole-rock MgO, nor between Nd isotope signature and concentration (Fig.S6). The dolostone sample 260–2200 and the halite-anhydrite sample 260–3001.2 are fairly radiogenic in ⁸⁷Sr/⁸⁶Sr_i (0.708755; 0.707792, respectively). The first is enriched in Sr (1024 ppm), the latter is Sr-poor (38 ppm). These host rocks were not analyzed for Nd isotopes due to their low Nd concentration (4 ppm).

Discussion

Geochemical correlations with STLIP chemostatigraphy

We compare our data on the TB sills to previously studied rocks from the STLIP, proposing correlations between sills, formations and intrusive magma types (Figs. 6, 7 and 8). The correlation between subvolcanic rocks and lavas can be hindered by accumulation or fractionation of mineral phases. The mineral-whole rock disequilibrium suggests some accumulation of mafic phases in the TB dolerites (Fig. 3). However, unlike what oberved for cumulitic rocks from layered mafic intrusions, the TB dolerites do not show extreme depletions in trace element concentrations, and in most cases their major element compositions do not appear controlled by accumulation of some mineral phase (cf. the compositions of rocks from the Freetown Layered Complex; Callegaro et al. 2017; cf. also STLIP sills in Egorova and Latypov 2013). This is also supported by similarities in major and trace element concentrations between microporphyritic (chilled margin) samples, and holocrystalline ones in homogeneous sills (e.g. 260B). We stress however that the proposed correlations are mostly based on isotopes and REE ratios, that are nearly insensitive to fractionation or accumulation of a basaltic mineral assemblage. For instance, 50% fractional crystallization (FC) of 5% olivine, 30% cpx, 25% plag from the most primitive of the dolerites at 1 kbar, anhydrous conditions calculated by the Magma Chamber Simulator (MCS; Bohrson et al. 2020; Heinonen et al. 2020), is negligible on the REE ratios used for correlations (FC path has the same dimension as the symbols used in Figs.6c and 7b).

Fig. 6

Classification diagrams for STLIP lava Sequences based on a Th/U and Gd/Yb ratios (Fedorenko et al. 1996; Ivanov et al. 2008); b ⁸⁷Sr/⁸⁶Sr_i vs. ΔTiO₂, defined as demarcation line between high- and low TiO2 (line equation is: TiO₂ = 3.45–0.0317 × Mg#, where Mg# is Mg/(Mg + 0.85Fe_tot); elements in atomic units; Ivanov and Balyshev (2005); c Gd/Yb vs. La/Sm (Lightfoot et al. 1993; Hawkesworth et al. 1995) and enlarged inset. The effect of 60% FC (mineral assemblage of 5% ol; 30% cpx; 25% plag) calculated at 1 kbar, anhydrous conditions by MCS (Bohrson et al. 2020; Heinonen et al. 2020) is negligible on REE ratios (calculated vector is smaller than symbols). Model mantle melting curves A, B and C are plotted in panel c

Fig. 7

a Within-sill fractional crystallization (FC) path (orange) calculated starting from sample 219/1410–15, was calculated with the MCS (Bohrson et al. 2020; Heinonen et al. 2020). The curve represents 60% fractional crystallization (5% ol; 30% cpx; 25% plag) from a basalt, shown also in Fig. 3 for major elements. The FC is calculated at 1 kbar at anhydrous conditions. Purple line starting from Tuklonsky picrite composition (Yao and Mungall 2021) marks the trajectory for FC of 16% olivine, 2% spinel and 7% clinopyroxene from a picritic melt, at 3 kbar, anhydrous conditions. Partition coefficients for the FC paths are from Bédard (2014). b Mantle melting curves (also plotted in panels a and c are calculated with the REEBOX PRO software (Brown and Lesher 2016). The full parameters used for the models are listed in Supplementary information 2. Tick marks are placed every 1% melting for curve C. Compositional fields for the STLIP lava formations in c are built from data from Lightfoot et al. (1990, 1993); Wooden et al. (1993); Hawkesworth et al. (1995)

Fig. 8

a REEBOX PRO curves are plotted on a Th/Yb vs Nb/Yb diagram. Mixing curves with the average UCC (Rudnick and Gao 2014) is also plotted. Tick marks every 5%. In b the dolerite compositions are plotted on the newly proposed LIP-printing discrimination diagram of Pearce et al. (2021). Both diagrams show the concurrence of a depleted (MORB-like) component, and an enriched (crustal-like) component in the pretrogenesis of these dolerites

None of the rocks in our dataset is alkaline or high-Ti (Fig. 3), confirming absence of alkaline STLIP products south of the Maymecha-Kotuy region and absence of high-Ti products south of the Kulyumber river (Krivolutskaya et al. 2020). A plot of ⁸⁷Sr/⁸⁶Sr_i vs. ΔTiO₂ (Ivanov and Balyshev 2005; Fig. 6b) helps discriminate between sequences. Most samples fall in the Upper sequence compositional field, with only 5 samples in the Middle sequence field. However, in a plot of Th/U vs. Gd/Yb (Fig. 6a), all the samples classify as Lower sequence, except for one dolerite from 277A1, showing affinity with the Middle sequence. REE ratios (Gd/Yb and Sm/Yb vs. La/Sm) and isotopes classify all three rocks from 277A1 within the Nadezhdinsky Fm. and Lower Talnakh compositional fields (Figs.5a; 6c; 7b). Dolerites from all the other sills overlap with the Upper Sequence lava formations, from Morongovsky to Samoedsky (Mr–Sm in the Figures). Sills 219C-E, 274C and 277B-E correlate more specifically with the Morongovsky-Mokulaevsky Fm. and the Norilsk 5A intrusives. We note that samples from sill 219D overlap with the compositional field of Tuklonsky (Tk) formation in Fig. 6c, but correlation with this formation is excluded by the Th/U vs. Gd/Yb diagram (Fig. 6a), and by isotope data (Fig. 5a).

Correlations with intrusive types hinge on Sr–Nd isotopes (note that no isotope analyses are available for 274 and 275 borheoles). Of the TB sills, 219B, 260A-B and 277B can be all classified as type 3 intrusives (dolerites and differentiated intrusions in Fig. 5a). Of these, 219B shows higher Sm/Yb, and lower La/Sm and ¹⁴³Nd/¹⁴⁴Nd_i at comparable ⁸⁷Sr/⁸⁶Sr_i with respect to the other dolerites from the same borehole. These signatures are typical of the Ambarninsky type II dolerites that Hawkesworth et al. (1995) correlated with the Kharaelakhsky Formation (Hr in Fig. 1c). Sill 277A1 can be classified as a 5B-Lower Talnakh-type intrusive (Fig. 5a-b). All the remaining sills, 219C–E, 277B to E, following the suborizonthal Sr-enrichment trend, are classified as Norilsk-type 5A intrusions.

We stress that these sills belong to the STLIP event beyond reasonable doubt (Vasiliev et al. 2000; Svensen et al. 2018; Burgess and Bowring 2015). Older mafic rocks, such as thoe of the Viluy province (Kiselev et al. 2014) or the volcanics found as xenoliths in the Maslov diatreme (Czamanske et al. 2000), show very different geochemistry (e.g. alkaline affinity, more fractionated La/Yb, high TiO₂). The TB dolerites plot within the field of the STLIP on the newly proposed classification diagram (LIP printing; Pearce et al. 2021; Fig. 8b). Some sills cropping out at the periphery of the TB, or along valleys cutting trough it, show U–Pb ages well within the STLIP event (Burgess and Bowring 2015). The only potential exception is sill 260B, correlated with the Usol’skii sill (Ivanov et al. 2008), interpreted to be the youngest in the area and dated by ⁴⁰Ar/³⁹Ar to the Early Triassic (ca. 244 Ma; Ivanov et al. 2013), although U–Pb dating never showed such young ages (Burgess and Bowring 2015).

Sill profiles and within-sill differentiation processes

Compositional variations along profiles from tholeiitic sills can offer information into syn- and post-emplacement differentiation processes. We organized samples according to their level in the sill to allow identification of compositional profiles for compatible elements or differentiation indexes (i.e., Mg#, CaO, Ni, and Cr). For the compositional profiles in Figs.S5–6 we refer to established nomenclature, i.e., I-, D-, S- and slash-shaped profiles (cf. Gibb and Henderson 1992; Latypov 2003; Galerne et al. 2010; Egorova and Latypov 2013), and those that did not fit in any established pattern are flagged as X-shaped profiles. The sampling density was not equal for each sill and in some case samples were not numerous enough to determine compositional profiles (e.g. for sills 277E, 274C and 275G). Also, since we use borehole samples, we only have access to one profile from each sill, and cannot exclude the existence of different types of compositional profiles along the same sill (Galerne et al. 2010). D-shaped profiles, with the least-differentiated rocks at the center and more evolved margins, are visible for sills 219B (4 samples over 215 m thickness; Fig.S7a-b) and 219C (5 samples over 148 m) and sill 277A1 (3 samples over min. 64 m). Sills 277D (4 samples over 193 m) and 219E (6 samples over 162 m) show a slash-shaped profile (Fig. S7c, d), i.e. gradually more differentiated rocks from base to top. Top samples from these sills are enriched in REE (Fig. 4a, b), with mildly negative Sr and Eu anomalies (0.80–0.90), likely reflecting crystallization from an evolved melt that fractionated plagioclase. A typical S profile (highest Mg# rocks at the top, followed by more evolved compositions) is shown by the thin sill 260B (7 samples over 97 m thickness; Fig. S7g, h). I-shaped profiles (uniform composition throughout the thickness) are visible for sills 277B (3 samples over 85 m) and 277C (3 samples over 85 m), although samples are available only from the upper half of these thin sills. Complex profiles (X-shaped) are shown by sills 219D (14 samples over 148 m; Fig. S7e, f), 260A (8 samples over 29 m) and sill 277A2 (6 samples over 100 m). The slash-shaped profiles in sill 277D (193 m) and 219E (162 m) are typical of sills that developed cumulates from the bottom upwards, as a result of closed-system crystallization and gravity-driven crystal settling followed by compaction of the cumulate pile and subsequent extraction of the interstitial melt to upper levels. This is evidenced also by the REE and IE patterns of the samples from these sills (Fig. 4). For example, the composition of the evolved samples of sill 277D can be reached after 60% fractional crystallization in situ (5% ol; 30% cpx; 25% plag) from a basalt (Fig. 3; liquid lines of descent modeled by MCS, starting from the composition of the most primitive dolerite in our dataset 219/1410-15). Incompatible trace elements are depleted in the bottom samples, and enriched at the top, as a result of the preferential partitioning of these elements into the melt during the early stages of crystallization. This is illustrated by La vs La/Yb (Fig. 7a), where the evolved roofs and cumulitic floors of the slash-shaped sills 219E and 277D-E are tied by horizontal FC paths (i.e. La/Yb is insensitive to FC) calculated with MCS. Coherently, the Eu/Eu* and Sr negative anomalies in the coarse-grained evolved rocks from the top portions of these sills reflect crystallization after plagioclase fractionation. In sill 219E (Fig.S8), olivine is present in the bottom samples (219/3930.6, 3933.5 and 3891.7) but not in the very coarse-grained, evolved top samples. Instead these show well-developed late-stage granophyric melt pockets (e.g. 219/3848.3). D-shaped profiles are visible only in relatively thick sills. These profiles may be the result of post-emplacement melt flow from the central part of the sill to the top and bottom portions, induced by thermal stresses (Aarnes et al. 2008) or by uneven sill shape (Galerne et al. 2010). Relatively evolved melt (enriched in Fe and incompatible elements) flows towards the bottom and top of the sill, leaving behind a depleted high-Mg# crystal mush rich in compatible elements. However, we do not have enough data to correlate this information with the textural analysis and mineral chemistry.

We may tentatively interpret the gabbroic enclave in sill 219D as an autolith (Figs. S2 and S8), a coarse-grained portion of the same magmatic body remobilized by convective flux (Sparks et al. 1984), activated during cooling, or after a new batch of hotter, more primitive melt entered the system. Alternatively, given its more evolved geochemistry, the enclave may be a pegmatitic segregation (e.g. Marsh 2002), albeit at an embryonic stage, although this would be at odds with the sharp contact with the fine-grained dolerite. Brecciated sedimentary (limestone/dolostones and halite/anhydrite) xenoliths of the host rock were also found within sill 219D (Fig. S2), and can be interpreted as rafts of host rocks enclosed in the magmatic body by stoping. Localized assimilation of these xenoliths may be responsible for the chaotic compositional profile (X-shape) seen for this sill (Fig. S7e, f), and for the enrichment in alkali shown by some samples of 219D (Fig. 3c).

We calculated by MCS the effect of fractional crystallization (FC) on the composition of 4 initial melts (219/1410–1415; 260/2225–2230; 277/2380–2385 and a Tuklonsky-like picrite; for clarity, only one is shown in Fig. 7a). Closed-system evolution is confirmed for some of the sills, with samples compositions aligned along calculated subhorizontal FC vectors on a plot of La/Yb vs. La concentration (Fig. 7). Modeled FC (5% ol; 30% cpx; 25% plag) describes a near-horizontal trend, unable to reproduce all the spread in La/Yb covered by the data. The differences in La/Yb can reflect different parental melts or crustal assimilation, as backed up also by positive trends between La/Yb or La/Sm and La. The high Th/Nb (from 0.16 to 0.51) compared to mantle-derived melts (e.g. MORB 0.053; Sun and McDonough 1989; cf. Figure 8a) shown by all the dolerites also suggests some crustal contribution.

In summary, closed-system differentiation processes can explain the within-sill variability in major and trace element concentrations. Instead, the differences in isotopes and REE ratios between different sills and in some case within the same sill need to be explained by other processes, such as different parental melts or open-system differentation. There is evidence for open-system evolution in some of the sills, e.g. (i) inverse zoning in some minerals suggesting injection of more primitive magma; (ii) autoliths suggesting turbulence created by injection of a new batch of magma; (iii) xenoliths within the sills (e.g. 219D); (iv) isotopic evidence of late-stage interaction with the evaporite-host rocks; (v) textures and mineralogical assemblages developed by magma-evaporite interaction at the contacts.

Constraints on the mantle source of the TB dolerites

We can try to deconvolve the effects of mantle melting, fractional crystallization (closed-system evolution) and crustal assimilation (open-system evolution) by analysis of the trace element and isotope features of the TB dolerites, although a precise reconstruction of their parental melts and of their mantle source is hindered by their fractionated (Mg# 63–34; Ni < 180 ppm) and contaminated (Th/Nb > 0.16) nature. Most TB dolerites geochemically resemble the lavas of the Upper sequence (Morongovsky to Samoedsky; hereinafter Mr–Sm) and the Norilsk-type 5A intrusive rocks. The dolerites from sill 277A1 are geochemically similar to the lavas of the Nadezhdinsky Fm and to the Lower Talnakh intrusives. Fedorenko et al. (1996) distinguished four primitive magma types in the classical Norilsk lava piles, among which the Tuklonsky picrites and the Upper sequence basalts. Due to their very uniform geochemistry and their positive εNd signatures, and the absence of picrites in this suite, they interpreted the Upper sequence (Mr–Sm) basalts as primitive mantle-derived melts. However, we favor the interpretation that the TB dolerites reflect various levels of differentiation from a picritic parental melt. A forward calculation using the MCS shows that it is possible to obtain the major element composition of the most primitive dolerites by fractionating 16% olivine, 2% spinel and 7% clinopyroxene from a picritic melt (Fedorenko et al. 1996). This is in agreement with the observation that Ni contents between 130 and 200 ppm can be derived by 15% FC of olivine from a mantle-derived primitive magma with 300–450 ppm Ni (Sharma 1997). Indeed, the positive correlation between MgO and Ni (not shown) suggests some control of olivine on these rocks. However, an attempt of reconstructing a primitive melt by addition of olivine to the TB dolerites with the PRIMELT3 MEGA software (Herzberg and Asimow 2015) failed because even the more primitive dolerites (MgO > 8 wt%) were flagged as having fractionated clinopyroxene. The TB dolerites are thus moderately evolved melts derived from a parental picrite, but both the trace element and the isotope chemistry discard the Tuklonsky picrite as a viable parental melt for most of the TB dolerites. Tuklonsky picrites have too high La/Yb ratios and are isotopically too enriched (Figs. 5a and 7a, b) to represent the parental melts of the Upper sequence-like TB dolerites. On the contrary, Tuklonsky picrites can be related to the 277A1 dolerites by crustal contamination (as suggested previously for Nd-like basalts by, e.g. Wooden et al. 1993; Fedorenko et al. 1996; Reichow et al. 2005; Yao and Mungall 2021).

None of the analyzed dolerites show the characteristic highly fractionated LREE/HREE (e.g. La/Yb > 6) or MREE/HREE (e.g. Gd/Yb > 2) indicative of low degree and deep (> 120 km) mantle melting shown by the Lower sequence high-Ti basalts and picrites (Sobolev et al. 2009a, 2011). All the analyzed rocks show typical arc-like signatures, such as high LILE concentration and negative Nb and Ti anomalies (Fig. 4), in stark contrast with melts typically derived from primitive mantle sources (e.g. OIBs). For instance, they show low Nb/Y and Zr/Y (Fig.S9), below the so-called (Icelandic) mantle plume array (Fitton et al. 1997), overlapping with depleted, MORB-like compositions or with continental flood basalts from the Nb-depleted southern Karoo LIP (high percentage melts of a subduction-modified mantle source; Luttinen 2018; Fig.S9). We stress that Nb depletion can reflect either a Nb-depleted mantle source or crustal contamination (Luttinen 2018). Indeed, these arc-like features were observed before in the Upper series lavas, and were interpreted in contrasting ways. Some authors reconized this as a source-related character (recycling of crust in the mantle source, Lightfoot et al. 1993; melting of the lithosphere, Hawkesworth et al. 1995; Fedorenko et al. 1996; melting of a hydrous, subduction-modified mantle source; Ivanov 2007; Ivanov et al. 2008; Wang et al. 2016). Others interpreted this character as reflecting evolution of mantle-derived melts in long-staging reservoirs at the base of the crust, where constant replenishment and assimilation could occur (Wooden et al. 1993; Arndt et al. 2003). More recently, Sobolev et al. (2011) partially reconciled these views by suggesting a) that the source of the main tholeiitic phase magmas of the STLIP contained up to 50 wt% pyroxenite derived from reaction of recycled oceanic crust in the mantle plume with the peridotitic lithosphere, b) that melting occurred at shallow levels (< 60 km), and c) that the resulting melts were strongly contaminated by the continental crust. We try to evaluate these models based on the evidence offered by the TB dolerites.

The most uncontaminated dolerites in our dataset show higher ¹⁴³Nd/¹⁴⁴Nd_i than the primitive mantle recalculated to 250 Ma (Reichow et al. 2016), arguing for a more depleted asthenospheric component (DMM-like) in their mantle source. This is well represented by the signatures of the Lower Sequence basalts (⁸⁷Sr/⁸⁶Sr_i 0.7028; ¹⁴³Nd/¹⁴⁴Nd_i 0.512755; Sharma 1997). Conversely, the relatively radiogenic ⁸⁷Sr/⁸⁶Sr_i signatures point toward an enriched (crustal-like) component.

Forward partial melting models were calculated with REEBOX PRO (Brown and Lesher 2016). We simulated melting of a mantle source with mixed peridotitic and pyroxenitic lithologies (Sobolev et al. 2011) where the pyroxenites represent recycled (continental and oceanic) crustal domains in the mantle. Model A (Fig. 7a–c): melting of a source with multiple lithologies, i.e. 68% peridotite with trace element composition of the depleted MORB mantle, DMM (Workman and Hart 2005); 2% pyroxenitic component representing recycled continental upper crust (MIX1G; trace element compositions of UCC as in Rudnick and Gao 2014); 30% harzburgite that does not reach melting, representing a lithological scenario compatible with the thermal destabilization of the cratonic lithosphere below Siberia described by Sobolev et al. (2011). Decompressional melting of this mantle source was modeled considering a pre-existing lithospheric thickness of 50 km (Sobolev et al. 2011) and mantle potential temperatures (Tp) from 1450 to 1600 °C. We set a columnar shape for the melting region because more suitable for axial upwelling in mantle plumes or the rupture or foundering of thick, continental lithosphere (Tegner et al. 2019). Most dolerites from 219 and 277 boreholes are best represented by melting at the highest Tp, with a final pooled melt composed by ca. 10% pyroxenitic (the pyroxenite melts completely) and ca. 90% peridotitic melt (the peridotite melts up to 27%). The last melting occurs at ca. 2GPa, within the spinel stability field. This melting model results in a picrite that already bears the characteristic arc-like trace element pattern shown by the dolerites (Nb–Ta and Ti through, Pb spike), except for Rb, Ba and Sr (Fig. 7c). Assimilation of the TB sediments or circulation of late-stage fluids can reconcile these discrepancies in mobile elements. We note that similar results (not shown) are obtained by melting in the same T–P conditions a lithologically homogeneous mantle source (100% peridotite) in which 2–3% upper crust is recycled and well homogeneized (following the approach suggested by Tegner et al. 2019 for mantle source modeling of LIP-related tholeiites showing arc-signatures). However, dolerites from the borehole 260 show too low La/Sm ratios to be produced by this melting scenario. More importantly, recycling of upper continental crust in the mantle source is discarded on isotopic basis. We cannot know the precise isotopic composition of the recycled crustal component, but even recycling a moderately radiogenic upper continental crust (e.g. ⁸⁷Sr/⁸⁶Sr_i 0.710; ¹⁴³Nd/¹⁴⁴Nd_i 0.5116), the mantle melts obtained are too isotopically enriched to represent the TB dolerites. If there is a recycled crustal component in the mantle source, this could have been a lower crustal component with low radiogenic signatures (akin to some lower crustal granulite xenoliths found in Udachnaya; e.g. ⁸⁷Sr/⁸⁶Sr_i 0.70570; Sr = 235 ppm; ¹⁴³Nd/¹⁴⁴Nd_i 0.5120; Nd = 14.3 ppm; Koreshkova et al. 2011).

Model B: A likely parental melt for most TB dolerites is obtained by melting similar lithologies (Tp of 1500 °C, pre-existing lithosphere of 60 km) in which the pyroxenitic component has lower crustal characteristics. In this case the melting path ends at 2.08 GPa, producing a pooled picritic melt composed by 91% peridotite-derived melt (maximum melting degree of 18%) and by 9% pyroxenite-derived melt (maximum melting degree of 63%), and a relatively low radiogenic isotopic signature (⁸⁷Sr/⁸⁶Sr_i 0.7041; ¹⁴³Nd/¹⁴⁴Nd_i 0.51242). The fraction of pyroxenite-derived melt agrees with independent estimations based on trace element in olivine in Upper sequence lavas (max 50%; Sobolev et al. 2011).

Model C: A plausible parental melt for TB dolerites from 260 borehole is produced by 11% melting of a uniform peridotitic mantle source (DMM trace elements; Workman and Hart 2005; Tp 1500 °C; preexisting lithosphere 80 km) enriched by 2% recycled lower crust. The produced melt has isotopic signatures (⁸⁷Sr/⁸⁶Sr_i 0.7039; ¹⁴³Nd/¹⁴⁴Nd_i 0.51248) compatible with the 260A-B sills, and the slighlty thicker lithosphere can be explained by a farther south position of the borehole. We note that mafic lower crust (Archean tholeiitic protholiths granulitized in the Proterozoic, as shown by xenoliths sampled by the Udachnaya kimberlite), could have plausibly foundered or delaminated during early interaction of the STLIP plume whith the lithosphere. Alternatively, Czamanske et al. (2000) suggested that subduction/underthrusting of an accretionary block with mafic volcanics could have occurred shortly before the STLIP event under the western edge of the Siberian platform, and this material could enrich the Siberian sublithospheric mantle.

Caveat: this is a forward model, and the chosen parameters are plausible, but arbitrary, i.e. different combinations of parameters might lead to similar results, and we do not have strong constraints on the isotopic flavors of the recycled components. Furthermore, as exemplified in Fig. 8a (Th/Yb vs. Nb/Yb), the dolerite trace element patterns can be interpreted as reflecting crustal enrichment in the mantle source (cf. mantle melting curves obtained by REEBOX PRO modeling), or crustal assimilation by parental melt produced by a depleted mantle (cf. mixing curves with UCC and local gneissic and leucogranitic Proterozoic basement; Czamanske et al. 2000), i.e. it is difficult to deconvolve the effects of mantle vs. crustal enrichment. However, in summary, considering these dolerites as parental mantle melts from a sublithosheric mantle enriched with recycled upper crust (model A), as suggested by Fedorenko et al. (1996) or Lightfoot et al. (1993) for the correlative Upper Sequence lavas would result in the best match for the trace element patterns, but would require (a) probably unreasonably isotopic depleted signatures for the recycled upper crustal component and (b) unreasonably low-MgO contents for picrititc melts produced at the modeled Tp and P conditions. The TB dolerites have fractionated olivine and clinopyroxene, as seen from attempts to reconstruct their primitive melt with the PRIMELT3 software and from MCS calculations. Rather, TB dolerites were more plausibly derived from picritic parental melts that fractionated to basalts by ponding in crustal magma chambers. It is also likely that they assimilated some crust in this context (see below), unless they fractionated below the lower crust (cf. Heinonen et al. 2016). Estimated pressures of crystallization of clinopyroxenes for Middle and Upper sequence rocks reveal crystallization at < 20 km depth, i.e. within the upper crust (Yao and Mungall 2021; Suvorov et al. 2006). Clinopyroxenes from sill 260B are an exception, as they reveal crystallization pressures up to 8 kbar (Neave and Putirka 2017). We hence favor a model that explains the arc-like signatures of the TB dolerites partly as inherited from recycled lower crust (Model B–C) present in the mantle source in the form of pyroxenites (Sobolev et al. 2011), and partly as acquired by assimilation of the crust (Wooden et al. 1993; Arndt et al. 2003). Alternatively, the parental melts of the TB dolerites could derive from a depleted asthenosphere (DMM-like, with the isotopic flavor of the Lower Sequence rocks; Sharma 1997; Model D), acquiring their crustal-like trace elements only by assimilation en-route to the surface. We find this interpretation less likely because (a) it would not explain the high percentage of pyroxenite-derived melts calculated by Sobolev et al. (2011); (b) in such case one would expect to see STLIP products with compositions in between the more isotopically depleted Lower Sequence and the more enriched Middle and Upper Sequences; and (c) this scenario would require higher percentages of assimilated crust, which are not favored by the moderately radiogenic ¹⁸⁷Os/¹⁸⁸Os signatures found in the STLIP main tholeiitic series (Horan et al. 1995). We also note from modeling attempts that adding even small amounts of water to the mantle source (100–200 ppm) results in produced melt compositions that do not fit our data (e.g. more fractionated MREE/HREE), hence the common presence of biotite in these rocks is more indicative of water input from interaction with the sedimentary host rocks, such as marlstones (Pang et al. 2013), rather than a hydrous mantle source (cf. Ivanov et al. 2008; Wang et al. 2016; Ware et al. 2018).

Crustal interaction: two different processes, two diverging isotopic trends

Trace element and isotopic geochemistry of sills from the Lower Tunguska show that they experienced some degree of crustal contamination. In Sr–Nd isotopic space, two trends are suggested by some of our data, and they are further supported by previously published data (Fig. 5a). Some samples from borehole 219 and 277 describe a trend towards higher radiogenic ⁸⁷Sr/⁸⁶Sr_i at nearly constant ¹⁴³Nd/¹⁴⁴Nd_i, with respect to the main data cluster. The leucocratic enclave 219D/2007.0B plots to the Sr-radiogenic end of this trend (Fig. 5a). Norilsk-type STLIP intrusives (Naldrett et al. 1992; Hawkesworth et al. 1995; Arndt et al. 2003; Yao and Mungall 2021), some lavas of the Upper Sequence (from Morongovsky to Samoedsky Formations), and one of the Kulyumber dolerites (Krivolutskaya et al. 2020) also plot along this trend. Conversely, no rocks from Taimyr (Reichow et al. 2016) follow it.

A different trend, towards high radiogenic ⁸⁷Sr/⁸⁶Sr_i and low radiogenic ¹⁴³Nd/¹⁴⁴Nd_i, is indicated by the dolerites from 277A1 (Fig. 5a). Lower Talnakh-type intrusions, lavas of the Nadezhdinsky Fm. and some dolerites from the Taimyr region show similar isotopic compositions (Hawkesworth et al. 1995; Arndt et al. 2003; Reichow et al. 2016; Yao and Mungall 2021). A few data from the literature plot in an intermediate position between more isotopically depleted samples and the 277A1 dolerites, allowing for complete “sampling” of the trend (Hawkesworth et al. 1995; Kiselev et al. 2014; Reichow et al. 2016). We note that intrusive rocks from the Kulyumber river area, south of Norilsk (Figs.1a and 5a; Krivolutskaya et al. 2020) also fall on a steep trend in Sr–Nd space, but with a different slope with respect to the one shown by the dolerites studied here. These two isotopic trends suggest involvement of at least two different contaminants and/or processes, one that is able to change both Sr and Nd isotopic signatures (steep trend), and one that is actively swaying Sr isotopes and concentrations, but is ineffective on the Nd isotopic system (sub-horizontal trend).

Contamination with granitoids and metapelites in upper crustal magma chambers

Forward mantle melting modeling reveals that part of the enriched trace element and isotopic signatures in the TB dolerites can be imparted directly by their mantle source, which included recycled lower crust. However, assimilation in crustal reservoirs appears to play an additional role, a scenario favored also by previous authors (Wooden et al. 1993; Arndt et al. 2003; Reichow et al. 2005, 2016; Yao and Mungall 2021). Specifically, addition of radiogenic Sr and LILE to the modeled parental melts is required. We start from the composition of the melts obtained by the REEBOX PRO calculations (Models B and C), and estimate crustal assimilation by Energy Constrained Assimilation and Fractional Crystallisation modeling (EC-AFC; Spera and Bohrson 2001; Heinonen et al. 2021; full parameters are given in the Supplementary information 2). As evidenced from geobarometry on clinopyroxenes in correlative STLIP lava flows of the Upper and Middle sequence (Yao and Mungall 2021), the magma reservoirs were located in the upper crust (5–20 km depth; Suvorov et al. 2006). The basement below the Tunguska Basin is sampled by the Maslov diatreme, 25 km south of Norilsk (Czamanske et al. 2000), which transported xenoliths of Proterozoic basement granitoids (gneiss, leucogranites, minor trondjemites) and metapelites. There is variability in isotopic compositions for the silicic rocks, but they all show very high radiogenic ⁸⁷Sr/⁸⁶Sr_251Ma (> 0.720) and low radiogenic ¹⁴³Nd/¹⁴⁴Nd_251Ma (< 0.5118). We modeled contamination using available compositions of gneiss and leucogranites, and metapelites (Czamanske et al. 2000). This composition is reasonably representative of the Siberian craton upper crust, given that similar lithologies crop out also further SE in the Aldan shield where the Archean and Proterozoic basement is exposed (Jahn et al. 1998). Very minor assimilation of this upper crustal lithologies allows reproducing the main cluster of TB dolerites isotopic signatures (Fig. 5a). In particular, 1% assimilation of granitoids with 50% crystallization, can reproduce the isotopic signature of sill 219B, whereas 2% assimilation of metapelites with 60% crystallization can reproduce the Nd systematics in the dolerites (Fig.S6). The evolved isotopic signatures of 277A1 rocks (i.e. low radiogenic ¹⁴³Nd/¹⁴⁴Nd_i and relatively high radiogenic ⁸⁷Sr/⁸⁶Sr_i) are reproduced by 5% assimilation of granitoids with 65% crystallization. The same signatures are seen in the Nadezhdinsky formation (Fedorenko et al. 1996), in the Lower Talnakh group intrusives (Hawkesworth et al. 1995; Arndt et al. 2003), in the West Siberian Basin basalts (Reichow et al. 2005) and in Taimyr (Reichow et al. 2016). Different plausible contaminants were suggested to explain these signatures, e.g. tonalitic upper crust (Lightfoot et al. 1990); the Bolgokhtokh granodiorite (Lightfoot et al. 1993; Wooden et al. 1993; Hawkesworth et al. 1995; Sibik et al. 2015); alkaline or acidic STLIP products (Reichow et al. 2005, 2016); granitoid crust (Arndt et al. 2003; Yao and Mungall 2021). In all these previous works (but also in Naldrett et al. 1992; Lightfoot and Keays 2005; Sobolev et al. 2009a), the Tuklonsky picrites were considered as the parental melts of the enriched Nadezhdinsky Fm. and Lower Talnakh-type intrusions. These picrites show enriched isotopic signatures and trace elements (¹⁴³Nd/¹⁴⁴Nd_i = 0.5122 and ⁸⁷Sr/⁸⁶Sr_i = 0.7057; Nd = 7 ppm; Sr = 175 ppm; Sobolev et al. 2011; Yao and Mungall 2021), probably reflecting the presence of recycled sediments in the lithosphere (Lightfoot et al. 1993; Hawkesworth et al. 1995). Starting from the Tuklonsky picrite ca. 4% assimilation of Proterozoic granitoids after 60% crystallization of the initial Tuklonsky picrite can explain the isotopic signatures observed (Fig. 5a), and the Nd concentrations in 277A1 rocks (Fig.S6). We do not have a strong preference for one starting picrite against the other.

What is evident from Fig. 5b is a decoupling between Nd and Sr in the TB dolerites. This is best explained by assimilation of Sr-rich lithologies, that are abundant in the Tunguska Basin sediments (marlstones, limestones, anhydrites; Pang et al. 2013), and is discussed in the following section.

From ECAFC model attempts, the lower crust (composition based on lower crustal granulite xenoliths from Udachnaya; Khoreshkova et al. 2011) appears as a less likely contaminant for the petrogenesis of most TB dolerites. Based on geobarometry, dolerites from sill 260B originated from melts that ponded in lower crustal chambers, but ECAFC modeling (starting from model C) discards their assimilation of lower crust. The overall picture emerging from assimilation modeling is that the plumbing system below the TB was complex and involved multiple magma chambers at different crustal levels, where the mantle-derived melts ponded and evolved before their emplacement in the basin.

Contamination in the volcanic basin: magma-sediment interaction involving carbonates and evaporites

Contamination trends involving Sr but not Nd such as those depicted by some of our dolerites were observed before in several studies of Siberian Traps intrusive rocks, and were interpreted as the result of contamination occurring at very shallow crustal levels, during emplacement of the sills in sedimentary host rocks (Hawkesworth et al. 1995; Arndt et al. 2003; Ripley et al. 2003, 2010; Lightfoot and Keays 2005; Keays and Lightfoot 2010; Yao and Mungall 2021). We also note that the subhorizontal Sr-trend is described by samples coming from sills and intrusions intruded in evaporitic or carbonatic host rocks (our sills 219C–E, 277E; intrusive rocks described in Hawkesworth et al. 1995; Arndt et al. 2003; Yao and Mungall 2021). Conversely, no rocks from Taimyr, where evaporitic deposits are absent, follow this trend (Reichow et al. 2016). All these observations support an origin of these signatures by interaction with the TB sediments.

Energy-constrained models involving assimilation of non-silicate lithologies are not available to date due to a lack of thermodynamic data (Heinonen et al. 2021), hence we evaluate potential assimilation paths by mixing models. In Fig. 5b, all the plotted compositions seem to converge towards a ⁸⁷Sr/⁸⁶Sr_i of 0.708–0.709. All our analyzed host rocks have isotpic signatures within this range (cf. Supplementary information 2). Most metapelites, skarn, anhydrite and marlstones analyzed by Pang et al. (2013) also show similar Sr isotopic values, and when analyzed, relatively radiogenic ¹⁴³Nd/¹⁴⁴Nd_i (0.5124–0.5122), comparable with those of the dolerites. We did not measure ¹⁴³Nd/¹⁴⁴Nd_i in the host rocks, but the concentration of Nd in evaporitic lithologies is much lower than that of a basaltic melt (see Fig. 4d), and because of this, assimilation of these lithologies, whatever their ¹⁴³Nd/¹⁴⁴Nd ratio, would be unable to change the Nd isotopic signature of the initial melt significantly. As a result, a contamination trend that affects only Sr is produced. In support of this, while there is a rough positive correlation between Sr and ⁸⁷Sr/⁸⁶Sr in the dataset (Fig. 5b), a plot of Nd vs. ¹⁴³Nd/¹⁴⁴Nd does not show any correlation (except for sill 277A1 rocks; Fig. S6).

Mixing lines can be calculated with different starting compositions, as illustrated in Fig. 5b. Up to 7% uptake of marlstone can account for the Sr enrichment. From mass-balance calculations, addition of 7% marlstone does not significantly change the major element composition of the magma. On the contrary, similar amounts of assimilated dolostone or anhydrite are expected to concomitantly enrich the magmas in Ca. Assimilation via melting or chemical dissolution of anhydrite would also enrich them in S, not only in Sr (Iacono-Marziano et al. 2017). For instance, Sr-rich dykes of the Central Atlantic Magmatic Province in Portugal, whose compositions were modeled with a simple mixing of up to 20% limestone—80% basalt, result in strong Ca and Sr enrichment in the most contaminated rocks (Callegaro et al. 2014). Conversely, the TB dolerites show no Ca enrichment and Sr does not reach 10³ ppm levels. During assimilation of carbonate rocks, whole-rock Ca can be buffered by FC of Ca-rich minerals (e.g. augite or plagioclase; Whitley et al. 2020), and might be recorded by enrichment in Ca in the whole rocks or phenocrysts (e.g. Wenzel et al. 2001, high Ca in olivine). We do not see these effects in the mineral chemistry. Furthermore, there is no evidence of melting of the country rocks at the contacts, at least for the sills sampled. Black et al. (2014) and Yao and Mungall (2021) calculated evaporite assimilation around 0.08–0.25 wt% from isotopic analyses of sulfur on STLIP intrusive rocks. Bulk assimilation of carbonatic and evaporitic lithologies might indeed be limited to these very low amounts also in the TB dolerites. Rather, there is evidence for metasomatic late-stage interaction with the host rocks, such as (a) the presence of hydrothermal veins (with carbonates, sulfates, scapolite) cutting through the dolerites, (b) enrichment of mobile elements derived from the host rocks (e.g. Sr, Ba, Na, K; Fig. 9 and S11); (c) petrographic and geochemical evidence of hydrothermal alteration in the leucocratic enclave from sill 219D (219/2007.0B; Figs.S2; S3; S8), which plots to the radiogenic end of the Sr-trend; (d) presence of small (up to 500 µm) evolved interstitial pockets with biotite and euhedral quartz, texturally suggesting a hydrothermal origin. Several studies have shown that interaction between magma and evaporitic/carbonatic host rocks can take place by various processes besides direct melting or bulk assimilation, e.g. brine-borne diffusive transfer, hydrothermalism, or chemical dissolution (Ripley et al. 2003; Li et al. 2009; Deegan et al. 2010; Pang et al. 2013; Schofield et al. 2014; Mata et al. 2015; Iacono-Marziano et al. 2017; Yudovskaya et al. 2018; Heimdal et al. 2019).

Fig. 9

Schematic N–S cross-section across the Tunguska basin illustrates the petrogenesis of the TB dolerites (horizontal and vertical scales non homogeneous). Melting of a mixed peridotitic and pyroxenitic source with recycled lower crustal domains (1) produced picrites that assimilated small amounts of granitoids and metapelites while ponding in upper crustal magma chambers (2). Further interaction with host rocks (marlstones and evaporites) after emplacement within the sedimentary basin enriched the dolerites in Sr, Cl, Na (3). The sills are color-coded suggesting the proposed correlations with the Norilsk chemostratigraphy

In-situ interaction with halite-dominated host rocks deserves separate attention. Assimilation of halite is not very effective in modifying Sr isotopic signatures of the dolerites, but it can strongly affect Na, K and Cl concentrations. Dolerites from borehole 260 and 275, and some samples from sill 219D show a trend of increasing alkali contents with decreasing silica in the TAS diagram (Fig. 3c). This is opposite to the positive trend shown by the other dolerites, reflecting evolution by fractional crystallization and crustal assimilation. We calculated simple mixing models between the more primitive, less alkali-rich 260 dolerite and halite. The trend shown by the samples is best reproduced by mixing of up to 7 wt.% with an impure salt composition (Pang et al. 2013), containing high alkali (31 wt.%), but also some silica (6 wt.%). We do not have data on whole-rock Cl concentrations, but preliminary element mapping on biotites from evolved interstitial pockets in the TB sills revealed high concentrations of Cl at the crystal rims, similar to Heimdal et al. (2019). This potentially reflects formation of a Cl-rich fluid and elemental transfer via a hydrous fluid, as observed by Pang et al. (2013) for aureoles surrounding mafic sills at Norilsk. Ivanov et al. (2008) put forward that derivation from an arc mantle source would explain many geochemical features observed for the STLIP, for example the Sr positive anomaly and the presence of primary micas in the basalts. Alternatively, Sr enrichment can be explained by crustal assimilation (e.g., Wooden et al. 1993). We note that the modeled trace element compositions of the parental melts produced by a mantle source enriched with recycled crust (Models B–C) lack the Sr and K positive anomalies observed for most TB dolerites. We have shown how assimilation of marlstone and late-stage interaction with evaporitic lithologies can enrich Sr and K in the basalts, hence Sr concentration and the presence of biotite are not good proxies of mantle source composition in the STLIP.

The TB dolerites show evidence of contamination with marlstone and additional interaction with evaporitic host rocks, potentially via circulation of fluids produced within the metamorphic aureoles around the sills. We find evidence of this signature in sills from all the investigated boreholes, located across all the Tunguska Basin. These new data are the first to demonstrate that interaction with evaporites did not occur exclusively around sills (Pang et al. 2013) and large intrusions (Naldrett et al. 1992; Arndt et al. 2003; Ripley et al. 2010; Yao and Mungall 2021) of the Norilsk area, but was a widespread process in the STLIP plumbing system, as previously suggested by Svensen et al. (2009, 2018).

Relative chronology and potential implications

The shallowest sill analyzed, 277A1, located in the most central area of the STLIP, is chemostratigraphically the oldest of all the considered sills, and 260B the youngest (potentially younger than the STLIP main phase; Ivanov et al. 2013). All the other sills are correlated with the Upper Sequence activity of the STLIP, and in particular, sills 219C–E, 274C and 277B–E, correlate best with the Morongovsky-Mokulaevsky Fm. and the Norilsk 5A intrusives. These sills represent not only the most voluminous magma-type in the analyzed boreholes (minumum 811 m of cumulative thickness), but also intrude exclusively evaporites (rock salt, carbonates and marls), and show clear signs of interaction with these host rocks. Svensen et al. (2018) calculated through thermal modeling that the contact metamorphic aureoles produced by cooling of the sills can generate 52–80 tons of CO₂ m⁻². On a basin-wide perspective, this means that only 0.7–1.2% of the total area in the Tunguska Basin (i.e. an area of 12–19000 km² out of a region in excess of 1.5 × 10⁶ km² encompassed by the sill complex) is sufficient to generate 1000 Gt CO₂. The observed shallow level magma-sediment interaction with limestones and evaporites supports these calculations with chemical and petrographic evidence. It should be noted that widespread occurrence of explosive vents rooted at the sills and tephra troughout the Tunguska basin (Svensen et al. 2009, 2018; Jerram et al. 2016; Polozov et al. 2016) supports the notion of discharge of thermogenic gases following sill emplacemement and associated magma-sediment interaction.

Based on magnetostratigraphic evidence, Pavlov et al. (2019) proposed that the Norilsk-type 5A intrusions are coeval or nearly coeval with the boundary between the Morongovsky and Mokulaevsky Fm. In turn, some of the Norilsk-type intrusions overlap in time with the extinction interval (Burgess and Bowring 2015). Because a significant portion of the TB sills correlate geochemically with Mr-Mk Formations and the Norilsk-type 5A intrusions, it is intriguing to think that they might also be coeval. This hyporthesis could be tested in future works by precise U–Pb dating. If confirmed, it might mean that this brief, voluminous pulse of magmatism was not restricted to the Norilsk area, but can be tracked all across the central Tunguska Basin, where there is widespread evidence for magma–sediment interaction and high potential for thermogenic gas production.

According to Burgess et al. (2017), there is a shift in style of STLIP activity from effusive (phase 1) to intrusive (phase 2), and the beginning of phase 2 coincides in time with the End-Permian extinction as dated at the GSSP in Meishan, China (between 251.941 ± 0.037 and 251.880 ± 0.031 Ma; Burgess et al. 2014). However, most dated sills in the TB postdate the extinction interval, and the extinction-triggering portion of the STLIP must be a currently unidentified intrusive subinterval. Three U–Pb zircon ages were published by Burgess and Bowring (2015) for the Norilsk intrusion (G22-105–2: 251.64 ± 0.10/0.12/0.30 Ma; G22-63–5: 251.907 ± 0.067/0.094/0.29 Ma; G22-65–0: 251.813 ± 0.065/0.092/0.28 Ma). The oldest age overlaps the extinction interval. From the geochemical correlations presented here, many of the sills buried below the Tunguska basin could be good candidates for the unidentified intrusive subinterval called for by Burgess et al. (2017), responsible for outstanding thermogenic gas discharge to the end-Permian atmosphere. However, in absence of high-precision geochronology data, this link remains speculative, because geochemical correlation might not necessarily mean age concurrence. For instance, intrusions from Taimyr interpreted to have a Norilsk-type 5A signature, were dated by ⁴⁰Ar/³⁹Ar to ~ 247 Ma (Reichow et al. 2016), i.e. younger than the ages shown by the intrusions in the Norilsk region (Burgess et al. 2017). It should be noted however that this age has a very large uncertainty and that these rocks are geochemically different from the sills presented here that correlate with Norilsk-type 5A intrusives (e.g. in Zr/Y–Nb/Y; Fig.S9). We also stress that, (1) no U–Pb ages were obtained by Burgess and Bowring (2015) from the tholeiitic lava Formations in Norilsk; (2) none of the 17 ages published by Burgess and Bowring (2015) was produced from samples of the sills studied here, and (3) full geochemical data are not available for the dated sills, except for sill A10-23.1 from Bratsk, previously studied by Black et al. (2014), which is compositionally very different from all the sills included in our dataset (it correlates with the evolved portion of the Tolstomysovskii sill of Ivanov et al. 2008; Fig. 3). As put forward by Burgess and Bowring (2015), the beginning of intrusive activity throughout the STLIP is unconstrained, and future U–Pb geochronology endeavors should focus on the sills buried beneath the surface in the Tunguska Basin, despite their mafic composition. High-precision U–Pb data were successfully obtained on mafic sills in sedimentary basins from other LIPs (Greber et al. 2020; Davies et al. 2021), where zircon and baddeleyite crystallized in evolved interstitial pockets similar to those observed for the TB dolerites.

Concluding remarks

Deep sills intruded into the Tunguska Basin are low-Ti, tholeiitic dolerites geochemically similar to the Upper Series lavas of the Siberian Traps LIP. Sill parental melts are most likely picrites derived by extensive melting of a mantle source with recycled crustal components below a thinned lithosphere (80–50 km), within the spinel stability field. The recycled components, most likely foundered or subducted lower crustal domains, were potentially represented by pyroxenites. The parental melts underwent open-system evolution in upper custal magma chambers, assimilating from 1 to 5% granitoids and metapelites from the Proterozoic basement. Some sills were emplaced within the Tunguska Basin sedimentary sequence as a single pulse, followed by closed-system crystal fractionation, producing slash-shaped or I-shaped compositional profiles. Others show evidence for magma-mixing in their crystal cargo (inverse zoning) or turbulence caused by replenishment (autholiths), resulting in D- or complex-shaped profiles. Most sills show evidence for interaction with volatile-rich host rocks (marlstone, carbonate, anhydrite, and halite). Widespread occurrence of explosive pipes throughout the Tunguska basin suggest that this interaction led to thermogenic gas production and discharge. On the basis of isotope and trace element geochemistry we conclude that a) sill 277A1correlates with the Lower Talnakh-type 5B intrusions and the Nadezhdinsky lava formation; b) sill 219B correlates with the Kharaelakhsky (Hr) lava Formation and the Ambarninsky type II dolerites; c) sill 260B correlates with the Usol’skii sill, and is probably younger than all the lava formations defined for the Norilsk area; d) all the remaining sills correlate with Upper sequence lavas and, more specifically, sills 219C–E, 274C and 277B–E correlate with Morongovsky-Mokulaevsky lava Formations and with Norilsk-type 5A intrusions. This latter magma-type is recorded in at least three boreholes across the Tunguska Basins, is the most voluminous, and shows evidence for interaction with volatile-rich host rocks. Future high-precision geochronology will test the hypothesis presented here that these sills represent the hitherto unidentified early intrusive phase of STLIP responsible for extensive thermogenic gas production and discharge impacting the end-Permian atmosphere.