Introduction

Porphyry Cu-Mo ± Re ± Au deposits are one of the most important global resources for base (Cu, Mo), precious and rare (Re, Se, Ag, Te, Au, PGE) metals, which are essential for high-tech applications, such as semiconductors, catalysts, high-temperature superalloys, and photovoltaics (Tabelin et al. 2021). However, the scarcity of these deposits and the poorly constrained element enrichment processes could lead to a supply risk in the near future (Grandell et al. 2016; Kesieme et al. 2019). Porphyry-style deposits form in subduction-related settings, in which hydrous melting of the mantle wedge forms oxidized and volatile-rich magmas, which may release metal-bearing aqueous fluids to the uppermost crust (< 5 km) (Richards 2011; Sillitoe 2010; Wilkinson 2013). Interaction of these fluids with the adjacent host rocks causes pervasive alteration halos surrounding the porphyry intrusions (Seedorff et al. 2005; Sillitoe 2010). Primary single-phase high-temperature (> 600 °C) magmatic fluids of intermediate salinity (5–15 wt% NaCl equiv.) can evolve towards a hypersaline liquid (> 26 wt% NaCl equiv.) by phase separation, which is related to the formation of the potassic alteration (Audétat 2019; Bodnar et al. 2014; Heinrich 2005; Kouzmanov and Pokrovski 2012). Later-stage intermediate-temperature fluids (250 – 350 °C) of moderate salinity (5–20 wt% NaCl equiv.) form shallower zones of sericitic alteration that partly overprint the earlier potassic alteration (Audétat 2019; Heinrich 2005; Kouzmanov and Pokrovski 2012). Acidic (pH < 3) and less-saline vapor-rich fluids (< 10 wt% NaCl equiv.) derived from phase separation and which are influenced by SO2 disproportionation during fluid upflow can mix with meteoric waters forming (advanced) argillic alteration zones at temperatures < 300 °C (Bodnar et al. 2014; Hedenquist and Arribas 2022; Seedorff et al. 2005; Sillitoe 2010). This complex fluid evolution is preserved by alteration assemblages and consecutive hydrothermal quartz-sulfide veins that are classified as A-, B-,C-, D, and E-type veins according to their consecutive formation, mineralogy, and textural characteristics by Sillitoe et al. (2010). The porphyry-epithermal transition and fluid evolution is reflected by the trace element and stable isotope composition of hydrothermal minerals (e.g., pyrite, molybdenite, quartz) and by fluid inclusion micro-analysis (Bodnar et al. 2014; Hutchison et al. 2020; Keith et al. 2022; Maydagán et al. 2015; Rottier and Casanova 2021; Shafieibafti et al. 2021).

Molybdenite is the major host for Re in porphyry systems (Barton et al. 2020), which can locally be associated with a rare rheniumdisulfide mineral (rheniite, ReS2), as known from several NE Greece porphyry occurrences (Barton et al. 2020; Voudouris et al. 2009; 2013a) or from sublimates of volcanic gases (Tessalina et al. 2008). Experimental studies have shown that Re is highly mobile and soluble in high-temperature, oxidized, and Cl-rich fluids whereas Re is immobile and insoluble in more reduced fluids (Xiong et al. 2006; Xiong and Wood 1999, 2002). However, Re systematics in natural molybdenite are more complex than predicted from the experimentally-derived prograde solubility of Re in sulfide-rich fluids. For instance, a Re enrichment has been observed in early high-temperature (Spencer et al. 2015; Stein et al. 2001; Ye et al. 2022; Zhang et al. 2021) and late lower-temperature (McFall et al. 2019; Ren et al. 2018; Shafieibafti et al. 2021) porphyry quartz-sulfide veins associated with potassic and sericitic alteration, respectively. Furthermore, Re shows an affinity to high-salinity liquids (Xiong et al. 2006), but also to low-salinity vapors (Tessalina et al. 2008) depending on the physicochemical fluid parameters (e.g., fO2, fS2, S- vs. Cl-rich) and its dominant species in aqueous solutions (ReCl4 vs. HReO4 vs. ReS2 vs. ReO2) (Bernard et al. 1990; Xiong et al. 2006; Xiong and Wood 2002). In addition, Re can also be mobile in supergene post-depositional environments as Re(OH)4 (e.g., Ciobanu et al. 2013; Xiong et al. 2006), and hence the main magmatic-hydrothermal processes and physicochemical fluid parameters that control the enrichment of Re in porphyry systems remain enigmatic.

In this study, we combine sulfide and quartz trace element and stable S isotope micro-analysis with fluid inclusion microthermometry from consecutive hydrothermal quartz-sulfide vein generations from the Maronia Cu-Mo ± Re ± Au porphyry, NE Greece. We present a comprehensive hydrothermal model that links the spatial and temporal fluid evolution to key processes that can control the Re enrichment in molybdenite from porphyry systems.

Geological overview and sample localities

The Rhodope metallogenic province in NE Greece represents one of the most prospective areas for base, precious, and rare metals in the European part of the Tethyan metallogenic belt, where numerous porphyry and epithermal prospects with an exotic mineralogy and element enrichment have been found (Moritz and Baker 2019; Richards 2015; Voudouris et al. 2019). The northward subduction and slab-retreat of the African plate beneath Eurasia since the Cretaceous initiated intrusive and extrusive calc-alkaline to shoshonitic magmatic activity along the ~ 50 km long NE-SW striking age-progressive Maronia Magmatic Corridor in NE Greece (MMC, Fig. 1a; Jolivet and Brun 2010; Moritz and Baker 2019; Perkins et al. 2018). The ~ 6 km2 wide Maronia pluton (29.8 to 29.6 Ma) consists of high-K, calc-alkaline to shoshonitic monzogabbro, monzodiorite to granite that formed by fractional crystallization and minimal crustal assimilation of primary mafic magmas derived from hydrous melting of a phlogopite-bearing metasomatized mantle (Perkins et al. 2018; Schaarschmidt et al. 2021b). The S-poor and Cl-rich magmas were emplaced at ~ 4.6 to 6.1 km crustal depth, followed by rapid uplift, and the formation of a porphyry microgranite. The release of magmatic fluids at 1–2 km depth caused the formation of the Cu-Mo ± Re ± Au porphyry mineralization (Melfos et al. 2002; Schaarschmidt et al. 2021b). The porphyritic microgranite (Fig. 1b) is associated with concentric zones of pervasive and superimposed porphyry alteration (potassic, sericitic) hosting an extensive hydrothermal quartz-sulfide vein network (Fig. 1c; Melfos et al. 2002, 2020). An epithermal overprint of the sericitic alteration is evidenced by extensive silicification and argillic alteration of the host rocks associated with epithermal vein emplacement (Fig. 1c). The mineralized samples of this study represent a lateral profile along the coastline of the Maronia porphyry (Fig. 1b, c) and consist of A-, B-, D-, and E-type quartz-sulfide veins (Fig. 2; Table 1) from the potassic, sericitic, and argillic alteration. The samples provide systematic insights into the spatial and temporal evolution of a Re-rich porphyry system and its epithermal transition.

Fig. 1
figure 1

Geographical location (a) of the NE-SW trending Maronia magmatic corridor (MMC) in the Rhodope massif in NE Greece. (b) Geological map showing the main lithological units with the sampling localities of the southern part of the Maronia pluton and the mineralized Maronia porphyry system. (c) Hydrothermal alteration map of the Maronia porphyry system. Figure 1b and c are modified after Melfos et al. (2020)

Fig. 2
figure 2

Photographs showing the field relations of the vein types from the Maronia porphyry system. (a) Sinuous A-type vein in potassic-sericitic altered microgranite (b) Sinuous A-type vein crosscut by straight B-type veins. (c) Massive B-type quartz veins containing molybdenite and pyrite. (d) Potassic-sericitic altered microgranite associated with thin sinuous A-type vein and crosscutting relationship between A- and B-type veins. (e) Sericitic altered monzonite crosscut by a massive pyrite D-type vein surrounded by a sericite alteration halo. (f) Sericitic overprinted stockwork showing crosscutting relationship of A- and D-type veins. (g) E-type vein crosscutting altered monzonite. (h) Epithermal quartz-carbonate-pyrite-chalcopyrite-sulfosalt vein crosscutting argillic altered microgranite. (i) Silicified epithermal sample showing pyrite (E1) in the quartz-carbonate matrix as well as a late stage arsenopyrite (E2)-sphalerite-galena-sulfosalt epithermal vein crosscutting the matrix. Abbreviations qz = quartz, mol = molybdenite, py = pyrite, fsp = feldspar, ser = sericite, kln = kaolinite, cal = calcite, ccp = chalcopyrite, ss = sulfosalt, gn = galena, apy = arsenopyrite, sp = sphalerite

