Introduction

Assimilation of sedimentary sulfur greatly improves the probability for the formation of economically important copper-nickel-platinum group element (Cu-Ni-PGE) sulfide deposits from high-temperature komatiitic and basaltic magmas (e.g., Keays and Lightfoot 2010; Ripley and Li 2013; Barnes et al. 2016). Yet, rocks formed from hot magmas with high assimilation potential can be sulfide-poor even though they show clear spatial association with comparatively sulfur-rich sedimentary rocks (Barnes et al. 2016). The predominant process of sulfur assimilation depends on the type of sedimentary rock and can vary as a function of temperature. Assimilation of anhydrite evaporite is relatively simple as it dissolves in magmas and can supply sulfur until the silicate melt becomes sulfate-saturated (e.g., Iacono-Marziano et al. 2017). In the case of carbonaceous material and sulfide-bearing siliciclastic sedimentary rocks, such as black shales, sulfur can be mobilized by devolatilization fluids (e.g., Ripley 1981; Thériault and Barnes 1998; Ripley et al. 2007; Virtanen et al. 2021) or it can be transported to the magma in the form of solid or molten sulfides when the surrounding silicate material undergoes partial or wholesale melting (e.g.Thériault and Barnes 1998; Queffurus and Barnes 2014; Samalens et al. 2017a; Virtanen et al. 2021). Some notable examples of sulfide deposits that formed due to the interaction between magma and black shale are found in Duluth Complex in Minnesota, USA (e.g., Ripley 1981; Thériault and Barnes 1998; Ripley et al. 2007), Pechenga in northwestern Russia (Barnes et al. 2001), and Kabanga in northwestern Tanzania (Maier and Barnes 2010). The nature of sulfur liberating processes in the siliciclastic rocks is complex and depends on the main sulfide species, the amount and type of carbonaceous material (CM), and the silicate assemblage, the latter of which controls the majority of the devolatilization and melting reactions. Constraining the phase assemblage at the time of magmatic assimilation is difficult due to the overprinting by contact-metamorphic and possible subsequent regional metamorphic reactions. Consequently, conditions prevailing in sulfide-bearing sedimentary rocks prior to the associated magmatic events are rarely constrained even though the likelihood of sulfur assimilation is dictated by these conditions. One notable exception is the Virginia Formation black shale adjacent to the Duluth Complex, which has been extensively characterized for its regional metamorphic and contact-metamorphic phases and reactions.

The Duluth Complex is a Mesoproterozoic composite mafic layered intrusion that hosts cumulative Cu-Ni(-PGE) sulfide resources of 4.4 Gt at 0.66% Cu and 0.22% Ni, which are some of the largest known on a global scale (Listerud and Meineke 1977; Miller et al. 2002). Based on the sulfur isotopic compositions, the adjacent sulfide-bearing Virginia Formation black shale is the dominant source of sulfur in the intrusion-hosted deposits, outweighing the contribution from magmatic sources (Mainwaring and Naldrett 1977; Ripley 1981; Arcuri et al. 1998). Troctolite-hosted disseminated Cu-Ni(-PGE) sulfide deposits are volumetrically the dominant type and typically present in the basal South Kawishiwi and Partridge River intrusions (Ripley 1981; Thériault and Barnes 1998; Ripley et al. 1999; Miller et al. 2002). Comparison of the well-studied regionally metamorphosed black shale precursor with the contact-metamorphosed counterpart (Ripley 1981; Rao and Ripley 1983; Ripley and Alawi 1988) together with heating experiments conducted with the black shale (Virtanen et al. 2021) indicate that sulfur was likely transported to these disseminated sulfide deposits by devolatilization fluids (see also Baker et al. 2001). The ability of the fluids to effectively mobilize sulfur was ascribed largely to the presence of reduced CM in the black shale (e.g., Ripley 1981; Virtanen et al. 2021). The CM is of biogenic origin and the sulfides, mostly in the form of disseminated pyrite, in the Virginia Formation have been thought to largely originate as a result of bacterial sulfate reduction (e.g., Arcuri et al. 1998).

Volumetrically minor and sparse massive sulfide lenses of the Duluth Complex are generally found within the noritic rocks that were subject to more thorough assimilation compared to the troctolites as evident from a clear change in the igneous mineral assemblage (olivine replaced by orthopyroxene) at the vicinity of abundant black shale xenoliths (Thériault and Barnes 1998; Miller et al. 2002; Queffurus and Barnes 2014; Samalens et al. 2017a). Importantly, many of these black shale xenoliths originate from a peculiar and volumetrically minor CM- and pyrrhotite-rich portion of the Virginia Formation called the Bedded Pyrrhotite Unit (BPU: e.g.Zanko et al. 1994; Queffurus and Barnes 2014; Samalens et al. 2017a). The BPU xenoliths acted as the main source for sulfur in the massive sulfide occurrences but here sulfur assimilation was mediated by partial melt rather than fluids (Queffurus and Barnes 2014; Samalens et al. 2017a). Compared to the normal black shale of the Virginia Formation, the precontact-metamorphic history of the BPU is poorly defined, especially for the CM and sulfur. Bacterial sulfate reduction has been suggested as the origin for sulfur based on the heavy δ34S values (Zanko et al. 1994), but it has not been studied in detail why the BPU differs from the rest of the Virginia Formation by having higher CM content and pyrrhotite as the main sulfide.

Here we characterize the structural state of the CM with Raman spectroscopy in different black shale units of the Virginia Formation. We use this data together with volatile chemistry and detailed petrographic observations to shed light on the formation of the CM- and sulfur-rich BPU. We also discuss how sulfide speciation affected sulfur mobilization during contact metamorphism and how retrograde hydration reactions contributed to the volatile contents of the contact-metamorphosed BPU. Results of this study contribute to the better understanding of sulfur assimilation processes during emplacement of mafic-ultramafic magmas and to refining of models for formation of magmatic sulfide ore deposits.

Geological setting

The Virginia Formation is a Paleoproterozoic thinly bedded black shale formation located in Minnesota, USA (Lucente and Morey 1983; Severson and Hauck 2008; Fig. 1). It is part of a continuous sequence of quartzites, iron-rich carbonates, and black shales that gradually formed during transition from a stable continental shelf to a deep basin-type depositional setting (Lucente and Morey 1983; Ojakangas et al. 2001). The stratigraphically upward change from the carbonate-dominated Biwabik Iron Formation (Fig. 1) to the organic carbon (Corg) bearing Virginia Formation black shale is gradational and marks a shift from shallower depositional regime to deeper anoxic environment (Ojakangas et al. 2001). The Virginia Formation is correlative with the black shales of the Rove Formation to northeast and the Thomson Formation to south (Hemming et al. 1995; Taranovic et al. 2015; Fig. 1). The Gunflint Formation that is known for well-preserved microfossils (e.g., Rasmussen et al. 2021) is stratigraphically below the Rove Formation and correlative with the Biwabik Iron Formation (Fig. 1). The current stratigraphic thickness of the Virginia Formation is at least 500 m and the whole formation is interrupted with NW-SE trending strike-slip faults caused by post-depositional crustal deformation (Lucente and Morey 1983). The Virginia Formation was subjected to subgreenschist to greenschist facies (up to ∼300–400 °C) regional metamorphism at 1786 Ma (Lucente and Morey 1983; Ripley et al. 2001a; Vallini et al. 2007). The regional metamorphic grade of the Virginia Formation increases from north to south (Hemming et al. 1995).

Fig. 1
figure 1