Table 1 Characteristics of the mineralization and vein types from the Maronia porphyry based on observations in this study and from Melfos et al. (2002, 2020)

Methods

Imaging techniques

Polished thin (n = 22) and thick (n = 18) sections from the various alteration zones and vein generations were petrographically examined to identify the different mineral phases, textures, and paragenetic relations. Back-scattered electron imaging and energy-dispersive X-ray spectrometry (EDS) were performed at the GeoZentrum Nordbayern by scanning electron microscopy (SEM). A Hitachi TM4000 tabletop was used to identify unknown mineral phases and to avoid effects of mineral inclusions with respect to the following microanalytical studies. A TESCAN Vega 2XMU SEM was used for cathodoluminescence (CL) imaging of hydrothermal quartz to distinguish different growth textures and paragenetic stages.

Microanalytical techniques

Quartz micro-analysis

The trace element composition (Li, Na, Mg, Al, Si, K, Ca, Ti, Fe, Ga, Ge, As, Rb, Sr, Y, Sb, Th, and U) of hydrothermal quartz (n = 335) in selected mineralized samples (n = 8) was determined by LA-ICP-MS at the geophysical, environmental, and mineralogical laboratories (HelLabs) of the Department of Geology and Geophysics, University of Helsinki. A Coherent GeoLas MV 193 nm laser was coupled to an Agilent 7900s ICP–MS operated with a plasma power of 1500 W. Helium (1.05 l/min) and Ar (0.81 l/min) were used as carrier gases and Ar acted as plasma (14.9 l/min) and auxiliary gas (0.9 l/min). A single spot ablation pattern with a frequency of 10 Hz and a fluence of 13 J/cm2 was used minimizing the quartz explosions (Audétat et al. 2015; Rusk 2012). The total analysis time for each spot was set to 120 s, including 40 s of gas blank analysis and 30 s for washout. A beam diameter of 60 to 90 μm was applied according to the size of the analyzed grain. SRM NIST 612 glass (Jochum et al. 2005) was used as the primary external standard and the stoichiometric Si content of quartz (46.7 wt% Si) was used for internal standardization, as the total trace element contents in hydrothermal quartz are typically < 0.5 wt% (Rottier and Casanova 2021). Analytical accuracy was monitored by the repeated analysis of the secondary standards SRM NIST 610 (n = 88), BHVO-2G (n = 39), and BCR-2G (n = 39) (Jochum et al. 2005), which was better than 6.0% for most elements. We note that fluid inclusions are common in hydrothermal uartz, which can affect the ablation signal, and therefore the trace element content of the measured quartz (Maydagán et al. 2015; Rottier and Casanova 2021). We followed the method described by Rottier and Casanova (2021) identifying fluid inclusion-derived elements due to a strong correlation with Na and careful evaluation of the time-resolved LA-ICP-MS depth profiles. Only Rb (R2 = 0.49) and Fe (R2 = 0.34) show a weak correlation with intensity spikes of Na, indicating a potential influence by fluid inclusions. Rare mica and Ti-Fe-rich inclusions (e.g., rutile and titanite) were identified by intensity spikes (electronic supplement). Intensity spikes were excluded from the time-resolved LA-ICP-MS patterns during data processing, and thus the presented data exclusively refer to the hydrothermal quartz composition.

Sulfide micro-analysis

The major element composition (Fe, Cu, Zn, Mo, Re, S) of pyrite (n = 428), arsenopyrite (n = 22) and molybdenite (n = 157) was determined in selected samples (n = 12) by wavelength-dispersive X-ray spectrometry (WDS) using a JEOL JXA-8200 Superprobe at the GeoZentrum Nordbayern. The quantitative analyses were carried out using a focused beam at an acceleration voltage of 20 kV and a beam current of 20 nA. Chalcopyrite (Cu, Fe, S), sphalerite (Zn), native Mo (Mo), and native Re (Re) were used as reference materials. Each EPMA spot was marked on a back-scatter electron image for the subsequent LA-ICP-MS measurements of the same spot.

The trace element composition of pyrite (n = 284), arsenopyrite (n = 12), and molybdenite (n = 80) was determined in selected samples (n = 13) by LA-ICP-MS at the GeoZentrum Nordbayern using a Teledyne Analyte Excite 193 nm laser coupled with an Agilent 7500c ICP–MS. The ICP–MS operated with a plasma power of 1190 W. Helium (0.9 l/min) and Ar (0.95 l/min) were used as carrier gases. Argon was used as plasma (14.9 l/min) and auxiliary gas (0.9 l/min). A single spot ablation pattern with a frequency of 20 Hz and a fluence of 4.04 J/cm2 was used. The total analysis time for each spot was set to 55 s, including 20 s of gas blank analysis prior to ablation. A beam diameter of 30 to 50 μm was applied according to the size of the analyzed mineral. The following primary external standards were used for the analysis: PO724 B2 SRM (Au), MASS-1 (Co, Ni, Cu, Zn, Ga, Ge, As, Se, Mo, Ag, Cd, In, Sn, Sb, Te, Tl, Pb, Bi) and FeNiS1 (Re). To monitor the instrument drift, the sulfide standards were analyzed several times during an analytical session. Analytical precision monitored by the repeated analysis of the primary standards (MASS-1 Po724, FeNiS1) was < 6.0% RSD for all elements except for Te (7.2%) and Cd (7.5%). Analytical accuracy of the secondary standard (UQAC-FeS1) was better than 10.0% for all elements, except for Pb (11.1%), As (12.2%), Zn (13.8%), Se (14.4%) and Ag (16.5%). Intensity spikes, as observed in the time-resolved LA-ICP-MS depth profiles (electronic supplement), were excluded prior to data reduction to avoid effects of mineral inclusions on the trace element data. Sulfur was chosen as the internal standard due to a lack of interferences in the analyzed element spectrum, since there is only a 16O18O interference on 34S of about 1% for all sulfides, which is produced from air (Sylvester 2001).

Fluid inclusion microthermometry

Fluid inclusion data (n = 80) of quartz were obtained from seven double polished wafers of ~ 200 μm thickness from the different vein types. Following the CL imaging to distinguish the different quartz generations, fluid inclusion microthermometry was carried out with a LINKAM THM-600/MS 90 heating-freezing stage, at the Department of Mineralogy, Petrology and Economic Geology at the Aristotle University of Thessaloniki, Greece. The stage was calibrated with liquid/solid reference organic substances (chlorobenzene, dodecane, napthalene) and with distilled water with known melting temperatures, yielding an accuracy of ± 2 °C during heating, and ± 0.5 °C during freezing. Although fluid inclusions along growth zones were noticed, they were too small to be analyzed and only inclusions grouped in individual clusters in intergranular quartz grains could be measured. They were identified due to similarities in size and shape and constant phase ratios based on petrographic observations and are interpreted as primary fluid inclusion assemblages (FIAs) along growth zones (Bodnar et al. 2014; Goldstein and Reynolds 1994). Due to the potential effect of post-entrapment processes, microthermometry of secondary fluid inclusions scattered along fractures or passing across grain boundaries in trails was avoided. The spreadsheet HokieFlincs_H2O-NaCl (Steele-MacInnis et al. 2012) was used for the estimation of the phase relations and properties. Bulk salinities of hypersaline fluid inclusions with halite and sylvite daughter crystals were calculated from dissolution temperatures in the system H2O-NaCl-KCl using the program AqSo1e (Bakker 2003).

Sulfur isotopes of pyrite, molybdenite and arsenopyrite

Pyrite (n = 8), molybdenite (n = 7), and arsenopyrite (n = 3) separates were extracted as powders from the different vein types by a diamond microdrill with a diameter of 0.8 mm. Approximately 150–250 µg of sulfide powder and 400–700 µg of V2O5 were placed into tin capsules and analyzed for their δ34S composition by a Flash EA IsoLink elemental analyzer interfaced to a ThermoScientific Delta V Advantage isotope ratio mass spectrometer (EA-IRMS) at the University of Münster. Reproducibility was determined through replicate analyses of the reference materials (IAEA-S-1, Ding et al. 2001) and duplicates of the sulfide powders, which was better than 0.1‰ (1σ).

In situ analyses of 34S/32S ratios of pyrite (n = 121) were performed at the MAGMA Lab of the TU Berlin using an Agilent 8900 ICP-MS/MS coupled to a Teledyne Analyte Excite 193 nm excimer laser. The interface set-up is identical to the one described in Börner et al. (2022), but included an additional cyclonic mixing chamber (electronic supplement). Analyses were performed using a 50 μm spot, a repetition rate of 10 Hz and a fluence of 2.0 J/cm2 either in spot or in line mode at a scan speed of 3–4 μm/s. Mass bias and drift were corrected using the PPP-1 pyrite (δ34S = 5.3 ± 0.2; Gilbert et al. 2014), which was measured repeatedly (~ 3 spots per hour, > 15 analyses per session). The Balmat pyrite (δ34S = 13.8–16.4‰, Gilbert et al. 2014) was used as a secondary standard and gave 16.7 ± 1.3‰, 2 SE, n = 21 (session 1, 11/10/2022), 16.7 ± 1.4‰, 2 SE, n = 18 (session 2, 12/10/2022), and 15.3 ± 1.3‰, 2 SE, n = 15 (session 3, 13/10/2022). Measurement uncertainty (internal precision) for a single analysis was typically 1.5‰ (1 SD) for spot and 1.3‰ (1 SD) for line analyses. The S isotope data is presented relative to the V-CDT reference material and reported in standard notation (δ34S).

Results

Hydrothermal mineralogy

Alteration and vein classification

The sulfide mineralization at Maronia occurs as disseminated grains (pyrite + molybdenite ± pyrrhotite ± chalcopyrite) in the altered microgranite and in the porphyry A-, B-, D-, and epithermal E-type veins according to the classification by Sillitoe (2010) (Figs. 2, 3, 4 and 5; Table 1). Based on the association with distinct alteration zones and their crosscutting relations the hydrothermal quartz-sulfide veins can be separated into an early to late formation as part of four distinct consecutive hydrothermal stages (Fig. 5). The earliest A-type veins (< 3 cm) are sinuous and quartz-dominated with minor pyrite + pyrrhotite ± magnetite ± chalcopyrite (Fig. 2a). The subsequent straight, massive, and thick (1–10 cm) quartz-dominated B-type veins host pyrite + molybdenite ± chalcopyrite and crosscut the A-type veins (Fig. 2b – d). Both A- and B-type veins are related to vein stockworks in the potassic alteration of the porphyritic microgranite. Straight D-type veins in the sericitic alteration are quartz-poor and host abundant pyrite + molybdenite. The D-type veins partly show sericite halos and rarely crosscut or reopen earlier veins (Fig. 2e, f). Two types of epithermal-style veins (E1 and E2) were distinguished in the argillic alteration zone (Fig. 2g, h) and the silicified host rocks (Fig. 2i). The E1-type veins consist of quartz + calcite + pyrite ± sulfosalts and are associated with disseminated euhedral pyrite that occurs in a silicified massive quartz-calcite gangue (Fig. 2h, i). By contrast, the E2-type veins consist of arsenopyrite + galena + sphalerite ± sulfosalts and crosscut the silicified quartz-calcite matrix that is related to stage E1 (Fig. 2i).

Fig. 3
figure 3

Microscopic and SEM-CL textures observed in the hydrothermal mineralization and hydrothermal quartz from the Maronia porphyry system. (a) Thin quartz-pyrite A-type vein crosscut by massive B-type quartz vein. (b) Granular quartz crystals of A-type veins. (c) CL-bright mosaic quartz and “cobweb” texture(Rusk 2012) with darker CL intensity on quartz rims. (d) CL-bright mosaic quartz crosscut by CL-grey quartz from later generations. (e) Massive granular quartz B-type vein associated with pyrite and molybdenite. (f) B-type vein consisting of granular quartz. (g) Cl-grey B-type quartz vein reopened by CL-dark quartz from D-type quartz resulting in a central suture. CL-grey quartz grows perpendicular to the vein wall. (h) Granular CL-grey quartz shows CL-zoning at the termination to open vugs. (i) Quartz-poor pyrite-rich D-type vein showing a distinct sericite halo. (j) Fine-grained sericite halo mantling fine-grained quartz. (kl) Fine-grained CL-dark quartz of D-type veins. (m) Epithermal silicified sample containing E1 pyrite in the euhedral quartz + calcite matrix. The euhedral CL-black quartz matrix is crosscut by an arsenopyrite (E2)-sphalerite-galena-sulfosalt epithermal vein. (n) Close up of quartz-calcite matrix of silicified epithermal sample. (o) Euhedral CL-black quartz crystals growing into open space in epithermal veins. (p) Strong CL and growth zoning of CL-black quartz in the porous matrix of silicified epithermal samples. Abbreviations qz = quartz, py = pyrite, mol = molybdenite, ser = sericite, apy = arsenopyrite, gn = galena, ss = sulfosalt, sp = sphalerite, cal = calcite, fsp = feldspar

Fig. 4
figure 4

Microscopic and BSE textures observed in the hydrothermal mineralization and sulfides from the Maronia porphyry system. (a) Magnetite- and pyrite-bearing A-type vein. (b) Chalcopyrite, pyrrhotite, and magnetite inclusion in pyrite from A-type veins. (c) BSE image of pyrrhotite inclusions in pyrite. (d) Blade-shaped molybdenite associated with cracked subhedral pyrite from B-type veins. (e) BSE image of tennantite-tetrahedrite, bournonite and molybdenite association in the massive quartz matrix of B-type veins. Small native Au inclusion occur associated with subhedral pyrite grains. (f) BSE image of blade-shaped molybdenite associated with rheniite. (g) Massive pyrite overgrowing rosette- to blade-shaped molybdenite. (h) Cataclastic fracture-rich subhedral pyrite. (i) BSE image of molybdenite inclusion in subhedral pyrite. (j) Epithermal vein consisting of euhedral arsenopyrite (E2), subhedral galena, and subhedral pyrite-marcasite intergrowth. E1-type pyrite is pore-rich and euhedral. (k) Subhedral zone pyrite-marcasite in arsenopyrite E2-type vein. (l) Pyrite-marcasite intergrowth associated with galena and arsenopyrite in E2-type veins crosscutting the quartz-carbonate matrix of silicified host rock. Abbreviations py = pyrite, mag = magnetite, po = pyrrhotite, ccp = chalcopyrite, mol = molybdenite, tn-tt = tennantite-tetrahedrite, bnn = bournonite, qz = quartz, cal = calcite, apy = arsenopyrite, gn = galena, rhn = rheniite

Fig. 5
figure 5

Simplified paragenetic sequence of mineral formation with respect to the different alteration zones and vein types in the Maronia porphyry system, based on observations in this study and from Melfos et al. (2002, 2020). The thickness of the bars reflects the abundances of the different minerals, as estimated by optical and electron microscopy. Abbreviations mag = magnetite, cal = calcite, qz = quartz, py = pyrite, ccp = chalcopyrite, mol = molybdenite, po = pyrrhotite, ss = sulfosalts, rhn = rheniite, apy = arsenopyrite, gn = galena, sp = sphalerite. Note: Sulfosalts include tennantite-tetrahedrite, meneghinite, bournonite, famatinite, chalcostibite, zinkenite, and boulangerite

Quartz textures

Hydrothermal quartz in the different vein types can be distinguished by CL-SEM imaging due to distinct CL textures and intensities (Fig. 3). A-type veins commonly show CL-bright anhedral quartz mosaics with a homogenous texture (Fig. 3a – d). The grain boundaries show a darker CL intensity and a “cobweb” texture (Rusk 2012). Quartz growth zoning and CL darker quartz veins crosscutting the quartz matrix were rarely observed in this vein-type (Fig. 3c, d). B-type veins are dominated by CL-grey quartz mosaics that are locally associated with euhedral quartz that grew into open spaces, and which shows alternating CL-grey and CL-dark growth bands (Fig. 3e – h). Quartz in the B-type veins often formed perpendicular to the vein-selvage. In thick and massive B-type veins (> 5 cm), an open-space central suture is occasionally visible due to the reopening of the vein followed by the precipitation of later-stage CL-dark quartz (Fig. 3g). D-type veins are dominated by subhedral quartz with dark CL intensity (Fig. 3i – l). E1- and E2-type veins host CL-black quartz, which is characterized by large euhedral crystal growth and strongly oscillating CL intensity (Fig. 3m – p). In summary, CL-intensity gets darker from the early- to late-stage veins alongside an increasing abundance of euhedral quartz with growth zoning (Fig. 3).