Geological sketch map showing the locations of the Virginia Formation and the Duluth Complex. The stars show the locations of sampled drill cores and one outcrop (DC70): yellow stars indicate geochemical samples and red stars indicate geochemical samples analyzed also with Raman spectroscopy. Drill cores B1-46 and B1-384 have been acquired ~ 450 m from each other and are indicated by a single star. For more detailed drill core locations, see the supplementary Plate I of Severson (1991)

The lower portion of the Virginia Formation black shale is dominantly composed of carbonaceous argillites, whereas carbonaceous siltstones and graywackes are more abundant in the upper portion (Lucente and Morey 1983). Carbonate concretions are present throughout the Virginia Formation (Lucente and Morey 1983). Quartz, chlorite, muscovite, and lesser albite are the main minerals in the black shale and typical minor constituents include K-feldspar, calcite, ilmenite, anatase, titanite, pyrite, chalcopyrite, sphalerite, and CM (Lucente and Morey 1983; Andrews and Ripley 1989). Carbonaceous material is present in variable amounts and the typical Corg content of the Virginia Formation is 0.3–2.6 wt% (Ripley and Alawi 1988; Ripley et al. 2001a). The CM is biogenic and experienced thermal maturation during diagenesis and early regional metamorphism as evidenced by the traces of preserved amino acids in the Virginia Formation (Niehaus and Swain 1972) and pyrobitumen in the Rove, Thomson, and Biwabik Iron Formations (Swain et al. 1958; Mancuso et al. 1989). The CM is characterized by δ13C values from − 33.8‰ to -27.1‰ (Andrews and Ripley 1989; Ripley et al. 2001b). Pyrite is the most abundant sulfide in the Virginia Formation, whereas chalcopyrite and sphalerite are present in lesser amounts (Lucente and Morey 1983; Virtanen et al. 2021). The bulk sulfides yield δ34S in the range of 4.5 − 19.1‰ (Ripley 1981; Thériault and Barnes 1998; Ripley et al. 2001b; Queffurus and Barnes 2014).

At ∼1096 Ma, the Duluth Complex formed as repetitious emplacement of sills causing contact metamorphism and partial melting of the Virginia Formation black shale (Park et al. 1999; Ripley et al. 2007; Sawyer 2014; Swanson-Hysell et al. 2021). The magma emplacement temperature has been estimated to be roughly 1200 ± 50 °C (Virtanen et al. 2022). The contact metamorphic aureole of the Virginia Formation can be divided into the Devolatilization zone and Partial melting zone (Fig. 2). The Devolatilization zone is characterized by increase in grain size and replacement of chlorite and muscovite by cordierite and biotite due to subsolidus dehydration reactions at ∼500–600 °C (Labotka et al. 1984; Ripley and Alawi 1988; Sawyer 2014). Regarding the minor phases, CM and pyrite have been replaced by graphitic material and pyrrhotite (Ripley et al. 2001a; Sawyer 2014). In the Partial melting zone, peritectic melting of biotite at ≥ 700 °C is evidenced by in situ neosomes, i.e., recrystallized partial melt pockets (Severson and Hauck 2008; Sawyer 2014). Orthopyroxene became stable in the contact aureole at 740–790 °C and complete melting of biotite, quartz, and K-feldspar melt followed at 800–850 °C (Sawyer 2014). The peak contact metamorphic temperature at the contact exceeded 900 °C as determined from the contact-metamorphosed granites adjacent to the Virginia Formation (Benkó et al. 2018). The partial melting caused shearing, which disrupted the original bedding and the post-depositional strike-slip faults (Severson and Hauck 2008; Sawyer 2014). Sulfides in the Devolatilization and Partial melting zones of the contact aureole cover a wider range of δ34S between 0.1‰ and 25.8‰ compared to the regionally metamorphosed black shale (Mainwaring and Naldrett 1977; Ripley 1981; Ripley and Al-Jassar 1987; Park et al. 1999; Queffurus and Barnes 2014). It should be noted, however, that the literature data for sulfur isotopes from the regionally metamorphosed black shale is scarce compared to the contact-metamorphosed black shale.

Fig. 2
figure 2

Schematic presentation of the contact between the Virginia Formation and the Duluth Complex after Severson et al. (1996). The metamorphic assemblages and estimated temperatures are from Lucente and Morey (1983), Labotka et al. (1984), Sawyer (2014), and Benkó et al. (2018). The figure is not in scale. The abbreviations are as follows: Qz = quartz, Chl = chlorite, Ms = muscovite, Ab = albite, Bt = biotite, Crd = cordierite, Kfs = K-feldspar, Pl = plagioclase, Opx = orthopyroxene

The contact metamorphism caused widespread devolatilization in the Virginia Formation (e.g., Rao and Ripley 1983; Park et al. 1999; Ripley et al. 2001a). The devolatilization fluids redistributed H2O, Corg, and sulfur within the contact aureole and transported these components to the Duluth Complex intrusions contributing to the formation of the Cu-Ni(-PGE) deposits (Ripley 1981; Park et al. 1999; Ripley et al. 2007; Virtanen et al. 2021). Sulfur mobilization in the contact aureole was efficient due to the presence of a reduced fluid that formed due to the breakdown of CM and silicate dehydration reactions (Andrews and Ripley 1989; Ripley et al. 2007; Virtanen et al. 2021). Advective fluid transport of elements has been envisioned where suitable faults and layer parallel structures were available, whereas slower diffusive transport likely dominated elsewhere (Park et al. 1999). Based on oxygen and hydrogen isotopes as well as altered sulfide compositions and trapped low-temperature fluid inclusions, fluids infiltrated the magmatic rocks of the Duluth Complex still at the post-magmatic stage (Park et al. 1999; Ripley et al. 2007 Gál et al. 2011; Benkó et al. 2015).

The BPU is locally present within the lowermost 60 m of the Virginia Formation as distinctive 1–30 m thick horizons of thinly laminated argillite enriched in CM, sulfur, and metals (Miller et al. 2002; Fig. 2). The BPU has been mostly characterized based on sporadic occurrences in the contact aureole and as xenoliths (Fig. 2) within the Duluth Complex (Severson et al. 1996; Hauck et al. 1997b; Severson and Hauck 2008; Queffurus and Barnes 2014; Samalens et al. 2017a), although rare occurrences are known outside the contact aureole (Williams et al. 2010; Samalens et al. 2017a). The BPU xenoliths were an important source of sulfur for some of the massive Cu-Ni(-PGE) occurrences of the Duluth Complex (Queffurus and Barnes 2014; Samalens et al. 2017a). The BPU contains approximately 5–20 vol% of sulfides, which are present as 1–3 mm thick laminae and lenses that either conform with the original bedding or crosscut the metamorphic recrystallized silicate fabric (Severson et al. 1996). Pyrrhotite is typically the dominant sulfide (∼90 vol%), but chalcopyrite, cubanite, pentlandite, sphalerite, galena, and marcasite are also present (Hauck et al. 1997b; Williams et al. 2010; Queffurus and Barnes 2014). The sulfides in the BPU show highly variable and generally higher δ34S (5.6–29.1‰) than sulfides within the normal black shale of the Virginia Formation (Ripley and Al-Jassar 1987; Severson 1994; Zanko et al. 1994; Arcuri et al. 1998; Thériault and Barnes 1998; Ripley et al. 2001b; Queffurus and Barnes 2014; Samalens 2017). The CM in the BPU is isotopically similar with the rest of the Virginia Formation having δ13C in the range from − 33.6‰ to -29.2‰ (Ripley et al. 2001b). The silicate mineral assemblages within the regionally metamorphosed and contact-metamorphosed BPU are similar to the normal black shale of the Virginia Formation.

Materials and methods

Sample materials

We refer to the portion of the Virginia Formation black shale that is free of sulfide laminae as normal black shale (NBS), to distinguish it from the BPU. We use suffix -CA when referring to the NBS or BPU that experienced contact metamorphism. The thin and thick sections for the Raman spectroscopic and petrographic characterizations come from the NBS (n = 3), BPU (n = 1), and from the BPU-CA within the Devolatilization zone (n = 2) (Fig. 3). The NBS samples come from the drill core B2 (Fig. 1) from the depths of 390 m (DC3), 443 m (DC8), and 478 m (VF-BS1). The sample LTV 7555-01 comes from the drill core LTV 7555 (Fig. 1) from the depth of 38 m. The sample DC70 represents the BPU-CA and comes from an outcrop at the contact with the Dunka Pit Cu-Ni deposit (Fig. 1). The sample B1-338-01 is from depth of 79 m from the drill core B1-338, which samples the BPU-CA in the footwall contact of the Babbitt Cu-Ni deposit (Fig. 1).

Fig. 3
figure 3

Photographs of representative thin section samples used in the Raman spectroscopy measurements. The three vertical lines in the LTV 7555-01 are laser ablation tracks made in an earlier study (Samalens et al. 2017b). The samples DC70 and B1-338-01 are from the Devolatilization zone. The Bedded Pyrrhotite Unit samples contain pyrrhotite laminae or lenses, which are indicated with white arrows. The length of the scale bar for each sample is 1 cm

We included the samples listed above for volatile chemistry analyses with two exceptions: sample DC70 for which suitable material was not available and sample VF-BS1 for which the necessary data already existed (Virtanen et al. 2021). We included one additional NBS sample (DC5) from the drill core B2 (Fig. 1) from the depth of 465 m. Contact-metamorphosed NBS-CA is represented by samples A4-18-25 (depth of 197 m) and A4-18-26 (199 m) that are from the drill core A4-18 (Fig. 1). For the BPU-CA, we included seven additional samples: sample B1-46-11 is from the depth of 93 m from the drill core B1-46 (Fig. 1); samples A4-18-03 (3 m), A4-18-04 (17 m), A4-18-09 (65 m), and A4-18-11 (87 m) are from the drill core A4-18 (Fig. 1); samples B1-384-14 (336 m) and B1-384-26 (325 m) are BPU xenoliths from the drill core B1-384 (Fig. 1).

The samples with DC prefix have been studied previously for sulfur isotopes (δ34S) and trace elements by Thériault and Barnes (1998), Queffurus and Barnes (2014), and Samalens et al. (2017a, b). The sample with LTV prefix has been studied for trace elements by Queffurus and Barnes (2014) and Samalens et al. (2017a, b). Sulfur isotopic (δ34S) data has been acquired previously from the samples A4-18-04, A4-18-09, B1-384-14, B1-384-26 by Thériault and Barnes (1998), Queffurus and Barnes (2014), and Samalens (2017).

Petrography

We conducted the petrographical observations and qualitative phase identification at the Helsinki Geoscience Laboratories (Hellabs) of the University of Helsinki. We used a Leica DM2700 P microscope equipped with a Leica DFC450 C camera for transmitted and reflected light microscopy. For qualitative phase identification, we carbon-coated the thin and thick section samples and analyzed them with a JEOL JSM-5900LV Scanning Electron Microscope (SEM) equipped with a SAMx energy-dispersive silicon drift detector (EDS). Control software of the SEM and the EDS detector are Digital Image Scanning System 5 and IDFix, respectively. For the EDS measurements, we set the working distance to < 15 mm and the signal measurement time to 10 s. We used acceleration voltage of 15 kV in the measurements, except for sulfide identification for which we used an acceleration voltage of 20 kV.

We supplemented the qualitative identification of aluminum silicate phases subject to polymorphism with Raman spectroscopy. We conducted the analyses at the Department of Chemistry, University of Helsinki, using a confocal NT-MDT Ntegra Raman spectrometer equipped with a 532 nm frequency doubled Nd: YAG laser and an Andor Newton CCD detector. For the analyses, we used a 100× magnification objective and set the grating to 1800 gr/mm.

Raman spectroscopy of carbonaceous material

We used Raman spectroscopy to characterize the CM in the NBS, BPU, and BPU-CA. The physical structure of CM is sensitive to temperature changes in diagenetic and metamorphic conditions (e.g., Beyssac et al. 2002; Lahfid et al. 2010; Henry et al. 2019). With increasing temperature, the disordered precursor CM is subject to gradual structural ordering in a process called graphitization until it finally transforms into perfectly crystallized graphitic material at around 640 ± 50 °C (Beyssac et al. 2002). The degree of structural ordering of the CM can be quantified from its Raman spectral features. Fully ordered graphite structure is identified from a single G-band at ∼1580 cm− 1 in the Raman spectrum, while structural defects give rise to various D-bands (e.g., Beyssac et al. 2002; Lahfid et al. 2010; Beyssac and Lazzeri 2012). Many studies have shown that the relative sizes of the G- and D-bands are a useful measure of the structural ordering of CM (e.g., Beyssac et al. 2002; Lahfid et al. 2010; Beyssac and Lazzeri 2012; Henry et al. 2019). When the precursor CM has been subject only to graphitization sensu stricto, the structural ordering can be used as a proxy for the peak (diagenetic or metamorphic) temperature conditions (e.g., Beyssac et al. 2002; Lahfid et al. 2010; Henry et al. 2019). As an example, the thermometry may not be valid for graphitized migratory bitumen, which forms from the residual CM during catagenic stage of graphitization. Several types of CM with different structural ordering can coexist in the same rock and concurrently experience subsequent temperature-dependent structural changes but only the structural ordering of the residual CM is representative for the full range of graphitization.

We measured the structural ordering of the CM with a Dilor LabRam II confocal micro-Raman spectroscope at the Department of Earth Sciences, ETH Zürich. The Raman spectroscope is equipped with a 532.11 nm diode-pumped solid-state laser and an Olympus BX40 microscope. All the measurements were conducted using the 100× magnification objective, which provides spatial accuracy of ∼1 μm. The laser power on the sample was limited to < 1 mW to avoid overheating and a notch filter was applied to block the Rayleigh scattering from the sample. The confocal hole size was 1120 μm with samples DC70 and B1-338-01, and 1128 μm with DC3, DC8, VF-BS1, and LTV 7555-01. A 1800 gr/mm grating was used to disperse the signal to a Peltier-cooled CCD detector with 1024 × 256 pixels.

Before each analysis session, we calibrated the detector with a silicon standard. We measured samples from approximately similar metamorphic conditions within one analytical session to avoid instrument-related uncertainties arising, for example, from the different confocal hole sizes. The measurement time was 30 s and the number of spectral accumulations was five. We measured the CM in confocal mode through transparent minerals to avoid contributions from possible structural degradation caused during polishing of the samples (see Beyssac and Lazzeri 2012). To account for possible sample heterogeneity and the effects related to crystallographic orientation, we took 13–15 measurements from each sample. From the BPU sample (LTV 7555-01), we measured CM spatially related to the groundmass and to the fracture zones (see Results).