Sulfide petrography

Pyrite and molybdenite are the most common sulfide minerals (Table 1; Fig. 5) that develop various textures in the different vein types. A-type veins commonly host pyrite with inclusions of chalcopyrite and pyrrhotite, which is rarely associated with magnetite (Fig. 4a – c). B-type veins are characterized by subhedral pyrite with small fractures that is commonly associated with platy or tabular molybdenite (Fig. 4d). Chalcopyrite is rare and unevenly distributed in the quartz gangue. Native Au occurs as micrometer-sized (< 2 μm) inclusions in disseminated pyrite, within the gangue quartz, and/or associated with the sulfosalts (Fig. 4e). Rheniite (ReS2) occurs as fibrous varieties or inclusions along cleavages or grain boundaries of larger platy or tabular molybdenite (Fig. 4f). D-type veins host pyrite that shows massive to cataclastic textures, which is associated with bladed to rosette-like molybdenite that lacks rheniite inclusions (Fig. 4g – i). E1-type veins contain pyrite that either forms euhedral or porous varieties that are partly inclusion-rich (Fig. 4j). Subhedral zoned pyrite is associated with arsenopyrite, galena, sphalerite, and sulfosalts in the E2-type veins (Fig. 4k, l).

Mineral chemistry

Quartz composition

Hydrothermal quartz from the different vein types show distinct mean and median trace element compositions. Arsenic and Sb contents increase systematically from the early porphyry veins towards the late epithermal veins, whereas the Ti contents decrease (Fig. 6a – c). The occurrence of rutile/titanite inclusions in hydrothermal quartz and the host rocks suggest a TiO2 activity of ~ 1, allow the estimation of quartz formation temperatures by Ti-in-quartz thermometry (Zhang et al. 2020). The calculated temperatures (Fig. 6d) show a systematic decrease from ~ 470 °C in the A-type veins to ~ 240 °C in the epithermal veins (Table 2). Quartz from different vein-types can be distinguished by distinct Al/Ti and Ge/Ti ratios (Fig. 7) (Rottier and Casanova 2021). The trace element chemistry of quartz correlates with variations in CL intensity, as indicated by decreasing CL intensities with increasing Ti, As, and Sb contents (Fig. 6a – c).

Fig. 6
figure 6

Box and whisker plots of LA-ICP-MS trace element data of hydrothermal quartz from the different mineralization stages at the Maronia porphyry. Ti-in-quartz temperatures in (d) are calculated based on the fluid inclusion entrapment pressure of the respective vein generations (Table 2), a TiO2 activity of 1, and the calibration from Zhang et al. (2020). Note For statistical presentation in the box and whisker plots values below the detection limits were replaced by values of half of the detection limit (Dmitrijeva et al. 2020)

Table 2 Fluid inclusion homogenization temperatures, melting temperatures, salinities, and entrapment pressures of the different hydrothermal vein types from the Maronia porphyry derived by fluid inclusion microthermometry
Fig. 7
figure 7

Bivariate classification diagram (Ge/Ti vs. Al/Ti) for different hydrothermal porphyry-epithermal vein types at Maronia based on Rottier and Casanova (2021) which verifies the textural and mineralogical vein classification (after Sillitoe 2010) of the investigated samples

Pyrite and arsenopyrite composition

Pyrite from the different vein-types shows increasing As contents from As-poor pyrite in the porphyry veins (A-, B-, D-type) towards arsenian pyrite in the E1-type and arsenopyrite in the E2-type veins (Fig. 8a). Trace elements in pyrite exhibit distinct mean and median values between the porphyry and epithermal stage and various elements like, As, Sb, Au, Ag, Te, Cd, Zn, Tl, Ni, and Pb increase from the porphyry to the epithermal stage (Fig. 8b – e). By contrast, Re is generally enriched in pyrite from the porphyry stage with the highest Re contents reported in pyrite from the B-type veins (2.3 ± 2.9 ppm, 1σ, n = 52; Fig. 8f). Arsenopyrite in the E2-type veins is enriched in most trace elements, except Se, Co, Ni, and Te relative to pyrite from the other vein-types. In addition, epithermal pyrite and arsenopyrite are characterized by higher As/Sb (> 100) and lower Se/Ge (< 100), Se/Te (< 50), and Co/As (< 1) values relative to pyrite from the porphyry stage (Fig. 8g, h, and Fig. 9a – d).

Fig. 8
figure 8

Box and whisker plots of LA-ICP-MS trace element data of hydrothermal pyrite (A-, B-, D- and E1-type) and arsenopyrite (E2-type) from the different mineralization stages at the Maronia porphyry. Temperature and salinity data in (a, g, h) are derived from fluid inclusions from the respective vein generations (Table 2). Note For statistical presentation in the box and whisker plots values below the detection limits were replaced by values of half of the detection limit (Dmitrijeva et al. 2020)

Fig. 9
figure 9

Bivariate plots of LA-ICP-MS trace element ratios of hydrothermal pyrite (A-, B-, D- and E1-type) and arsenopyrite (E2-type) from the different mineralization stages at the Maronia porphyry. The reference data for porphyry pyrites are from Dexing (Reich et al. 2013), Aqua Rica (Franchini et al. 2015) and Koloula (Keith et al. 2022). Reference data for high sulfidation epithermal systems are from Sangihe Island (King et al. 2014) and Agua Rica (Franchini et al. 2015) and for low sulfidation epithermal systems are from Vatukoula (Börner et al. 2021). The temperature and salinities are based on fluid inclusion data (Table 2)

Molybdenite composition

The Mo (32.1 to 61.1 wt%) and Re (0.06 and 23.9 wt%) contents in molybdenite from EPMA of the B- and D-type veins exhibit a strong negative correlation (R2 = 0.76; Fig. 10a). The average LA-ICP-MS Re contents in molybdenite from Maronia are significantly higher (3918 ± 3854 ppm, 1σ, n = 81) compared to the global Re average in molybdenite from other porphyry systems worldwide (779 ± 2581 ppm, 1σ, Re, n = 1718; Fig. 10b; Barton et al. 2020). Rhenium contents in single molybdenite grains from B- and D-type veins are variable (Fig. 4f), but based on an unpaired two-tailed t-test the mean Re content in molybdenite from the B-type veins (6631 ± 4308 ppm, 1σ, n = 36) is statistically distinct (p < 0.05) compared to molybdenite from the D-type veins (Figs. 10c and 1746 ± 1184 ppm, 1σ, n = 45). Lead (Fig. 10d) and Tl (not shown) are enriched in molybdenite from the D-type relative to the B-type veins.

Fig. 10
figure 10

Bivariate plots of EPMA (a) and box and whisker plots of LA-ICP-MS (bd) trace element data of hydrothermal molybdenite from the different mineralization stages at the Maronia porphyry. The global mean Re content of molybdenite from different porphyry style deposits in (b) and (c) are derived from Barton et al. (2020). Temperature and salinity in (c and d) are derived from fluid inclusions in the respective vein type (Table 2). The δ34S values in (c) are from molybdenite separates in the respective vein types (Table 3). Note For statistical presentation in the box and whisker plots values below the detection limits were replaced by values of half of the detection limit (Dmitrijeva et al. 2020)

Table 3 Compilation of δ34S isotope data of in situ measurements of pyrite and of sulfide separates from pyrite, molybdenite, and arsenopyrite

Fluid inclusion petrography and microthermometry

Quartz-hosted fluid inclusion data were collected from A-, B-, D- and E-type veins (Table 2). Three types of primary fluid inclusions were observed: (1) Type 1 inclusions are two-phase (L-V), low-salinity (< 0.7 wt% NaCl equiv.), liquid-rich inclusions that homogenize during heating into the liquid phase by the disappearance of the vapor phase (occupying 10 vol%, Fig. 11a), (2) Type 2 inclusions are multiphase (L-V-S), hypersaline (47.2–34.5 wt% NaCl equiv.), and liquid-rich that homogenize into the liquid phase by the disappearance of the solid phases (mainly halite ± sylvite) and then of the vapor phase (occupying 20 vol% of the inclusion, Fig. 11b). Some inclusions that homogenized by halite dissolution may have been affected by post-entrapment modification (Audétat 2023; Audétat and Günther 1999) and were excluded from further microthermometry. Anhydrite, hematite, and chalcopyrite are common daughter minerals in the type 2 inclusions from the A- and B-type veins (Fig. 11b) but are rare in the type 2 fluid inclusions of the D- and E-type veins. These solid phases did not dissolve during heating up to 600 °C. (3) Type 3 inclusions are two-phase (V-L), vapor-rich (occupying 80 vol%), and moderately saline (20.4–8.5 wt% NaCl equiv.) that homogenize into the vapor phase by the disappearance of the liquid (Fig. 11c).