The PeakFit software (version 4.12) was used to analyze the measured Raman spectra. Each spectrum was smoothed to remove noise from the measurement. A linear baseline was visually fitted to the spectra to remove background signal from the measured bands. Measurements with intense fluorescence signal were excluded due to the ambiguities in determining a proper baseline for the background signal. We deconvolved the spectra using the Voigt functions (see e.g., Beyssac et al. 2002; Lahfid et al. 2010). Five bands (G, D1, D2, D3, and D4) were required to deconvolve the Raman spectra of the NBS samples DC3, DC8, and VF-BS1, whereas three to four bands (G, D1, D2 ± D3) were required with the CM in the BPU sample LTV 7555-01. Only the G-band was required for spectra from the BPU-CA samples DC70 and B1-338-01.

We performed thermometry calculations based on the quantified Raman spectra of the CM measured from the NBS samples. We applied the temperature calibration of Lahfid et al. (2010), for which we calculated the RA1 parameter, which is a ratio of the bands in the Raman spectra: (D1 + D4)/(G + D1 + D2 + D3 + D4). The peak metamorphic temperature was calculated following the linear calibration of Lahfid et al. (2010): T= (RA1-0.3758)/0.0008. This thermometer is calibrated for temperatures in the range of 200–320 °C and has a typical uncertainty of 50 °C. As CM in the BPU formed due to processes other than simple graphitization (see discussion) and because CM in the BPU-CA is fully graphitized, these samples were excluded from the thermometry calculations.

Geochemistry

The bulk rock geochemistry for powdered (< 80 μm) samples was conducted at the Service d’Analyse des Roches et des Minéraux, Centre de Recherches Pétrographiques et Géochemiques, Nancy. Sulfur and Corg analyses were performed with a Horiba EMIA-320V2 Carbon-Sulfur analyzer. For sulfur measurements, untreated sample powder was induction-combusted in an oxygen-atmosphere and the amount of sulfur in the combustion products was measured by infrared absorption. The amount of Corg was measured with the same procedure, but the sample powder was pretreated with HCl to remove possible carbonate-bound C. For the measured sulfur contents of 0.1–0.5 wt%, the coefficient of variation (CV) is < 20%, for sulfur contents of 0.5–1.0 wt%, the CV is < 10%, and for sulfur contents of > 1.0 wt%, the CV is < 5%. For the measured Corg contents of < 0.5 wt%, the CV is < 15% and for Corg contents of > 0.5 wt%, the CV is < 5%.

The amount of mineral-bound H2O (H2O+) in the bulk rock was determined by Karl Fischer (KF) titration. The sample powder was induction-heated at 1000 °C and N2 was used to transport the pyrolysis products to the titration vessel with KF reagent. The amount of combusted H2O+ was determined with coulometric titration. The CV for the H2O+ measurements is < 5%.

Results

Phase identification and structural observations

We describe the details of the identified phases and structures that are relevant to understanding the distribution Corg, H2O, and sulfur between the different black shale types of the Virginia Formation. The studied NBS samples represent typical Virginia Formation black shale consisting of modally subequal quartz, chlorite, and muscovite with lesser albite (Fig. 4). The grain size is mainly 1–40 μm. The CM is present as evenly distributed micron-scale dispersed material in the mineral interstices (Fig. 4a). Euhedral pyrite is the most abundant sulfide and present both as dispersed individual grains and as lenticular aggregates generally composed of hundreds of grains (Fig. 4). The pyrite aggregates are conformable with the bedding and associated with CM accumulations (Fig. 4a). Chalcopyrite is present as dispersed grains in the groundmass and together with sphalerite they are present as inclusions between chlorite cleavage planes.

Fig. 4
figure 4

Reflected light photomicrograph (a) and backscattered electron image (b) of the normal black shale sample VF-BS1 from outside the contact aureole. The bedding is in the horizontal direction in a and in vertical orientation in b. Pyrite is present as individual grains and as a lenticular grain aggregate, which is associated with carbonaceous material. The abbreviations are as follows: Qz = quartz, Chl = chlorite, Ms = muscovite, Ab = albite, Py = pyrite, and CM = carbonaceous material

The groundmass texture, grain size, and mineralogy of the BPU are mostly similar to the NBS samples. The CM in the groundmass is dispersed in the interstitial space and occurs as micron-scale flakes as in the NBS. Sulfides are present as sparse dispersed individual grains of pyrrhotite and minor chalcopyrite, whereas pyrite is absent. The sample also contains fracture zones parallel to the original bedding (Fig. 5a–d). The fracture minerals show misalignment relative to the original bedding and they are distinctively coarser than the groundmass (1–40 μm) with quartz grains reaching > 200 μm and pyrrhotite several hundreds of µm in length (Fig. 5a–d). In the fracture zone, quartz is morphologically diverse. Subhedral equant grains are common and locally enclosed in the largest sulfides. Anhedral quartz is characterized by a variety of elongated shapes that have rod-like (Fig. 5a–c), angular (Fig. 5d), and wavy (Fig. 5d) cross-sections indicating replacement of precursor clasts. These replacement grains or grains aggregates are locally partly composed of sulfides and often partly or fully coated with CM (Fig. 5a–d). Quartz is also present as a few microns thick layers on sulfides and CM (Fig. 5a–c). In the fracture zones, the CM is coarser than in the groundmass exceeding 100 μm in size and it commonly hosts TiO2 phase and titanite inclusions (Fig. 5a–d). Sulfides associated with the fracture zones are mostly pyrrhotite with lesser chalcopyrite, galena, and marcasite. Anhedral quartz flanks the largest sulfides and, together with calcite, is present as inclusions in the sulfides (Fig. 5a–b). Chalcopyrite is typically present at the rims and fractures of the largest pyrrhotite grains (Fig. 5a).

Fig. 5
figure 5

a Reflected-light photomicrograph and b backscattered electron (BSE) image of the same view of the Bedded Pyrrhotite Unit sample LTV 7555-01 from outside the contact aureole. The left side of the images is dominated by a fracture zone enriched in sulfides, carbonaceous material (CM), TiO2 phase, titanite, and quartz (Qz) that are larger than the groundmass on the right side. Locally, CM occupies the center of the authigenic grains. c and d show BSE images from the same sample showing the variable morphology of CM, pyrrhotite, and Qz in the fracture zones. Sulfides and Qz are authigenic and generally coated with CM. The bedding is in the vertical orientation in ad. The abbreviations are as follows: Qz = quartz, Chl = chlorite, Ms = muscovite, Ab = albite, Po = pyrrhotite, Ccp = chalcopyrite, Cal = calcite and Ttn = titanite

The main silicate minerals in the groundmass of the BPU-CA sample DC70 from the Devolatilization zone are cordierite, quartz, plagioclase, K-feldspar, and biotite (Fig. 6a). Graphitic material, sillimanite, rutile, and apatite are the most common minor phases. Quartz, plagioclase, and K-feldspar grains are generally 50–100 μm in diameter, whereas cordierite is present as mm-sized poikilitic grains enclosing the other main silicates and graphite (Fig. 6a-b). Cordierite is generally in contact with K-feldspar, which locally forms nearly monomineralic rims around the cordierite grains (Fig. 6a). Graphitic material is ubiquitous and present as flakes that reach ∼300 μm in the longest dimension (Fig. 6). The flakes are often intergrown with cordierite (Fig. 6a), biotite (Fig. 6c-d), and pyrrhotite (Fig. 7a). Anhedral pyrrhotite with grain size of 20–150 μm is the dominant sulfide in the groundmass (Fig. 6a). It is often rimmed and intergrown with graphitic material flakes (Fig. 7a) and less commonly with rounded rutile grains.

Fig. 6
figure 6

Backscattered electron images (a-d) of the Bedded Pyrrhotite Unit sample DC70 from inside the contact aureole. The original bedding is in the vertical orientation. a Groundmass of the sample consists of quartz, plagioclase, K-feldspar, cordierite, and phlogopitic biotite. Cordierite is poikilitic and preferentially rimmed with K-feldspar (lightest grey phase). b Muscovite is exclusively present between cordierite and pyrrhotite laminae. c Large phlogopitic biotite often hosts anhedral cordierite inclusions and is typically present in contact with pyrrhotite laminae. Minute quartz inclusions are present in biotite and in the anhedral cordierite grain indicated with an arrow. d Thin rim of quartz is present between pyrrhotite lamina and phlogopitic biotite. The biotite grain contains a sillimanite inclusion. The abbreviations are as follows: Qz = quartz, Crd = cordierite, Kfs = K-feldspar, Pl = plagioclase, Gr = graphitic material, Po = pyrrhotite, Ms = muscovite, and Phl = phlogopitic biotite

Fig. 7
figure 7

Reflected-light photomicrographs of the contact-metamorphosed Bedded Pyrrhotite Unit sample DC70 from inside the contact aureole showing groundmass pyrrhotite intergrown with graphitic material (a) and chalcopyrite on the flanks and fractures of the pyrrhotite lamina (b). The data related to the laser pits and tracks were published in Samalens et al. (2017b). The abbreviations are as follows: Gr = graphitic material, Po = pyrrhotite, Ccp = chalcopyrite

Characteristic mm-scale pyrrhotite laminae traverse the BPU-CA samples parallel to the bedding (Fig. 3). Chalcopyrite is present at the rims and fractures of the pyrrhotite laminae (Fig. 7b). The sulfide laminae are separated from the typical groundmass minerals by a zone of phengitic muscovite (Mg-bearing), phlogopitic biotite, quartz, K-feldspar, and sillimanite, whereas cordierite is not in direct contact with the sulfide laminae (Fig. 6b-d). Anhedral cordierite together with typically rounded quartz, plagioclase, and K-feldspar are found as inclusions in biotite at the vicinity of sulfide laminae (Fig. 6c-d). Muscovite and biotite cleavage planes do not conform with the bedding (Fig. 6b-d) and elongated sulfides with thicknesses from < 1 μm to < 10 μm are commonly observed between the cleavage planes (Fig. 6c). Less than 10 μm thick rim of quartz is locally present in between the pyrrhotite laminae and micas (Fig. 6b and d). A zone depleted in cordierite and enriched in granoblastic quartz and feldspars, is present in between the hydrous silicates and the typical groundmass.

Raman spectroscopy of carbonaceous material

We present the representative first-order Raman spectra of the CM for each sample in Fig. 8 and the calculated structural parameters for the NBS samples in Table 1. The Raman spectra of the biogenic CM in the NBS samples are typical of disordered CM, with well-developed disordered bands (D1-, D2-, D3-, and D4-bands) (Fig. 8). This means the CM is aromatic, but also turbostratic (i.e., the CM did not reach the stacked 3D ordering of graphite). The structural parameter RA1 for these samples is homogeneous and clustered in a narrow interval of 0.62 to 0.65 (Table 1). The CM from the BPU sample (LTV-7555-01) is structurally distinctive from the CM in the NBS as it displays a narrower D1 band of higher amplitude, while other disordered bands show lesser contribution to the spectra (Fig. 8). The structural parameter RA1 for these samples is 0.68 and no difference is observed between the CM in the groundmass and in the fracture zones with both subsets showing the same Raman spectral features and yielding the same RA1. The Raman spectra of the CM from the BPU-CA samples (DC70 and B1-338-01) differs drastically from the other samples. It generally shows a narrow G-band and a lack of disordered bands (Fig. 8). This is the characteristic spectral signature of a perfectly graphitic material (see e.g., Beyssac et al. 2002).

Table 1 Raman structural parameters for the carbonaceous material in the Virginia Formation black shale samples. The samples with perfectly graphitic material (DC70 and B1-338-01) are not shown as the Raman parameter is consistently zero for them due to the lack of disordered bands
Fig. 8
figure 8

Representative first-order Raman spectra of carbonaceous material measured from the Virginia Formation black shale. The black vertical stippled lines indicate the ideal locations of the graphite-band (G) and the various defect bands (D1, D2, D3, and D4). The temperatures for the normal black shale samples outside the contact aureole were determined with the geothermometer based on the RA1 parameter of Lahfid et al. (2010). The average RA1 (n = 14–15) is shown for the disordered carbonaceous materials. See text and Table 1 for more information of the geothermometry

Carbon, sulfur and H2O contents

The NBS has variable volatile content on a decameter-scale within a single drill core (Table 2). In the NBS samples, the Corg content is 0.51–1.97 wt%, sulfur content is 0.32–1.08 wt%, and H2O content is 3.14–4.05 wt% (Table 2; Fig. 9). Compared to the predominant NBS, the BPU sample has much higher Corg content of 4.25 wt% and lower H2O content of 2.59 wt% (Table 2; Fig. 9a). The sulfur content of the BPU is 1.07 wt%, which is in the range of sulfur contents measured from the NBS (Table 2; Fig. 9b). The two contact-metamorphosed NBS-CA samples from the Devolatilization zone have Corg contents of 0.75–2.08 wt% and sulfur contents of 0.64–0.86 wt%, which are similar to the contents measured from rest of the NBS (Table 2; Fig. 9). The spatial variability of Corg content is large considering that the NBS-CA sample locations are separated by less than two meters in the A4-18 drillcore (Table 2). The NBS-CA samples are slightly less hydrous than the NBS samples with H2O contents of 1.81–2.02 wt% (Table 2; Fig. 9). Compared to the NBS-CA volatile contents published in literature, our samples are relatively enriched in volatiles (Fig. 9). Compared to the NBS, the volatile contents of the contact-metamorphosed BPU-CA samples are very variable with Corg of 0.44–4.42 wt%, sulfur contents of 2.29–7.37 wt%, and H2O contents of 0.93–3.26 wt% (Table 2; Fig. 9). Typically, the Corg and sulfur contents of the BPU-CA are much higher than in the NBS, whereas the H2O contents are overlapping or lower (Table 2; Fig. 9). Compared to the NBS-CA, which is the prevalent contact-metamorphosed type of black shale, the BPU-CA is enriched in volatiles (Fig. 9).

Table 2 Organic carbon (Corg), sulfur (S) and H2O data from the Virginia Formation black shale samples
Fig. 9
figure 9

H2O vs. organic carbon (Corg) (a) and sulfur (S) (b) for Virginia Formation black shale samples from outside and inside of the contact aureole. Xenolith samples from the Bedded Pyrrhotite Unit are marked with arrows. The reference data are from Ripley and Alawi (1988) and Virtanen et al. (2021)

Discussion

Formation of the Bedded Pyrrhotite Unit

The presence of chlorite and muscovite suggests that the Virginia Formation experienced subgreenschist to greenschist facies metamorphic conditions in the 1786 Ma regional metamorphic event (Lucente and Morey 1983; Vallini et al. 2007) with the peak temperature reaching roughly 300–400 °C (Ripley et al. 2001a). Our Raman measurements indicate peak temperatures in the range of 300–340 ± 50 °C (Fig. 8; Table 1) supporting the previously interpreted peak regional metamorphic temperatures. This pronounced conformity between the two methods of temperature estimation indicates that the disordered CM in the NBS represents original biogenic CM that was subjected to graphitization sensu stricto. Importantly, this suggests that the material had already reached the main oil and gas window of maturation (e.g., Ma et al. 2016) prior to the formation of the Duluth Complex at ∼1096 Ma (Swanson-Hysell et al. 2021).