Fig. 11
figure 11

(ac) Representative photographs of quartz-hosted type 1 (L-V), type 2 (L-V-S) and type 3 (V-L) fluid inclusions assemblages in the porphyry and epithermal veins at the Maronia porphyry (Table 2). (d) Histogram showing the homogenization temperatures of the fluid inclusion types from the different vein types. (e) Homogenization temperature vs. salinity plot of the different fluid inclusion types from the different vein types. NaCl saturation curve (halite + liquid + vapor curve) and critical curve of the systems H2O-NaCl are also shown. (f) Pressure-temperature cross sections of the NaCl-H2O phase diagram illustrating the evolution pathways of the magmatic-hydrothermal fluids at Maronia. The parental pluton intrusion depth is based on amphibole thermobarometry and the subsequent shallower emplacement of the porphyritic microgranite is inferred from field relations (Schaarschmidt et al. 2021b)

Two discrete fluid inclusion assemblages (FIAs) were distinguished: (1) FIA-1 consisting of liquid-rich inclusions with a vapor bubble (type 1), which only occur in the E-type veins (Fig. 11a). (2) FIA-2 with coexisting aqueous multiphase hypersaline liquid-rich (L-V-S, type 2) and vapor-rich low-density (V-L, type 3) inclusions occurring in the A-, B-, and D-type veins (Fig. 11c). The fluid inclusions of FIA-2 have a negative-crystal shape, but also exhibit round, elliptical, or irregular shapes, similarly to those of FIA-1. Their maximum diameter is 74 μm, with the majority ranging from 10 to 20 μm. Mean homogenization temperatures systematically decrease from the A- (~ 489 °C), B- (~ 396 °C), D- (~ 344 °C) towards the E-type veins (~ 282 °C; Table 2). The homogenization temperatures of the coexisting type 2 and 3 inclusions in the same sample are comparable. Eutectic temperatures of the type 1 and 3 inclusions (-21.2 to -23.5 °C) suggest the presence of NaCl and KCl in the fluids. Fluid salinities range from 47.2 to 18.2 wt% NaCl equiv. in the A-type and 43.0 to 8.5 wt% NaCl equiv. in the B- and D-type veins (Fig. 11e; Table 2). Type 1 fluid inclusions in the E-type veins generally show the lowest salinities (< 0.7 wt% NaCl equiv.). Mean fluid minimum entrapment pressures for coexisting type 2 and 3 fluid inclusions decrease from the A- (~ 389 bar), B- (~ 227 bar), D- (~ 118 bar) towards the E-type veins (~ 61 bar, Fig. 11f; Table 2). The absence of liquid CO2 and/or clathrate melting suggests only insignificant contents of CO2 (< 3.5 mol %), which agrees with previous Raman spectrometry results from fluid inclusions at Maronia (Melfos et al. 2002).

Sulfur isotopes of pyrite, molybdenite, and arsenopyrite

Pyrite separates show similar positive δ34S values in the A- (~ 3.4‰) and D-type (~ 3.3‰) veins, which are higher than those of the B-type (~ 2.6‰) and lower than those from the E1-type veins (~ 4.6‰) (Table 3). The arsenopyrite separates from the E2-type veins have a δ34S composition (~ 4.3‰) comparable to pyrite from the E1-type veins (Fig. 12a). The δ34S composition of the molybdenite separates from the B-type (~ 3.2‰) and D-type veins (~ 3.3‰) are similar to the pyrite separates from the same vein-types (Fig. 12a; Table 3).

Fig. 12
figure 12

(a) Compilation of δ34S data of sulfide separates and in situ δ34S LA-ICP-MS measurements at the Maronia porphyry (Table 3) compared to previously published δ34S value for pyrite and molybdenite from Maronia (Melfos et al. 2002), other porphyry systems (Hutchison et al. 2020), and different sulfur reservoirs such as arc melt inclusions (Gurenko 2021; Kawaguchi et al. 2022), arc lavas (Alt et al. 1993; de Hoog et al. 2001; Ueda and Sakai 1984; Woodhead et al. 1987), MORB (Labidi et al. 2012, 2014), meteoric water(Geyh et al. 2000) and seawater (Tostevin et al. 2014). (b) Box-and-whisker diagram of in situ δ34S values for pyrite from the porphyry- and epithermal veins. (c) SO2 disproportionation model adapted fromHutchison et al.(2020) showing the expected δ34SH2S trajectory A for the δ34SSO2 of 11.0‰, trajectory B for an initial δ34SSO2 of 5.0‰, and trajectory C for the evolution towards higher SO4/H2S ratios during the fluid evolution to more reducing conditions (Hutchison et al. 2020)

The in situ δ34S composition of pyrite varies between − 4.1 and 14.8‰ with an average of 4.5 ± 3.3‰. A much larger δ34S variance was observed by the in situ measurements compared to the data of the pyrite separates (Fig. 12a, b). Based on an unpaired two-tailed t-test the mean in situ δ34S value (2.9 ± 2.5‰, 1σ, n = 49) of pyrite from the B-type veins is significantly (p < 0.05) lower than in the A- (5.6 ± 4.1‰, 1σ, n = 27), D- (6.0 ± 3.7‰, 1σ, n = 19), and E1-type (5.1 ± 2.4‰, 1σ, n = 26) veins, which however are statistically similar (p > 0.05). The mean δ34S value of the pyrite separates and the in situ measurements are comparable for the B- and E1-type veins. By contrast, pyrite separates from the A- and D-type veins display about 2–3‰ lower mean δ34S values than the in situ data (Fig. 12a, b). The intra-grain (i.e., core-rim) δ34S variability in pyrite systematically decreases from the early A-type to the late E1-type veins (Table 3).

Discussion

Paragenetic framework and fluid evolution

The various superimposed alteration zones and the associated mineralized vein-types provide a framework for the temporal and spatial evolution of the Maronia porphyry system and its epithermal transition (Fig. 5). The sericitic overprint on the potassic alteration at Maronia indicates telescoping of the system (Melfos et al. 2002, 2020) in response to cooling of the magma body and the downward retreat of the lithostatic-hydrostatic transition. Textural evidence, such as crosscutting relations and change from sinuous A-type to straight B-type veins (Figs. 2 and 3), suggest a consecutive formation and the change from ductile and lithostatic to brittle and sub-lithostatic conditions (Monecke et al. 2018; Sillitoe 2010). Fractured and cataclastic pyrite in the D-type veins (Fig. 4h) is consistent with brittle hydrostatic formation conditions, which also resulted in the reopening of the B-type veins and the precipitation of D-type CL-dark quartz along their central suture line (Fig. 3g). The different CL intensities and textures of hydrothermal quartz (Fig. 3) with distinct trace element contents and ratios (Figs. 6 and 7) support consecutive vein generations (Monecke et al. 2018; Rottier and Casanova 2021). The increasing As and Sb contents in quartz from the early- to late-stage veins (Fig. 6) combined with the occurrence of As-Sb sulfosalts and arsenopyrite in the E-type veins (Table 1; Fig. 4j, k) is consistent with an evolution from a As- and Sb-poor fluid in the porphyry stage to a As- and Sb-rich fluid at epithermal conditions, as also observed by the pyrite chemistry (Fig. 8) (Heinrich 2005; Rottier and Casanova 2021). The association of the different vein types with distinct alteration assemblages (Fig. 5) record decreasing fluid pH conditions during the different veining stages, as reflected by the change from feldspar- (potassic, pH ~ 7) to sericite- (sericitic, pH = 4.5–7) and illite/kaolinite-stable (argillic, pH = 4–6) alteration (Corbett and Leach 1998; Hedenquist and Arribas 2022).