Our BPU sample contains a slightly more structurally ordered CM compared to the NBS (Fig. 8). This could be a sign of recrystallization of the BPU-hosted CM at a slightly higher temperature compared to the CM in the NBS as the CM geothermometers would yield higher temperatures with higher structural ordering (e.g., Beyssac et al. 2002; Lahfid et al. 2010). The volatile chemistry, however, shows that the BPU has much higher Corg content than the NBS (Fig. 9), which cannot be explained by closed-system graphitization of the precursor CM of the NBS. It is clear that CM in the BPU is concentrated in the fracture zones (Fig. 5), unlike in its NBS counterpart where the CM is homogeneously distributed (Fig. 4a). The observed features in the BPU indicate that the CM did not form by simple graphitization of the original biogenic CM during regional metamorphism. For this reason, we cannot reliably apply the thermometer of Lahfid et al. (2010) for our BPU sample.

Fig. 10
figure 10

A conceptual model illustrating the sequential formation of the carbonaceous material, quartz, and sulfide-bearing fracture zones in the Bedded Pyrrhotite Unit. The figure is not in scale

The CM in the BPU exhibits morphological similarity with migratory bitumen (occurrence in cavities and as coatings on mineral grains, see Fig. 5a–d) found commonly in oil deposits (e.g., Sanei 2020; Rasmussen et al. 2021). We consider it possible that the CM in the BPU represents oil that accumulated into the microfracture zones and was subsequently converted to pyrobitumen during the 1786 Ma regional metamorphic event. Near the Virginia Formation, pyrobitumen has been identified in the Rove and Thomson Formation black shales as well as in the Gunflint and Biwabik Iron Formation (Swain et al. 1958; Mancuso et al. 1989; Rasmussen et al. 2021), which implies that oil formed prior to the main stage of the regional metamorphic event. As in our BPU sample, sulfides and quartz have been observed to be often closely associated with CM in oil deposits (e.g., Payzant et al. 1986; Parnell et al. 1994). Similar authigenic quartz grains coated and filled with CM that we find in the BPU (Fig. 5), have been described from the Gunflint Formation (Rasmussen et al. 2021), which is the stratigraphic counterpart of the Biwabik Iron Formation on the northern flank of the Duluth Complex (Fig. 1). Rasmussen et al. (2021) interpreted that the CM in the Gunflint Formation represents infiltrated oil that both wetted and became trapped within partly silicified microfossils. Oil also has the ability to mobilize the usually immobile Ti in diagenetic conditions as evidenced by TiO2 phase and titanite inclusions commonly found in pyrobitumen residuals (e.g., Galvez et al. 2012; Fuchs et al. 2015; Liu et al. 2019). In the fracture zones of our BPU sample, we observe sulfide laminae rimmed with quartz and Ti phase inclusions in the CM (Fig. 5b), which support the suggested formation of the CM via oil infiltration prior to the peak stage of regional metamorphism. The oil infiltration hypothesis is strongly supported by the textural features (fracture fillings, mineral grain coatings, see Fig. 5a–d) and the notably different structural ordering between the NBS and BPU-hosted CM (Fig. 8).

Based on the observation that the authigenic quartz and the sulfide minerals are present in the same horizons of the BPU as those infiltrated by oil, the processes could be interdependent. Textural evidence indicates that the authigenic quartz and sulfide grains in the fracture zone precipitated towards pore space that was originally filled with clasts of precursor sedimentary material (Fig. 5a‒d), possibly similar to the lenticular carbonate clasts described from the Virginia Formation (Lucente and Morey 1983). Authigenic quartz and sulfides are common in sedimentary rocks, where they replace calcareous fossils or evaporitic clasts (e.g., anhydrite), which readily dissolve into the pore water during diagenesis or hydrothermal processes (e.g., Birnbaum and Wireman 1985; Bustillo et al. 1999; Machel 2001; Pujatti et al. 2021; Rasmussen et al. 2021). In the case of the BPU, such biogenic or abiogenic clasts could have originated from the same source as in the stratigraphically underlying cherty-calcareous Biwabik Iron Formation (White 1954; Lucente and Morey 1983; Morey 1992). Indeed, the BPU horizons have been currently identified only within 60 m of the Biwabik Iron Formation (Miller et al. 2002). We suggest that the characteristic features of the BPU, namely fractures zones enriched in CM, quartz, and sulfides, formed as a consequence of porosity, which controlled the water and oil infiltration as well as the dissolution of the original clasts (Fig. 10). The process was likely initiated by dissolution of the original clasts and subsequent precipitation of silica in the pore space (Fig. 10). Sulfide precipitation followed silicification as evidenced by quartz that forms the outer rims of most of the authigenic grains aggregates and locally shows euhedral growth to the pore space now occupied by the sulfides (Fig. 10, see also Fig. 5a‒b). As the infiltrated oil replaced the pore water, silica precipitation came to a halt as has been commonly observed in oil deposits (e.g., Worden and Morad 2000; Rasmussen et al. 2021).

The reason why the BPU is enriched in sulfur and why the main sulfide phase is pyrrhotite rather than pyrite as in the rest of the Virginia Formation has not been previously discussed in detail. Understanding sulfide enrichment and mineralogy in the BPU is important because it is the main source of sulfur in some of the massive Cu-Ni(-PGE) sulfide deposits of the Duluth Complex (e.g.Queffurus and Barnes 2014; Samalens et al. 2017a). As discussed above, local high porosity was the main reason for sulfur enrichment to the BPU. Based on previous studies, the sulfides of the BPU show heavy sulfur isotopes with δ34S of 5.6–29.1‰, which indicates formation either via bacterial sulfate reduction (Zanko et al. 1994; Arcuri et al. 1998; Ripley et al. 2001b; Queffurus and Barnes 2014) or thermochemical sulfate reduction (see e.g., Worden and Smalley 1996; Machel 2001). In both processes, pyrite is a typical end product that forms through preceding sulfide species such as hydrotroilite, mackinawite, greigite, and marcasite, whereas pyrrhotite is rare and at best scarce due to kinetic limitations (e.g., Morse et al. 1987; Schoonen and Barnes 1991; Raiswell and Canfield 1998; Rickard and Morse 2005). Therefore, we consider diagenetic origin of pyrrhotite in the BPU unlikely. Experiments have shown that pyrite decomposes into pyrrhotite during hydropyrolysis (in the presence of H2O) and pyrolysis of CM with the decomposition temperature generally in the range of 250–500 °C and 350–700 °C, respectively (e.g., Chen et al. 2000; Gai et al. 2014; Huang et al. 2015; Ma et al. 2016; Zheng et al. 2022). Because the BPU was clearly more porous (filled with oil ± H2O) and contained less mature CM (i.e., prone to thermal breakdown) than the NBS, it was more susceptible to hydropyrolysis and to lower temperature conversion of pyrite to pyrrhotite. We suggest that pyrite decomposition occurred in the BPU, while it did not take place in the NBS at the regional metamorphic peak temperature of 300–340 ± 50 °C (Fig. 4; Table 1). Importantly, pyrrhotite was already the main sulfide in the BPU during the formation of the Duluth Complex. As pyrrhotite has been shown to be stable during devolatilization in the Virginia Formation black shale (e.g. Ripley 1981; Andrews and Ripley 1989; Virtanen et al. 2021), liberation of sulfur from the BPU to the magma required partial melting of the silicates (e.g.Queffurus and Barnes 2014; Samalens et al. 2017a).

Contact metamorphism of the Bedded Pyrrhotite Unit

The BPU-CA is the most sulfur-rich type of black shale identified in the Virginia Formation (Fig. 9b) and similarly to the regionally metamorphosed BPU, it is enriched in Corg compared to the NBS (Fig. 9a). In the studied samples, the pyrrhotite laminae of the BPU-CA are clearly more extensive compared to the BPU (Fig. 3). We first discuss the contact-metamorphic framework of the BPU-CA and then the possible reasons for the observed high sulfur, Corg, and H2O concentrations.

The BPU-CA experienced the typical dehydration reactions described for the Virginia Formation black shale (e.g., Labotka et al. 1984; Ripley and Alawi 1988; Sawyer 2014), where the preceding muscovite and chlorite are replaced by cordierite + K-feldspar + plagioclase + quartz + biotite (Fig. 6a). The presence of perfectly ordered graphitic material in our samples (Figs. 6 and 8) suggests that the temperature reached 640 ± 50 °C (Beyssac et al. 2002), which is compatible with the earlier studies suggesting that prograde biotite is stable at 500–700 °C (Labotka et al. 1984; Sawyer 2014). Whereas the NBS-CA rocks that experienced the prograde dehydration reactions are typically depleted in H2O, sulfur, and Corg compared to the precursor NBS (Rao and Ripley 1983; Ripley and Alawi 1988; Ripley et al. 2001b), this is clearly not the case for the BPU-CA (Fig. 9). Excluding the BPU xenoliths, the BPU-CA has H2O contents that are generally higher than in the NBS-CA and even reach the concentrations of the NBS and BPU (Fig. 9). This indicates that the BPU-CA experienced hydration that post-dates the prograde contact metamorphism. The hydration stage is evident based on the clear hydration of prograde cordierite, which is replaced by muscovite, biotite, and sillimanite (Fig. 6b-d). The replacing minerals can form via well-known hydration reactions (e.g., Tracy and Dietsch 1982; Ogiermann 2002):

$$\begin{aligned} &Cordierite\hspace{0.17em}+\hspace{0.17em}K-feldspar\hspace{0.17em}+\hspace{0.17em}H_{2}O \\& \to Biotite\hspace{0.17em}+\hspace{0.17em}Sillimanite\hspace{0.17em}+\hspace{0.17em}Quartz\end{aligned}$$
(1)
$$\begin{aligned} & K-feldspar\hspace{0.17em}+\hspace{0.17em}H_{2}O \\&\to Muscovite\hspace{0.17em}+\hspace{0.17em}Quartz\hspace{0.17em}+\hspace{0.17em}K_{2}O \left(aq\right)\end{aligned}$$
(2)
$$\begin{aligned} &Cordierite\hspace{0.17em}+\hspace{0.17em}Quartz\hspace{0.17em}+\hspace{0.17em}H_{2}O\hspace{0.17em}+\hspace{0.17em}K_{2}O \left(aq\right)\\ & \to Muscovite\hspace{0.17em}+\hspace{0.17em}Biotite \end{aligned}$$
(3)

It is clear that these hydration reactions must have occurred after the peak contact metamorphic conditions, which strengthens the view that the hydrothermal circulation below the Duluth Complex extended to the post-magmatic stage (Park et al. 1999; Ripley et al. 2007 Gál et al. 2011; Benkó et al. 2015). The BPU xenoliths that were separated from the Virginia Formation by the magmatic rocks were largely unaffected by the retrograde hydration as evidenced by their lower H2O content (Fig. 9).

There are several possibilities for the high sulfur and Corg contents of the footwall and xenolith samples of the BPU-CA. First, the sulfur and Corg contents could reflect diagenetic enrichment in sedimentary sequence. Based on the one BPU sample that did not experience contact metamorphism, it seems plausible that there could have been large enough Corg enrichment in the diagenetic or regional metamorphic stage to explain the high Corg contents in the BPU-CA (Fig. 9). Although the sulfur contents of the BPU-CA are higher (Fig. 9) and the pyrrhotite laminae clearly thicker and more extensive than in our precursor BPU sample (Fig. 3), this apparent difference could simply reflect the unavailability of samples from the seemingly rare precursor BPU. Nevertheless, we consider it highly unlikely that replacement of similar clasts, as envisioned for the studied precursor BPU sample, could have formed such extensive laminae as observed in the BPU-CA. Shales, however, can contain evaporitic interbeds (Testa and Lugli 2000), which arguably could have been subject to similar replacement reactions as the clasts and form the observed extensive sulfide laminae.

Secondly, it is also possible that the high sulfur and Corg contents of the BPU-CA are partly related to the retrograde hydration reactions, which occurred in contact with the sulfide laminae. The contact-metamorphic devolatilization fluids in the Virginia Formation were mixtures of mainly H2O, CH4, CO2, and H2S (e.g., Ripley 1981; Andrews and Ripley 1989; Virtanen et al. 2021). Precipitation of the retrograde hydrous minerals in the BPU fixed H2O from the fluid, which should have simultaneously lowered the solubilities of sulfur and Corg in the fluid and promote sulfide and graphite precipitation (e.g., Rumble et al 1986; Luque et al. 1998; Reed and Palandri 2006; Galvez and Pubellier 2019; Zhang et al. 2021). The sulfides within the cleavage planes of the retrograde biotite suggest at least some chemical reactivity of the sulfides during the retrograde hydration (Fig. 6c). The high sulfur and Corg contents in the BPU xenoliths, which elude the retrograde hydration observed in the BPU-CA, seem to contradict this hypothesis (Fig. 9). It should be considered, however, that the magmatism at the vicinity of the Virginia Formation was periodically active for ∼25 M.y. from 1109 Ma to 1083 Ma (Swanson-Hysell et al. 2021). Hence, the BPU xenoliths within the lowermost intrusion of the Duluth Complex that formed at ∼1096 Ma (Swanson-Hyssell et al. 2021) could have been subjected to earlier cycle(s) of contact metamorphism. Any direct signs of the putative preceding retrograde hydration reactions would have been erased from the BPU xenoliths when they were subjected to the last peak conditions of contact metamorphism. Moreover, many of the BPU xenoliths experienced extensive partial melting and melt extraction (Queffurus and Barnes 2014; Samalens 2017; Samalens et al. 2017a), which further complicates the interpretation of their phase assemblages as well as sulfur and Corg budgets prior to the magmatic assimilation.

The suggested diagenetic and hydrothermal processes are not mutually exclusive and could have both contributed to the observed H2O, sulfur, and Corg contents of the BPU-CA. Although our data does not permit us to constrain the origin of sulfur and Corg in the BPU, we can be reasonably certain that pyrrhotite was the dominant sulfide phase of the BPU before the Duluth Complex magmatism. This is important because unlike the pyrite-bearing NBS, the BPU was not subjected to sulfur-loss during devolatilization due to the stability of pyrrhotite (e.g. Ripley 1981; Andrews and Ripley 1989; Virtanen et al. 2021). The high melting temperature of pyrrhotite (∼1190 °C: Kullerud et al. 1969), required extensive partial melting of the BPU (Queffurus and Barnes 2014; Samalens 2017; Samalens et al. 2017a) or efficient solid-state shattering of xenoliths (see Deegan et al. 2022) to liberate sulfur to the assimilating magma, which subsequently triggered the formation of the associated Cu-Ni(-PGE) deposits of the Duluth Complex. Consequently, the sulfide deposits that formed due to the assimilation of the BPU are mostly hosted in the norites that show signs of more extensive assimilation compared to the troctolites (Thériault and Barnes 1998; Queffurus and Barnes 2014; Samalens et al. 2017a).