Titanium is the main CL activator in quartz and the correlation of Ti and CL intensity can be used as a quantitative measure for the quartz precipitation temperature using different Ti-in-quartz thermometers (e.g., Huang and Audétat 2012; Zhang et al. 2020). Based on the abundant titanite and rutile inclusions in quartz and their occurrence in the altered host rocks (Table 1; electronic supplement; Melfos et al. 2020; Schaarschmidt et al. 2021b), we estimate a TiO2 activity of ~ 1 for our calculations. We use the Ti-in-quartz calibration by Zhang et al. (2020) and suggest that potential effects of fast quartz precipitation rates on the temperature estimates were minor at Maronia, as indicated by a good agreement of the temperatures with those derived from fluid inclusion micro-thermometry (Table 2). The calculated temperatures are also in the temperature range expected for the different vein types in the different alteration zones of porphyry-epithermal systems (Monecke et al. 2018; Sillitoe 2010). Thus, the CL-bright mosaic quartz in the early A-type veins is suggested to have been formed at temperatures of ~ 470 °C, followed by a continuous decrease in temperature, resulting in the precipitation of the CL-black quartz in the late E-type veins at ~ 240 °C (Fig. 6d; Table 2).

The co-existence of type 2 liquid-rich hypersaline and type 3 vapor-rich moderately saline fluid inclusions in quartz of the A-, B-, and D-type veins are indicative of fluid phase separation (Table 2; Fig. 11c, e) and/or the exsolution of magmatic fluids in the two-phase field at shallow to intermediate depth (< 4 km; Fig. 11f; Audétat 2019; Monecke et al. 2018). As the type 2 and 3 fluid inclusions were entrapped during phase separation, the homogenization temperatures do not require any pressure correction (Bodnar et al. 2014; Roedder 1984), which systematically decrease from ~ 500 °C in the A-type veins to ~ 350 °C in the D-type veins (Table 2). The type 1 fluid inclusions in the E-type veins indicate mixing of an ascending low-salinity vapor with meteoric water, resulting in the formation of a low-salinity fluid at ~ 280 °C (< 0.7 wt% NaCl equiv., Table 2; Fig. 11e). The entrapment pressures at lithostatic conditions for the A- and B-type veins and at hydrostatic pressures for the D- and E-type veins are consistent with a decreasing mineralization depth from < 2.2 to ~ 0.5 km (Fig. 11f). The spatial proximity (< 500 m) between the early A- and B-type veins and the late E2-type veins that host a low- to intermediate-sulfidation mineral assemblage of pyrite + arsenopyrite + sphalerite + galena ± sulfosalt (Einaudi et al. 2003) indicate that epithermal fluid conditions are not only restricted to distal areas of the porphyry intrusion (Fig. 1b; Hedenquist and Lowenstern 1994). We interpret this to reflect the downward retreat of the crystallization front in the porphyry magma body, also known from other porphyry systems (Kouzmanov and Pokrovski 2012; Sillitoe 1994, 2010). The epithermal overprint of the porphyry system is accompanied by fluid cooling and the influx of meteoric fluids (Sillitoe 1994; Simmons et al. 2005), as indicated by the type 1 fluid inclusions at Maronia, possibly induced along structural pathways in proximity to the porphyry stock (Fig. 1b) (Melfos et al. 2020).

Summarizing, the formation of the Maronia Cu-Mo ± Re ± Au porphyry-epithermal mineralization occurs in four consecutive hydrothermal stages (Fig. 5), defined by changes in physicochemical fluid composition (e.g., temperature, salinity, pH) induced by phase separation and meteoric water influx during telescoping, contemporaneous uplift, and cooling of the porphyritic microgranite (Fig. 13).

Fig. 13
figure 13

Simplified model (not to scale) for the hydrothermal evolution of the Maronia porphyry systems as reflected by the magmatic evolution (Schaarschmidt et al. 2021b), alteration assemblages (Melfos et al. 2002, 2020), fluid inclusion data, sulfur isotope data of molybdenite, arsenopyrite, and pyrite as well as the trace element composition of hydrothermal quartz, molybdenite, pyrite and arsenopyrite. Changes of the vein mineralogy are presented in Table 1. (a) Early porphyry stage 1: Formation of potassic alteration as well as A-type veins form high salinity liquids during phase separation of the primary magmatic fluids. (b) Intermediate Porphyry stage 2: Formation of potassic alteration and precipitation of Re-rich molybdenite and rheniite, as well as native Au in B-type veins from high salinity oxidized liquids during phase separation with the lowest average δ34S composition at Maronia. (c) Late porphyry stage 3: Lithostatic-hydrostatic transition during telescoping of the porphyry systems by cooling of the intrusion. Formation of the sericitic alteration and relatively Re-poor molybdenite in D-type veins by intermediate temperature fluids during phase separation. (d) Epithermal stage 4: Waning magmatic activity results in cooling and mixing of the magmatic-hydrothermal fluids with metalliferous meteoric water causing argillic alteration and silicification of the host rocks, as well as the formation of E-type veins. Modified after Heinrich (2005), Sillitoe (2010), Kouzmanov and Pokrovski (2012) and Keith et al. (2022). Note: Dashed lines represent inferred and/or weak fluid contributions

Tracing phase separation and metal sources by pyrite chemistry

Pyrite from the porphyry stage (A-, B-, and D-type veins) at Maronia is depleted in elements such as As, Au, Te, Ag, Sb, Tl, Zn, and Pb relative to pyrite from the epithermal stage (E1- and E2-type veins), which is consistent with the trace element distribution of pyrite from other porphyry and epithermal environments (Franchini et al. 2015; Keith et al. 2022; Maydagán et al. 2013; Sykora et al. 2018). This can be linked to changes in fluid temperature, as these elements typically concentrate in pyrite that forms from fluids at lower temperatures (< 350 °C) such as at epithermal conditions (Falkenberg et al. 2022; Maydagán et al. 2015; Metz and Trefry 2000; Sykora et al. 2018). The Se/Te ratio in pyrite records fluid cooling (Keith et al. 2022), and verifies the transition from high-temperature porphyry (> 350 °C, Se/Te > 50) to low-temperature epithermal (< 350 °C, Se/Te < 50) conditions at Maronia (Figs. 8h and 9a). In addition, the Se/Ge ratio in pyrite has been suggested to be a tracer for mixing of different fluids, such as magmatic (Se/Ge ~ 100) and meteoric (Se/Ge < 100) fluids. The chemical similarity of Se and Ge indicates that hydrothermal fractionation processes like phase separation or cooling has only a negligible effect on the Se/Ge ratios (Brugger et al. 2016; Keith et al. 2022; Klose et al. 2021). At Maronia the Se/Ge ratios are consistently high (128 ± 36, 1σ, n = 69) in the A-, B-, and D-type veins (Fig. 9c), which show different fluid inclusion homogenization temperatures (Table 2), thus confirming that the Se/Ge ratio is independent of fluid temperature. By contrast, the Se/Te ratio varies over almost two orders of magnitude due to its strong temperature sensitivity (Fig. 8h). Hence, porphyry pyrite from Maronia precipitated from a magmatic fluid resulting in its high Se/Ge ratios, while the lower Se/Ge (35 ± 53, 1σ, n = 24) ratios in pyrite from the E-type veins implies that the magma-derived fluid mixed with variable proportions of meteoric water (Fig. 9c). This is supported by the associated type 1 fluid inclusions and the comparable Se/Ge ratios in pyrite from the E-type veins at Maronia and from other meteoric- and/or seawater influenced systems such as the Vatukoula Au-Te deposit, Fiji (Fig. 9c, Börner et al. 2021).