Perspective on the assimilation processes from sulfur isotopes

To provide direct evidence of the magmatic assimilation of sulfur from the black shale, we compiled available sulfur isotopic data from literature for the Virginia Formation (n = 125) and Duluth Complex (n = 645), which clearly show the presence of sedimentary sulfur in the sulfide deposits (Mainwaring and Naldrett 1977; Ripley 1981, 1990; Ripley and Al-Jassar 1987; Mogessie and Stumpfl 1992; Severson 1994; Zanko et al. 1994; Hauck et al. 1997a; Thériault and Barnes 1998; Arcuri et al. 1998; Park et al. 1999; Ripley et al. 2001b; Tharalson 2013; Queffurus and Barnes 2014; Samalens 2017; Smith et al. 2022; Fig. 11). The Virginia Formation is highly heterogeneous in terms of δ34S with the values ranging from 0.1‰ to 29.1‰ (Fig. 11). The NBS (including samples DC3, DC5, and DC8 considered in this study) and NBS-CA, have lower median δ34S values (7.2‰) compared to the BPU-CA (15.3‰ with the data including samples DC70, A4-18-04, A4-18-09, B1-384-14, and B1-384-26 considered in this study), which likely reflects the differences in the diagenetic and regional metamorphic processes that affected the sulfides as discussed above.

Fig. 11
figure 11

Sulfur isotopic data (δ34S in ‰ relative to the Canoyn Diablo Troilite) from the Virginia Formation, Duluth Complex, and Mid Ocean Ridge Basalts (MORB). The data for the Virginia Formation and Duluth Complex is from the following sources: Mainwaring and Naldrett (1977); Ripley 1981), 1990); Ripley and Al-Jassar (1987); Mogessie and Stumpfl (1992); Severson (1994); Zanko et al. (1994); Hauck et al. 1997a); Thériault and Barnes (1998); Arcuri et al. (1998); Park et al. (1999), Ripley et al. (2001b), Tharalson (2013), Queffurus and Barnes (2014), Samalens (2017), and Smith et al. (2022). Duluth Complex samples for which description of rock type was lacking, were excluded from the data attributed to the different rock types. Some samples were analyzed in multiple studies in which case duplicate measurements with identical δ34S values within 0.1‰ were included only once in the dataset. The reference data for the MORB mantle is from Labidi et al. (2013). The abbreviations are as follows: NBS = Normal Black Shale, CA = contact aureole/contact-metamorphosed, BPU = Bedded Pyrrhotite Unit, Troc. = troctolite, MS = massive sulfide, n = number of measurements, avg. = average, mdn. = median, Q1 = first quartile, Q3 = third quartile, max. = maximum value, min. = minimum value, outl. = outlier

The median δ34S value for the Duluth Complex in the reviewed literature data is 9.0‰, which is much higher compared to the median value of the mantle (-0.9‰: Labidi et al. 2013) and relatively close to the median value of the Virginia Formation (10.6‰, Fig. 11). Based on a subset of Duluth Complex samples for which the rock type was clearly defined (n = 346), a systematic increase in δ34S values can be observed from troctolitic rocks with median value of 7.5‰, to noritic rocks with median value of 12.4‰, and to massive sulfides (generally hosted in noritic rocks) with median value of 14.4‰ (Fig. 11). The sulfur isotopic data parallels the estimated amounts of assimilation of black shale (excluding selective fluid assimilation). This is generally constrained to be 0–5 wt% for the troctolitic rocks based on Sr isotopic data (Grant and Molling 1981; cf. Grant and Chalokwu 1992) as well as Nd, Pb, and Os isotopic data (Ripley et al. 1999;). For the noritic rocks, assimilation was likely in excess of 20 wt% based on oxygen isotopes (Ripley et al. 2007) and thermodynamic modeling (Virtanen et al. 2022). The NBS and NBS-CA have similar, although slightly lower, median δ34S values with the troctolitic rocks of the Duluth Complex, whereas the noritic rocks and massive sulfides have higher median δ34S values, which are closer to the median δ34S values of the BPU-CA (Fig. 11). These data are compatible with the interpretation that devolatilization fluids caused breakdown of pyrite in the NBS, which selectively mobilized sulfur and transported it to magma to form the main component of the troctolite-hosted sulfide deposits (Ripley 1981; Rao and Ripley 1983; Ripley and Alawi 1988; Thériault and Barnes 1998; Virtanen et al. 2021). As discussed above, pyrrhotite was already present in the BPU prior to the formation of the Duluth Complex. Therefore, sulfur was conserved during devolatilization, and silicate partial melting was required to liberate sulfur to the magma as suggested in previous studies (Queffurus and Barnes 2014; Samalens 2017; Samalens et al. 2017a). The extensive assimilation of the BPU caused the formation of the noritic rocks associated with the heavy δ34S sulfide occurrences. This example underlines that understanding the diagenetic and regional metamorphic processes of sedimentary rocks is important in the interpretation of isotopic data and constraining magma-sediment interaction processes.

Based on the available sulfur isotopic data, it could be argued that the troctolite-hosted sulfide deposits require an additional sulfur source with heavier δ34S values compared to the NBS such as the BPU (Fig. 11). However, it should be noted that the literature δ34S data for the regionally metamorphosed NBS is currently scarce (n = 13) and comes only from a single drill core (B2 in Fig. 1). Although, this drill core contains ca. 450 m section of the Virginia Formation, six of the 13 samples come within an interval of 11 m and show highly heterogeneous δ34S values ranging between 5.1–19.1‰ (Ripley 1981; Thériault and Barnes 1998; Ripley et al. 2001b; Queffurus and Barnes 2014). Due to the currently limited amount of available data from the NBS, the median δ34S value could be biased and should be treated with caution. For now, we consider that the NBS contains portions that have suitably high δ34S values to supply sulfur to the troctolite-hosted sulfide deposits of the Duluth Complex. If the BPU is indeed required to explain the high δ34S values in the troctolites, some yet undefined process is required to transport pyrrhotite-bound sulfur to the magma without significant assimilation of the silicate partial melt.

Conclusions

Our Raman data combined with volatile chemistry and petrological observations revealed the distinct origin of Corg in the different types of black shales within the Virginia Formation. The NBS with relatively low Corg content contains typical graphitized sedimentary CM with aromatic but turbostratic structure, whereas the CM in the Corg -rich BPU is pyrobitumen that indicates infiltration of oil during the diagenetic or early regional metamorphic conditions. Oil in the BPU infiltrated to the same pore space, where dissolution of precursor sedimentary clasts facilitated quartz and sulfur enrichment. During the formation of the Duluth Complex, the NBS contained pyrite as the main sulfide, whereas pyrrhotite dominated in the BPU. Sulfur mobilization in the NBS was possible during the subsolidus devolatilization stage of the black shale, while partial melting was required to release sulfur from the BPU to the magma. Hydrothermal fluid circulation in the contact aureole continued until the retrograde metamorphic stage and redistributed volatiles to suitable structures, such as to the vicinity of the sulfide laminae in the BPU-CA.