In addition to fluid cooling, phase separation is an omnipresent ore-forming process in porphyry-epithermal environments (Melfos et al. 2002; Pokrovski et al. 2013; Reed and Palandri 2006). Phase separation is recorded in pyrite by ratios of trace elements with different vapor-liquid partition coefficients (Keith et al. 2022; Román et al. 2019; Schaarschmidt et al. 2021a). For instance, As (Kdvapor−liquid = 0.03–9) has a stronger affinity to the vapor phase than Sb (Kdvapor−liquid = ~ 0.09) during vapor-liquid separation, whereas Cl-complexed elements like Co preferentially fractionate into the Cl-rich liquid phase (Brugger et al. 2016; Kouzmanov and Pokrovski 2012; Pokrovski et al. 2013). This results in higher As/Sb and lower Co/As ratios in vapor-rich fluids and their respective precipitates relative to the liquid phase and the original fluid (Keith et al. 2022; Schaarschmidt et al. 2021a). Thus, the negative correlation between Co/As and As/Sb reflects changes in fluid salinity related to vapor-liquid fractionation during phase separation (Fig. 9b). This indicates that the lower As/Sb (< 100) and higher Co/As (> 1) values in the porphyry pyrite from Maronia are related to precipitation from hypersaline liquids (Fig. 9b), as recorded by the type 2 fluid inclusions. By contrast, the high As/Sb (> 100) and low Co/As (< 1) values in pyrite from the E-type veins (Fig. 9b) suggest a formation from a low-salinity vapor-rich fluid. Hence, volatile elements, such as As, Sb, Te, and Tl (± Au) that are concentrated in pyrite of the E-type veins were mostly derived from vapor-rich fluids that mixed with meteoric waters, as indicated by their low Se/Ge values (< 100, Fig. 8c) and the related low-salinity type 1 fluid inclusion of the E-type veins (Cooke and McPhail 2001; Grundler et al. 2013; Pokrovski et al. 2013, 2014). Accordingly, Cl-complexed and non-volatile elements like Pb and Zn (Brugger et al. 2016; Pokrovski et al. 2013) must instead be sourced from metalliferous meteoric waters (Börner et al. 2022; Keith et al. 2022; Simmons et al. 2005). A similar process is likely responsible for the Au deposition in the epithermal stage at Maronia, as the Au contents in pyrite increase with decreasing Se/Ge ratios (Fig. 9d). By contrast, the occurrence of native Au in the high-temperature (~ 400 °C) porphyry B-type veins in the potassic alteration (Fig. 4e), and bulk ore contents of up to 12 ppm Au in the sericitic alteration (Melfos et al. 2020), cannot easily be explained by such a mixing process. Gold occurs as an Au(HS)2 complex in S-rich fluids at temperatures between 250 and 450 °C, similar to those observed in the B- and D-veining stages at Maronia. At temperatures > 450 °C AuCl2 complexes would possibly dominate (Pokrovski et al. 2014). More importantly, the Au(HS)2 species shows a solubility minimum at ~ 400 °C (Pokrovski et al. 2014), which is similar to the fluid inclusion homogenization temperatures of the B-type veins at Maronia (Table 2). Hence, we conclude that the native Au formation was probably triggered due to the low Au solubility combined with the additional destabilization of the Au(HS)2 complex by H2S loss to the vapor phase and to sulfide minerals during intense host rock sulfidation. Based on the integrated study of quartz-hosted fluid inclusions and genetically related pyrite, we conclude that trace element ratios in pyrite record changes in fluid conditions, such as cooling (Se/Te), phase separation (Co/As, As/Sb), and mixing of magmatic and meteoric water (Se/Ge), whereas fluid mixing primarily controls the epithermal Pb, Zn, and possibly the Au precipitation at Maronia (Fig. 13).

Sulfur sources and S isotope fractionation

Hydrothermal sulfides precipitate from the H2S component of fluids, where the related δ34SH2S composition is controlled by the magmatic sulfur source (Hutchison et al. 2020; McDermott et al. 2015), isotope fractionation between SO42− and H2S during the magmatic SO2 disproportionation (Kusakabe et al. 2000; Marini et al. 2011), boiling-induced isotope fractionation during the oxidation of H2S to SO42− (Börner et al. 2022; McKibben and Eldridge 1990; Schaarschmidt et al. 2021a), and mixing of S derived from different reservoirs with distinct δ34S composition (e.g., meteoric water, sediments, magmatic sulfur; Marini et al. 2011). The A-type vein pyrite separates that formed from the high-temperature magma-derived fluids have a mean δ34S value of ~ 3.4‰. Hence, the δ34SH2S value of the magmatic fluid at Maronia can be estimated using the established formation temperatures (~ 500 °C) and experimentally derived temperature-dependent equilibrium fractionation coefficients between FeS2 and H2S (αH2S−FeS2 = 1.0007; Ohmoto 1972), which yielded an initial δ34SH2S value of ~ 2.7‰. The δ34S values of the pyrite separates were preferentially used in this context instead of the in situ data, since they are representative for the fluid conditions and sources of the entire vein stage. However, the initial δ34SH2S values strongly depend on the initial δ34SSO2 composition of the magma-derived fluid, as indicated by S isotope fractionation during SO2 disproportionation (Marini et al. 2011):

$$4S{O}_{2}\hspace{0.17em}+\hspace{0.17em}4{H}_{2}O \leftrightarrow {H}_{2}S\hspace{0.17em}+\hspace{0.17em}3HS{O}_{{4}^{-}} + 3{H}^{+}$$
(1)

Up to 99% of the SO2 dissolved in a magma-derived fluid will disproportionate, as the fluid cools from > 700 °C to ~ 400 °C (Kouzmanov and Pokrovski 2012) and the δ34S composition of the resulting H2S component of the fluid will significantly decrease (Kusakabe et al. 2000; Ohmoto and Lasaga 1982). In order to achieve the proposed initial δ34SH2S value of ~ 2.7‰ by SO2 disproportionation at 500 °C and at a SO42−/H2S ratio of ~ 3 according to Eq. (1), a δ34SSO2 value of ~ 11.0‰ must be assumed (Fig. 12c, trajectory A). However, the expected range for the magmatic δ34SSO2 value is rather between 0 and 5‰, as demonstrated by the isotope fractionation models (Hutchison et al. 2020; McDermott et al. (2015) and the bulk S isotope composition (δ34S∑S = ~ 5.0‰) of calc-alkaline volcanic rocks and melt inclusions (e.g., Kawaguchi et al. 2022; Marini et al. 2011). Hence, our results for the initial δ34SSO2 value of the magma-derived fluid at Maronia is unrealistically high. Considering an initial magmatic δ34SSO2 value of ~ 5.0‰ and SO2 disproportion and a SO42−/H2S ratio of ~ 3, yields δ34SH2S values of about − 3.4‰ at 500 °C, -5.8‰ at 400 °C, and − 9.6 at 300 °C according to the δ34SH2S−SO2 fractionation model by Hutchison et al. (2020) (Fig. 12c, trajectory B). These temperatures overlap with common fluid temperatures of potassic, sericitic, and argillic alteration zones (Sillitoe 2010), and with our fluid inclusion temperature constraints (Table 2). Thus, with respect to the temperature-dependent H2S-FeS2 isotope fractionation (Ohmoto 1972; Ohmoto and Lasaga 1982), pyrite that formed from the H2S portion after SO2 disproportionation from a magma-derived fluid should have a δ34S value of -2.7‰ at 500 °C, -5.1‰ at 400 °C, and − 8.4 at 300 °C. The δ34S values of the pyrite separates are generally in the positive range and only few of the lowest δ34S values from the in situ δ34S data reach comparably low values (Fig. 12a; Table 3). Therefore, the influence of SO2 disproportionation on the δ34S composition of pyrite from Maronia was minor and it must have been controlled by other processes.

Alternately, phase separation can lead to strong S isotope fractionation due to the preferential partitioning of 34S into SO42− during H2S oxidation (Börner et al. 2022; McKibben and Eldridge 1990; Ohmoto and Lasaga 1982). The FIAs from the A-, B-, and D-type veins at Maronia provide evidence for phase separation (Fig. 11), as also supported by the trace element data in pyrite (Fig. 9). In addition, phase separation can lead to significant δ34S variations in individual pyrite grains on the micro-scale (Börner et al. 2022; McKibben and Eldridge 1990; Schaarschmidt et al. 2021a). The significant intra-grain core to rim δ34S variability in pyrite which decreases from the A- (up to 9.5‰), B- (up to 8.7‰), to the D-type (up to 6.8‰) veins may therefore record the effect of fluid phase separation on the S isotope system, as also supported by the fluid inclusion and pyrite trace element constraints, but cannot explain the general positive and high δ34S values. Instead, the Se/Ge values (< 100) in pyrite from the E-type veins and the related type 1 fluid inclusions provide evidence for mixing of the magma-derived fluid with meteoric water, which has variable, but also high δ34S values (-5 to 15‰; Fig. 12; Geyh et al. 2000), explaining the positive δ34S values of pyrite from the E-type and possibly D-type veins. However, a meteoric water influence during the formation of the A- and B-type veins that were emplaced at lithostatic pressures in the high-temperature porphyry stage is unlikely and cannot cause the positive δ34S values reported in the related pyrites. However, intense fluid-rock interaction during the potassic alteration can buffer the δ34SH2S of the fluid and its precipitates to positive δ34S values similar to those of the host rocks (Fig. 12a), as indicated by the δ34S composition of arc lavas (5.0 ± 3.9‰, 1σ, n = 125; Ueda and Sakai 1984; Woodhead et al. 1987; Alt et al. 1993; de Hoog et al. 2001) and silicate mineral hosted melt inclusions in arc lavas (6.0 ± 2.5‰, 1σ, n = 141; Gurenko 2021; Kawaguchi et al. 2022). In addition, decreasing fluid fO2 and pH as indicated by the change from potassic to sericitic and argillic alteration may also be caused by the reduction of aqueous SO42− to H2S via oxidation of Fe2+ in the associated host rocks and/or fluid to Fe3+, which would stabilize the reduced S species (Hutchison et al. 2020; Sulaksono et al. 2021). This would reduce the fluid SO42−/H2S ratio and also shift the δ34SH2S to higher values after SO2 disproportionation (Fig. 12c, trajectory C; Hutchison et al. 2020). Therefore, we summarize that the expected negative δ34S signature of the sulfides at Maronia, that formed following magmatic SO2 disproportionation and phase separation, are overprinted by fluid-rock interaction and mixing with meteoric water resulting in a fluid with a relatively low SO42−/H2S ratio and high δ34SH2S value (Figs. 12 and 13).

Hydrothermal controls on Re enrichment on the deposit scale

The high Re contents in molybdenite (6630 ± 4310 ppm Re, 1σ, n = 36) and associated rheniite in the B-type veins at Maronia (Fig. 4f) suggest that the related porphyry fluid was particularly enriched in Re. By contrast, molybdenite in the D-type veins lacks rheniite and has lower Re contents (1750 ± 1180 ppm Re, 1σ, n = 45; Fig. 10c) implying a lower Re budget in the later stage fluid, which is also recorded by the lower Re contents in pyrite (Fig. 8f). Irrespective of the vein-type, molybdenite from Maronia and other porphyries in NE Greece (Berzina et al. 2005; Voudouris et al. 2009, 2013a) is significantly enriched in Re compared to the global Re average in porphyry molybdenite (Figs. 10b and 779 ± 2581 ppm, 1σ, Re, n = 1718, Barton et al. 2020). This is evidence for a highly efficient Re precipitation process from Re-rich porphyry fluids, which may have been exsolved from a region-scale Re-rich magmatic systems possibly linked to the contribution of subducted sediments to the magmatic evolution (Sun et al. 2004; Xue et al. 2023).

At magmatic-hydrothermal supercritical conditions, Re is mainly transported as ReCl4 in Cl-rich fluids at 400 °C and as ReCl3+ at 500 °C (Xiong and Wood 1999). By contrast, at subcritical conditions, such as in shallow volcanic environments, vapor transport of HReO4 or Re(OH)4 is possible (Xiong et al. 2006), as evidenced by extremely Re-rich molybdenite and rheniite sublimation from volcanic gases (Korzhinsky et al. 1994; Tessalina et al. 2008). However, the elevated pressure conditions during the formation of B- and D-type veins at Maronia makes a sublimation related formation of molybdenite and rheniite from a vapor-rich fluid unlikely. Furthermore, cogenetic pyrite of the B- and D-type veins shows trace element signatures indicative of sulfide precipitation from a hypersaline liquid phase, as also supported by the fluid inclusion data (Table 2). Importantly, the E-type veins, which are related to a vapor-rich fluid, show the lowest Re contents in pyrite (Fig. 8f). Hence, we conclude that vapor-rich fluids did not play an important role in the transport of Re from source to sink at Maronia. Instead, it seems to be likely that most Re was transported as a Cl-complex (i.e., ReCl4 or ReCl3+) (Xiong et al. 2006). In sulfide-rich fluids, the Re solubility is strongly controlled by fO2 and Re is essentially immobile under reduced conditions (Barton et al. 2020; Xiong et al. 2006). Thus, large amounts of Re can be transported by an oxidized and Cl-rich fluid (Xiong et al. 2006), similarly to the fluid which was possibly responsible for the potassic alteration at Maronia. The high fO2 of the fluid is indicated by an oxidized alteration mineral assemblage consisting of magnetite, low Fe/(Fe + Mg) biotite, and low Fe/Al titanite (Melfos et al. 2020; Schaarschmidt et al. 2021b), as well as the common occurrence of oxidized daughter minerals (hematite, anhydrite) in the hypersaline type 2 fluid inclusions (Fig. 11b). By contrast, the absence of oxidized daughter minerals in the fluid inclusion assemblages of the subsequent sericitic alteration, suggest more reduced conditions hindering Re mobility and solubility (Einaudi et al. 2003; Seedorff et al. 2005; Sillitoe 2010). As the δ34S composition of molybdenite separates from B- (~ 3.2‰) and D-type (~ 3.3‰) veins are similar, we propose that the fluids were either derived from the same source with comparable S isotope composition and/or that their formation was controlled by similar processes and thus cannot explain the difference in the Re endowment. However, the bulk δ34S data of the molybdenite separates do not take processes into account that are only preserved by intra-grain δ34S variations, as observed in the in situ δ34S values of co-genetic pyrite (Fig. 12a; Table 3). Importantly, the in situ δ34S composition of pyrite from the B-type veins is statistically distinct compared to pyrite from the D-type veins (Fig. 12b), which is indicative of distinct formation conditions. This could either be explained by (1) a higher portion of magmatic H2S derived from SO2 disproportionation and/or (2) higher SO42−/H2S ratios in the fluid, causing higher fO2 conditions in the B-type compared to the D-type veins (Fig. 12c). Following our previous conclusions, we suggest that the involvement of an oxidized magma-derived fluid was key for the high Re budget and the formation of the Re-rich molybdenites during the B-type stage.

A large database approach by Barton et al. (2020) did not reveal any systematic temperature effect on Re enrichment in molybdenite. However, experimental results (e.g., Hurtig and Williams-Jones 2014; Xiong and Wood 2002) and empirical evidence from natural high-temperature porphyry molybdenite (Spencer et al. 2015); Stein et al. 2001; Ye et al. 2022; Zhang et al. 2021) in combination with our detailed paragenetic and geochemical observations from Maronia suggest that fluid cooling can contribute to efficient Re precipitation between ~ 400 °C and 350 °C in the porphyry environment. This is supported by the small but notable temperature differences recorded by the fluid inclusion homogenization temperatures from B-type (~ 400 °C) to D-type (~ 350 °C) veins and by the variance of temperature sensitive trace elements such as Tl and Pb in molybdenite (Fig. 10f, g). Comparable observations were made for the Re-rich Pagoni-Rachi porphyry in NE Greece (Voudouris et al. 2013b) and the SarCheshmeh porphyry Cu-Mo deposit in Iran (Shafieibafti et al. 2021), where Re precipitation was also related to fluid cooling below 400 °C. Thus, we propose that the especially the combination of reduction and cooling of an oxidized, Cl-rich, and Re-rich magmatic fluid can induce the efficient Re precipitation at Maronia, which might also be an important process for other Re-bearing porphyries worldwide.

Conclusions

The consecutive vein generations and alteration zones in the Maronia Cu-Mo ± Re ± Au porphyry system and its epithermal transition provide a well-constrained temporal and spatial framework for the fluid evolution and related Re enrichment in porphyry environments. The presented data indicate a continuous fluid evolution from near-neutral high-temperature (~ 500 °C) conditions towards moderately acidic intermediate temperature (~ 350 °C) fluids accompanied by phase separation, forming the A-, B- and D-type veins in the potassic and sericitic alteration zone. The late-stage epithermal E-type veins that are associated with an argillic alteration assemblage and silicified host rocks, formed by mixing of a vapor-rich fluid with metalliferous meteoric waters at lower fluid temperatures (~ 280 °C). The mostly positive δ34S values in the hydrothermal sulfides from Maronia are the results of intense fluid-rock interaction and mixing with meteoric water, overprinting the initial S isotope signature derived from magmatic SO2 disproportionation and phase separation. Importantly, these complex δ34S systematics are not preserved by the δ34S data of sulfide separates, emphasizing the importance of single grain in situ δ34S analysis for reconstructing ore-forming processes in detail.

At Maronia, we can show that the particularly high Re budget in the magma-derived fluid that forming the B-type veins is linked to their oxidized and Cl-rich nature at relatively high SO42−/H2S values enhancing Re solubility. The combination of reduction and cooling of this Re-rich fluid below ~ 400 °C, is a key process triggering the efficient Re precipitation needed to form Re-rich molybdenite and rheniite in the potassic alteration zone at Maronia. The high Re enrichment in molybdenite at Maronia and in other systems in NE Greece, compared to the global average from other porphyry systems worldwide, implies that in addition to the local deposit-scale hydrothermal processes a larger regional-scale enrichment process, such as the magmatic evolution might be important for Re enrichment, which is, however, a matter of future research.