Introduction

Silicic magmas in continental igneous systems chemically evolve through differentiation and/or contamination before they solidify to form plutons at different crustal levels and/or erupt on the Earth’s surface (e.g., Hildreth 2004; Bachmann et al. 2007; Cashman et al. 2017). The geochemistry of the resulting plutonic and volcanic rocks integrates a whole range of igneous processes that can occur from beneath the Moho to the surface. Hence, constraining parameters such as the physicochemical state of magma bodies at different crustal depths, crustal- and mantle-derived contributions in magmas and timescales of magma differentiation and storage, are instrumental to understand the behaviour of lithospheric magmatic systems and the mechanisms generating large volcanic eruptions and related mineral resources (e.g., Hildreth and Moorbath 1988; Thompson et al. 2002; Caricchi et al. 2014).

In trans-lithospheric magmatic systems emplaced in the continental crust, the role of crustal contribution in both the mantle source regions and the magmatic column during magma ascent and emplacement remains debated. Main discussions concern the processes that produce evolved magmas (see review by Moyen et al. 2021), from nearly pure crystal fractionation of mantle-derived melts in a closed system (e.g., Bowen 1928; Geist et al. 1995) to predominant crustal melting (e.g., Price et al. 2005; Bindeman et al. 2008; Matsumoto and Nakagawa 2010; Zheng and Gao 2021; Clemens et al. 2021), to “hybridization” of mantle- and crustal-derived magmas as intermediate options (see Jacob et al. 2021 and references therein). In addition, the source of silicic magmas (whether dominantly reworked crustal material or fractionated mantle-derived mafic melts) has fundamentally different implications on the heat and mass transfer in the crust (e.g., Moyen et al. 2021), and ultimately on the eruptive behaviour and style.

A major obstacle to understand the genesis of evolved magmas is the scarceness of exposed magma bodies from different crustal levels, making it difficult to assess the compositional variability as well as the temporal and spatial evolution of the entire crustal-scale magmatic column from early intrusions in the lower crust up to erupted products at the surface. To date, extensive geochronological and isotopic data have been obtained for crustal sections representing fossil intra-oceanic subduction zones, notably former island arc systems such as the Kohistan arc and Talkeetna arc (Rioux et al. 2007; Bouilhol et al. 2011; Jagoutz et al. 2019). However, in tectonic settings where the mantle-derived magmas emplaced in pre-existing continental crust, such data are either lacking or were obtained separately from distinct crustal levels and not integrated throughout the entire crustal column (e.g., the intra-continental Famatinian arc in Argentina; Walker et al. 2015). Here, we study the post-collisional Permian (ca. 280–285 Ma) Sesia Magmatic System (southern Ivrea-Verbano Zone, Southern Alps, NW Italy), an ideal natural laboratory to investigate the importance of crustal contribution on the evolution of mantle-derived magmas emplaced in an extending high-grade, late Neoproterozoic to late Paleozoic continental crust (e.g., Pin 1990; Zurbriggen et al. 1997; Schaltegger and Brack 2007).

The question of crustal contamination in the Sesia Magmatic System (SMS) has dominantly been addressed from a whole-rock geochemical perspective (e.g., Voshage et al. 1990; Sinigoi et al. 1991, 2011, 2016; Pinarelli et al. 2002), thermal modelling (e.g., Barboza and Bergantz 1996, 1997, 1998) and recently by garnet xenocrysts (Devoir et al. 2021) so far. However, there was no attempt to reconcile observations from both methods or obtain a self-consistent dataset at the scale of the entire magmatic system, i.e., from both the lower and upper crustal levels. In addition, the lifetime of the entire magmatic system and the composition of the crust at the time of mafic magma emplacement, which are key parameters constraining crustal fertility (e.g., Annen and Sparks 2002; Annen et al. 2006; Karakas and Dufek 2015), are still debated. In this study, we combine (1) a new, self-consistent dataset of zircon Hf isotope measurements in representative samples of both lower mafic and upper silicic igneous rocks that were previously dated by high precision CA-ID-TIMS zircon U–Pb single grain analysis (Karakas et al. 2019); with (2) a two-dimensional thermal model of polybaric mantle-derived magma differentiation and associated crustal melting (modified after Karakas et al. 2017 and Karakas et al. 2019). We use our results, in combination with previously published data, to constrain the evolution of the SMS and clarify the relative contribution of mantle- and crustal sources, and how these vary in space and time throughout the formation of such crustal-scale magmatic systems.

Geologic background

The Ivrea-Verbano Zone and the adjacent Serie dei Laghi (northern Italy, Fig. 1) constitute one of the best-exposed, most complete sections of continental crust on Earth (e.g., Foutain 1976; Handy 1987; Voshage et al. 1990; Bertotti et al. 1993; Schaltegger and Gebauer 1999; Schaltegger and Brack 2007). Field observations and geobarometric studies have revealed increasing pressures from the silicic to mafic units (Henk et al. 1997; Demarchi et al. 1998), suggesting that the entire system represents a tilted coherent, up to 25 km thick crustal-section (Handy et al. 1999).

Fig. 1
figure 1

Geologic map of Ivrea-Sesia crustal-scale magmatic system and sample locations, after Peressini et al. (2007), Brack et al. (2010), Quick et al. (2009) and Karakas et al. (2019). The samples were collected systematically from different igneous units for single grain zircon U–Pb CA-ID-TIMS analysis (Karakas et al. 2019): cumulus body with pyronexite composition (sample GP3); two samples from the Mafic Complex (3 g and 5 g); three samples from the Dioritic unit (IV-01, IV-02, and IV-06), four samples from granitoids (IVZ-03, IVZ-07, IVZ-11, and IVZ-12); one sample from an upper crustal intruded dyke (ROC-A), and a sample from the erupted rhyolite (MBR; Karakas et al. 2019)

Full size image

Within the Ivrea-Verbano Zone, a thick sequence of high-grade metamorphic rocks, including migmatitic metapelite, metapsammite, calc-silicate and amphibolite, is commonly grouped as the Kinzigite Formation (e.g., Schmid 1967; Sills and Tarney 1984). The metamorphic grade in these rocks increases from amphibolite to granulite facies approaching the Insubric line, i.e., the tectonic boundary separating the Ivrea-Verbano Zone from the Central Alps (Fig. 1; Zingg et al. 1990). Slivers of mantle peridotite, including the famous Balmuccia body, were tectonically emplaced in the Kinzigite Formation in pre-Permian time (Quick et al. 1995). During the Permian extensional event, the northern and southern parts of the Kinzigite Formation were intruded by a 8-km thick gabbro body, referred to as the Mafic Complex (MC; Rivalenti 1975; Rivalenti et al. 1984; Voshage et al. 1990; Sinigoi et al. 1991; Sinigoi et al. 1994, 2011). In the south, the latter is subdivided in Upper and Lower Mafic Complex, separated by a zone wherein bands of granulite facies paragneiss are common (the “Paragneiss bearing belt” of Sinigoi et al. 1996). Rocks of the Mafic Complex are dominated by gabbroic cumulates, which show a decreasing intensity of foliation and increase in intercumulus components approaching the top of the body, where the rocks were historically termed “diorites”. The Mafic Complex is characterized by an arcuate structure defined by layering and foliation that is best-explained by the “gabbro glacier” model (Quick et al. 1992, 1994), where the incremental magma intrusion into the gabbro unit over time was accompanied by the downward and outward deformation of the crystallizing material.

The “Serie dei Laghi” constitutes the upper part of the crustal section and consists of lower amphibolite facies gneisses and schists intruded by granite bodies (Boriani et al. 1988; Tavazzani et al. 2017, 2020 and references therein). The three main granite bodies, known as Baveno, Roccapietra and Valle Mosso, are displaced by the Pogallo and Cremosina faults (Fig. 1). The whole section is overlain by a bimodal association of volcanic rocks, dominated by rhyolitic tuffs. The presence and distribution of megabreccia units lead to interpret the volcanic field as a deeply eroded caldera, at least 13 km across (Quick et al. 2009).

The relationships between igneous units at different levels in the crustal section were largely debated in the last 40 years. The discussion involved the possibility that the granites were unrelated to the Mafic Complex, and result from an event which caused granulite-facies metamorphism and melting in the lower crust. This possibility was precluded by SHRIMP U–Pb zircon age data, which demonstrate that the Mafic Complex, the main granite bodies and volcanics formed within the same temporal lapse (Quick et al. 2009). The contemporaneity of MC, granites and volcanic rocks lead to group these igneous units in the “Sesia Magmatic System” (SMS, Sinigoi et al. 2016) narrowed down to less than 4 Myr by more precise, CA-ID-TIMS dating (Karakas et al. 2019).

Despite this, the petrogenetic relationships between them remain unclear. Based on overlapping, Nd isotopic compositions, Voshage et al. (1990) proposed that the granites of the SMS crystallized from melts fractionated from the Mafic Complex. Alternatively, based mainly on Sr isotope compositions, Sinigoi et al. (1994) and Pinarelli et al. (2002) proposed that the granites were produced by hybridization between mantle-derived magmas and crustal melts, possibly produced by heat advection supplied from the Mafic Complex. Likewise, based on a balance of Nd and Sr isotopes and trace elements, Sinigoi et al. (2016) inferred that granitic magmas of the SMS may be a mixture of residual melt fractionated from the Mafic Complex and crustal melts formed atop of the advancing front of the mafic intrusion. The importance of crustal melting in generating silicic magmas in the Sesia Magmatic System was questioned by Barboza et al. (1999) and Barboza and Bergantz (2000), mainly based on the observation that the thermal impact on the roof of the MC is limited to less than 2 km of migmatite, and that the metamorphic gradient in the Kinzigite Formation pre-dates the intrusion of the Mafic Complex. This is in agreement with field observations (e.g., Zingg 1990; Sinigoi et al. 2011) and recent geochronological data on the Kinzigite Formation showing that the onset of crustal melting pre-dates the emplacement of the SMS by >20 Myr (Ewing et al. 2014; Guergouz et al. 2018; Kunz et al. 2018).

Methods

LA-MC-ICP-MS zircon Hf isotope analysis

We measured the Hf isotopic compositions of zircon grains from the 11 samples studied by Karakas et al. (2019) for zircon trace elements and high-precision CA-ID-TIMS U–Pb dating, plus one additional sample from a granitic dyke intrusion. The locations of the samples are marked in Fig. 1. We note that the individual zircon grains used for Hf isotope analysis are different than those used in the ID-TIMS U–Pb analysis of Karakas et al. (2019), although all of the grains belong to the same hand samples. Moreover, the zircons investigated here for Hf isotopes were previously screened for U–Pb ages by LA-ICP-MS (as reported in Karakas et al. 2019), which allowed us to unambiguously identify magmatic vs. inherited/xenocrystic grains. In total, we investigated 6 samples from the lower crustal mafic complex, 5 from upper crustal granitoids and one from the silicic volcanic rocks.

Zircon Lu–Hf isotopic analyses were carried out by laser ablation—multi-collection—inductively coupled plasma—mass spectrometry (LA-MC-ICP-MS) at ETH Zürich. We have used a RESOlution (ASI/Applied Spectra) excimer ArF (193 nm wavelength) laser ablation system equipped with the dual-volume S-155 ablation cell (Laurin Technic), attached to a Nu Plasma 2 (Nu Instruments), multi-collector sector-field mass spectrometer. We used a laser spot size of 50 μm, a repetition rate of 5 Hz and an energy density of ca. 3.5 J∙cm–2. The carrier gas consisted of high-purity helium (ca. 0.35 L∙min−1) and argon sample gas from the MC-ICP-MS (ca. 1.0 L∙min−1). The MC-ICP-MS was optimized for maximum sensitivity on Hf isotopes and peak alignment throughout the investigated mass range. We acquired intensities for the following isotopes (corresponding Faraday cup indicated between brackets): 171Yb (L4), 173Yb (L2), 175Lu (Ax), 176(Yb + Lu + Hf) (H1), 177Hf (H2), 178Hf (H3), 179Hf (H4), 180Hf (H5), 181Ta (H6).

The data were processed offline with the Igor Pro-based Iolite v2.5 software (Hellstrom et al. 2008), using an in-house data reduction scheme. Instrumental mass bias for Yb and Hf isotopes were corrected to the natural abundance ratios 173Yb/171Yb and 179Hf/177Hf and the isobaric interferences of 176Yb and 176Lu on 176Hf were corrected using the natural abundance ratios of 176Yb/173Yb and 176Lu/175Lu, with all natural abundance ratios taken from Chu et al. (2002). The mass bias correction factor obtained for Yb isotopes was applied to Lu isotopes. Accuracy and external reproducibility of the method were controlled by repeated analyses of reference zircon standards GJ-1 (Morel et al. 2008), Plešovice (Slama et al. 2008), 91500 (Wiedenbeck et al. 1995), and Temora (Woodhead and Hergt 2005). The 176Hf/177Hf ratios obtained on all reference materials are within uncertainties identical to the recommended values (Fig. 2A, see also Supplementary Table S1).

Fig. 2
figure 2

Plots of εHf (calculated at crystallization age, i.e., 283 Ma for the samples) as a function of 176Yb/177Hf ratio for all the reference materials (A) and samples (B) measured in this study. In A the reference εHf for the various zircon reference materials were obtained using the recommended 176Hf/177Hf ratios and crystallization ages from Fisher et al. (2014). Error bars correspond to 2 S.E. uncertainties

Full size image

The initial 176Hf/177Hf ratios were calculated for each analysis using a 176Lu decay constant of 1.865∙10–11 (Scherer et al. 2001), the measured 176Lu/177Hf ratio and an average crystallization age of 283 Ma as constrained by high-precision CA-ID-TIMS dating (Karakas et al. 2019). The initial εHf were calculated using the parameters for the chondritic uniform reservoir (CHUR) recommended by Bouvier et al. (2008). The quoted uncertainties on initial isotopic compositions of unknowns include the analytical uncertainty (2 S.E.) and the average intra-session reproducibility (2 S.D.) of initial isotopic compositions of reference materials (typically ~1 εHf unit), propagated by quadratic addition. The data are reported in Supplementary Table S1.

Thermal modelling of crustal-scale magmatism

We employed a two-dimensional transient thermal model to quantify the crustal and mantle-derived melt evolution in the lower and upper crustal magma bodies and their influence on the melting of the surrounding crustal material over time (Karakas and Dufek 2015; Karakas et al. 2017). We assumed that the crust had a thickness of about 30 km, subdivided into a 23 km-thick, upper-mid crust consisting of amphibolite-facies metapelitic rocks, and a 7 km-thick, lower crust consisting of granulite-facies metapelitic rocks, in agreement with field observations (e.g., Quick et al. 2009; Brack et al. 2010). Initially, the crust is assumed to have a steady-state geotherm with assuming 0 °C at the surface and constant heat flux on the mantle-crust boundary corresponding to 70 mW/m2 heat flux on the surface, an average continental heat flux value shown in global compilations (Davies 2013). We assumed that radioactive elements are distributed in the upper 12 km of the crust, with exponential decay in concentration until this depth, which is consistent with geochemical data on the mid- to lower crustal rocks (Schnetger 1994). The heating as a result of radioactive element decay is calculated using the averaged values of U, Th, and K (Turcotte and Schubert 2002).

We assumed that mafic magmas start to intrude into the lower crust at the mantle-crust boundary 4 Myr before the caldera-forming eruption, as constrained by zircon U–Pb ID-TIMS ages of the lower crustal samples (Karakas et al. 2019). The magmas are intruded as a series of dykes and sills, where movement of the crustal material to both sides accommodates intrusion of dykes by the half-width of the dyke, while the emplacement of the sills is compensated by the downward movement of the crustal material according to the thickness of the sill. The model is divided in two stages of 2 Myr each. In the first stage, we focus on the evolution of the magmas in the lower crust that interact with the crustal material increasing the temperatures. We allow this lower crustal system to grow for 2 Myr (0–2 Myr). In the second stage, we only focus on the upper crust, which was pre-heated by the lower crustal intrusions and is characterized by the emplacement of magmas of intermediate-composition (extracted from the lower crustal system) at depths between 12 and 2 km (median 7 km depth). This stage lasts for an additional 2 Myr (2–4 Myr) (Karakas et al. 2019). In both lower and upper crustal intrusion scenarios, we assume that the frequency, volume, and timing of individual intrusion events are stochastic. We constrain the maximum and minimum limits of frequency and volume of the intrusions such that the total 2D volume intruded in the system matches the volume of igneous bodies within the total time frame constrained in the system (4 Myr from U–Pb dates, Karakas et al. 2019). The model is set to focus on large-scale (km) and long timescale (million years) evolution of the system, and does not focus on the details of the melt extraction processes that occur in much smaller (cm) and shorter (< year) scales. At the end of the simulations, we reach lower and upper crustal intruded volumes comparable to field observations of a 30 km × 7 km of lower crustal body overlain by a 5 km × 5 km upper crustal system (e.g., Brack et al. 2010).

The heat transfer from the intrusions to the crustal lithologies is mainly dictated by conduction (Biot number <<1, Carrigan 1988), where the effect of convection in the crust is negligible over the time (a few million years) and length scales (>10 km) we consider in the system. We solved the transient heat conduction equation with the latent heat as the source term. A predictor–corrector algorithm and underrelaxation was employed for numerical stability (Voller and Swaminathan 1991). The stochastic behaviour of the volume and timing of the intrusions were determined by the random number generator algorithm (L’ecuyer 1988). The details of numerical assumptions are included in Karakas and Dufek (2015) and Karakas et al. (2017).

The major difference between the current numerical model and the previous investigations (Karakas et al. 2019) is the quantification of the thermal and compositional evolution of lower and upper crustal systems at different spatial scales. For this purpose, we have defined in the model several, spatially restricted sub-zones of the magma bodies, referred to as “boxes” (Fig. 4A, B). Three boxes are defined in the upper crustal region (A–C) and five boxes are defined in the lower crustal region (I–V). The location of the boxes was selected such that we can observe the evolution of magmas in the center, on the edges, and in the middle parts of the magma bodies. In our model, the magma intrudes at different parts of the system over time and the material moves, while interacting with the crust around it. Hence, with defining these boxes, we are able to assess the smaller-scale changes in certain parts of the system and better link thereby the dynamics of the system as a whole with the isotopic record at the hand sample scale.

The melt fraction-temperature relationships and volumes of crustal lithologies and intruded magmas are constrained using the estimated parental compositions, pressures, and initial water concentrations based on the field observations. For upper/middle crustal lithologies (0–23 km), we have used experimental results of Vielzeuf and Holloway (1988) on metapelite lithologies, representative of the most common protolith of the Kinzigite Formation. For the lowermost 7 km of the crust (23–30 km), we used the melt fraction-temperature relationship suggested by Sinigoi et al. (2011), i.e., using a restitic metapelite lithology. We note that the crustal composition that we use from Vielzeuf and Holloway (1988) represents a fertile end-member for the crustal lithology (Kinzigite Formation). There is evidence that the Kinzigite Formation already underwent high-grade metamorphism and melt depletion by the time of the intrusion of the Mafic Complex (Ewing et al. 2014; Barboza et al. 1999; Barboza and Bergantz 2000). On the other hand, the lower crustal, granulite-facies rocks contain zircons crystallized from the cooling of an anatectic melt and whose U–Pb ages (ca. 290–250 Ma) overlap with, or even post-date the intrusion of the Mafic Complex (Guergouz et al. 2018; Kunz et al. 2018). Mass balance constraints on the Zr budget in some samples suggest that these zircons crystallized when the system still contained up to 25–45 wt.% melt (Mintrone et al. 2020). It is therefore likely that parts of the lower crust were not entirely melt-depleted when the Mafic Complex intruded. From the above discussion, we decided to use a fertile metapelite lithology as the upper- to mid-crustal composition (up to a depth of 23 km) for sake of simplicity and being consistent with the field observations.

For lower crustal intrusions, we used rhyolite-MELTS simulation results from a primitive, mafic starting composition, found as an enclave from the roof of the Mafic Complex (Sinigoi et al. 1994). The intermediate composition for the upper crustal intrusions is determined from rhyolite-MELTS simulations using as starting composition the granodiorite sample IVZ616-11 of Karakas et al. (2019). The melt fraction-temperature relationships of all lithologies used in the model are explained in detail in Karakas et al. (2019).

We quantify the melt fraction distribution for crustal and residual magmas, as well as the crustal melt contribution relative to the total amount of melt (crustal + residual from basalt fractionation), over the lifetime of the system. We calculate melt fractions at each individual pixel (40 m × 40 m) and then investigate the entire lower and upper crustal systems, while also gathering spatial and temporal information by quantifying the melt volumes in the upper and lower crustal “boxes”. The crustal melt fraction refers to the melted amount of crustal material due to heat conduction (vol. crustal melt / vol. crustal material in the pixel), while residual melt fraction is defined as vol. residual melt / vol. intruded melt in the pixel. More information on the assumptions for the density, heat capacity, and heat conductivity of different lithologies are included in Karakas et al. (2017).

We investigate snapshots of melt evolution in the system at regular time intervals, and quantify the instantaneous information from the system. The crustal melt contribution is defined as:

$${\text{crustal melt contribution}}_{{{\text{instantaneous}}}}^{{{\text{box}}}} \left( \% \right) = \frac{{{\text{volume}}_{{\text{crustal melt}}}^{t} }}{{{\text{volume}}_{{\text{total melt}}}^{t} }}*100$$
(1)

i.e., the percentage of crustal melt over total melt in the specified box at the specific time instant [t].

In our upper crustal simulations, we calculate the additional crustal melting that can occur in the upper crust, which adds to the contribution of lower crustal melting. To compromise between geological representativeness and computational efficiency, we took this into account by adding the space- and time-integrated average crustal melt contribution (%) from the lower crustal calculations to that calculated in the upper crust. Hence the corrected crustal melt contribution (%) in the upper crust is given by:

$${\text{crustal melt contribution}}_{{{\text{corrected}}}}^{{\text{box, upper crust}}} = {\text{crustal melt contribution}}_{{t\left( {{\text{end}}} \right)}}^{{\text{box, upper crust}}} + {\text{average crustal melting}}_{{t\left( {{\text{end}}} \right)}}^{{\text{overall system, lower crust}}}$$
(2)

where, \({\text{average crustal melting}}_{t(\mathrm{end})}^{\text{overall system, lower crust}}\) refers to the average, spatially integrated crustal melt contribution (%) in the lower crustal calculations at t(end) = 2 Myr.

To compare numerical simulations with the Monte Carlo mixing models, we have added an algorithm to our thermal model to track the history of each individual crustal melt volume rather than the instantaneous information gathered following the method described above. To retrieve this information, we mark all the crustal lithologies (each pixel) at the beginning of the simulation and then track the temperature and melt fraction history of each pixel. The ratio of maximum melt fraction value for each pixel by the total melt fraction is used as a proxy for the crustal melting productivity (%) over the lifetime of the system.

Inverting zircon Hf isotopic compositions into crustal melt/total melt ratios

To provide a direct comparison proxy with results of the thermal modelling, we inverted the Hf isotopic composition of each analyzed zircon into a crustal melt/total melt ratio using a two-component mixture model between mantle and crustal components. The mixture model uses a Monte Carlo Markov Chain routine with 103 individual models randomly varying the end-member Hf concentrations and isotopic compositions (Storck et al. 2020). The results of such models are very sensitive to the choice of the mantle and crustal end-member values, which we justify in the following.

As discussed in Schaltegger and Brack (2007), the nature of the mantle source of Permian magmatic rocks is uncertain, but could range between two end-members: an asthenospheric, isotopically depleted mantle (DM), similar in composition to the source of MORB; and an isotopically enriched lithospheric mantle. Accordingly, we have set the range of Hf isotopic compositions for the “primitive” mantle melt between εHf = +4, corresponding to the most primitive zircon εHf measured in this study; and εHf = +12, corresponding to the composition of the New Crust model of Dhuime et al. (2011) at 283 Ma. The lower bound represents a conservative estimate for enriched lithospheric mantle, as previous estimates proposed a minimum εNd value of +5 (Stille et al. 1996; Pinarelli et al. 2002) corresponding to εHf = +9 based on the global Hf–Nd array of Vervoort et al. (1999). The upper bound is consistent with the isotopically most depleted mantle-derived material of the Permian mafic complexes, i.e., the (260 Ma-old) melanocratic dykes of the Mont Collon area (western Alps) with εNd of +7 (Monjoie 2004) corresponding to εHf of +12. Concentrations were set as to vary between 1 and 6 ppm Hf, obtained from the representative range of Zr contents in samples of the Mafic Complex excluding outliers (30–200 ppm; Rivalenti et al. 1975; Sinigoi et al. 1994; 2011) and a canonical Zr/Hf ratio of 35 (Rudnick and Gao 2014).

In crustal rocks, Hf is mainly hosted in detrital zircons, whose isotopic variability is huge in this context. Indeed, the inherited zircon grains measured in this study show a range of εHf calculated at 283 Ma of 0 to −38 (Fig. 3, lower panel) and down to −55. Ewing et al. (2014) obtained Hf isotopic compositions of detrital zircon from metapelites of the Kinzigite Formation and reported similar ranges of εHf (0 to −70). However, the data stem from individual detrital/inherited cores, which would have been dissolved in melt in the high-temperature conditions of the lower crustal mush (> 900 °C maintained for several Myr; see modelling here and in Karakas et al. 2019). The Hf isotopic composition of such melt would homogenize to converge to the average εHf of the assimilated zircons. Therefore, we have used as bounds the range of average εHf calculated from detrital zircon populations in several samples of the Kinzigite Formation, i.e., −6 to −20 (Ewing et al. 2014). These values are consistent with the distribution of εHf observed in our inherited zircons, whereby the main population ranges in εHf between −3 and −15 with a few grains with much lower εHf (−35 to −55). The range also compares well with average εHf of detrital zircon from meta-sediments elsewhere in western Europe calculated at 283 Ma, i.e., −9 to −22 (from data of Chelle-Michou et al. 2017; Couzinié et al. 2019). The Hf concentrations of the crustal component was set as to vary between 4 and 21 ppm, which is the range of bulk-rock Hf concentrations from the samples of the Kinzigite Formation studied by Ewing et al. (2014). We acknowledge that the crust would not be assimilated as bulk metapelite, but this rather attempts at modelling conservatively the wide range of Hf contents observed in the local lower crust.

Fig. 3
figure 3

Crustal-scale sketch of the Ivrea-Sesia magmatic system with expected sample locations in the reconstructed stratigraphy. Distribution histograms of single grain εHf283 Ma values of all measured zircons are shown for each sample as a function of relative stratigraphic position. Individual sample locations are reported in Fig. 1. Kernel density distributions of εHf measurements are shown for each sample with 2σ error. Number of measurements are included with the parameter n. Bottom panel shows the distribution of εHf283 Ma values in 24 inherited zircon grains (LA-ICP-MS age >300 Ma)

Full size image

Results

Zircon εHf values of the SMS

In total, 226 Hf isotopic analyses on 183 zircon grains were obtained from 12 igneous rock samples (Fig. 1 and Supplementary Table S1). We do not observe any correlation between εHf and the 176Yb/177Hf ratio for reference materials and samples (Fig. 2). This suggests that isobaric interferences of 176Yb are accurately corrected and thus that the observed variations in εHf have a geological significance. Figure 3 shows the distribution of the obtained zircon εHf values across the entire crustal-scale igneous section. Positive zircon εHf values occur only in the ultramafic pyroxenite unit, ranging from +0.5 ± 1.1 to +3.6 ± 1.3 (all uncertainties are quoted at 2σ). All of the other samples show negative zircon εHf values, from the lower crustal mafic units (gabbros and diorites) to the upper crustal silicic units (granites and rhyolites), and their ranges overlap to a large extent. Lower crustal gabbros and “diorites” show a very similar range of zircon εHf between −4.2 ± 1.0 and +0.2 ± 1.0, with all samples being characterized by a nearly normal distribution centered at ca. −2.5 (Fig. 3). Compared to lower crustal rocks, upper crustal silicic rocks cover a wider range in zircon εHf values (−9.5 ± 1.3 to −0.2 ± 1.1), present in all individual samples, and thus display a slight shift to more negative values on average (median value at ca. − 4.1). The Roccapietra dyke also shows a closely overlapping range of zircon εHf with those of the upper crustal units (−4.9 ± 0.9 to −0.8 ± 1.4). The zircon εHf distribution of the rhyolites shows a sharp, single peak at ca. −4 (−4.2 ± 0.7 to −3.9 ± 1.1), which differs from the observations from other silicic units but is likely due to the small number of zircons suitable for Hf isotope analysis in this sample (n = 4). The inherited zircons (U–Pb age >300 Ma) spread over a wide range of εHf from 0.7 ± 1.1 down to −55.5 ± 0.8, with a main population characterized by εHf between −5 and −12, consistent with the data reported by Ewing et al. (2014) (Fig. 3, lower panel—see above and Supplementary Table S1).

Thermal modelling results

The results of the thermal model are presented in the form of two-dimensional plots of temperature and melt fraction (Fig. 4A, B) as well as the space–time evolution of the crustal contribution (%) in lower and upper crustal systems (Fig. 4C, D). The two-dimensional temperature and melt fraction profiles at 2 Myr show that the temperatures towards the center of the lower crustal magma body increased over 900 °C range, together with the melt fractions, which go up to 0.4–0.5 in the lower crust. However, volumetrically, the dominant melt fraction is < 0.3 and mostly composed of residual, fractionated mantle-derived melt. In the upper crust, the magmatic system is much smaller (Fig. 4A) and reaches temperatures of only 700–750°C. Accordingly, the melt fraction range in the upper crustal magma body is typically lower than the lower crustal system (although some locally increased melt fraction near a recent intrusion can happen from time to time).

Fig. 4
figure 4

Two-dimensional temperature and melt distribution in the crustal-scale magmatic system, calculated from thermal model. A Temperature profile of the lower and upper crust after 2 Myr of magma emplacement in the lower crust (0–2 Myr), followed by additional 2 Myr of magma emplacement in the upper crust (2–4 Myr). The boxes I–V (lower crust) and A–C (upper crust) show the corresponding areas of melt fraction calculations in C and D. B Two-dimensional distribution of melt fraction in the lower and upper crust, showing crustal and residual contributions to the total magma generated in the system as a result of dyke and sill intrusions that changed the crustal thermal structure as shown in panel A. C, D Temporal and spatial variations of the crustal contribution (%) (as defined by Eq. 1) and melt fractions in the lower (C) and upper (D) crustal magma bodies, at the scale of the box locations determined in A (upper panels, “Boxes”) and the entire lower/upper crust (lower panels, “Systems”). In upper panels, instantaneous crustal contribution (%) refers to the percentage of crustal melt in the specified box at each specific time step. Lower panels display both (1) the instantaneous crustal vs. residual melt fractions (colored histograms); and (2) the crustal melt contribution (%) normalized to the final amount of melt in the system after the 2 Myr duration of the model. In upper crustal calculations in D, crustal melt contributions are plotted both considering only upper crustal assimilation (left) and adding the average crustal melt contribution (%) of the lower crust (20%) (see Eq. 2)

Full size image

We further assessed the magma evolution in five areas in the lower crust and three areas in the upper crust, distributed in different parts of the magma body (Fig. 4A, marked boxes). Although the system likely experienced thermal convection to some extent, we note that the boxes, and hand specimens that represent an even smaller scale, are not representative of the entire averaged system, but they are the products of processes that have occurred in different portions of the magmatic system at different times. We display our observations in the lower and upper crustal sections separately, showing the evolution of melt fraction and/or crustal melt contribution spatially averaged over the entire lower and upper crustal “systems” (Fig. 4C, D upper panels) and at the scale of the “boxes” (Fig. 4C, D upper panels) defined in Fig. 4A.

In the lower crust, snapshots of melt fractions at different successive time steps evaluate the temporal evolution of the system. The overall system (Fig. 4C upper panel) remains at low melt fractions (<0.35 for mantle-derived, <0.2 for crustal melts) over the 2 Myr duration of the model, which creates a thermal buffer and keeps the entire system in a mush state over this long timescale. Overall, the instantaneous crustal melt contribution is between 20 and 40% within most of the lifetime of the lower crustal system and reaches about 20% at 2 Myr (red dashed- and solid lines, Fig. 4C upper panel). However, the crustal melt contribution differs significantly both temporally and spatially, from 0 to up to 80%. The systematic increase of melt fractions and temperatures towards the center of the system (Fig. 4B, C, boxes III and IV) is due to the incremental intrusion of primitive magmas from the central portion of the system over time. The increase of crustal melt contribution (%) towards the top of the magmatic system (Fig. 4B, C, box V) is explained by the change of the crustal lithology in our model from refractory granulite-facies metapelitic rocks to more fertile amphibolite-facies metapelites at the depth of 23 km (Fig. 4B).

In the upper crust, the model generates volumetrically a much smaller system compared to the lower crust, as constrained by field observations. We observe that the instantaneous melt fractions of the entire upper crustal system is <0.2 for the mantle-derived contribution, and <0.05 for the crustal lithology (Fig. 4D upper panel). Our results show that the crustal melt contribution in the upper crust is lower than that of the lower crustal system (<10% on average in the upper crust, while the average is ca. 20% for the lower crust; Fig. 4C, D upper panels). An early time snapshot displays local, transient crustal melt contributions up to >90% in the overall system (yellow dashed-line). However, we stress that the total melt volumes at these early times are, on average, very low (melt fractions for both crustal melt and residual mantle melt are both <0.01) and such high crustal melt contributions are not representative of the overall evolution of the magmatic system (Fig. 4D, upper panel). Although not as pronounced as for the lower crust, there are also local, instantaneous variations of the crustal melt contributions between boxes I, II and III, mainly between 0 and 20% but showing transient values as high as 50–70% (Fig. 4D lower panel).

Monte Carlo mixing simulations and comparison with numerical model

Monte Carlo mixing simulations based on Hf isotopic data

The mixing model with Hf isotopic and concentration ranges ratios (mantle:crustal contribution, Supplementary Table S2) shows a large range of possible crustal melt contributions (%), ranging from 6 to 27% in gabbros, 13–29% in diorites, and 15–84% in upper crustal silicic rocks (min and max of median values calculated for each measurement, see Fig. 5 and Supplementary Table S2). Kernel density distributions of lower crustal samples (Fig. 5A) follow overlapping distributions, characterized by peak in the range of 10–30% crustal melt contribution for the gabbro and diorite samples. The pyroxenite sample shows a distinct distribution skewed towards the mantle end-member, with crustal melt contributions of 5–15%. Comparison of these calculations with the pixel-by-pixel tracer calculations in the numerical model (see “Methods”) shows how the field and statistical observations (kernel density plot) are reproduced by the numerical estimations (bar plots). In the lower crust (Fig. 5A histograms), the space- and time-integrated, average crustal melt contribution shows 20% (purple bar). The crustal contribution decreases from 7% at the lower edge of the system (Box I) to 4% at the center (Box III), and increases to 13% at the top of the system (Box V). The Monte Carlo simulations for the assumed Hf concentration ratios show a kernel distributions peak at around 18% crustal melt contribution for the lower crustal samples, which compares well with the averaged (20%) crustal melt contribution of the thermal model for the lower crust.

Fig. 5
figure 5

Results of Monte Carlo simulations of mixing between crust-derived and residual melts gathered with thermal model results (histograms) for upper (A) and lower (B) crust. The kernel density distribution for each sample is plotted with lines of different colours. Histogram values show the average crustal melting using the tracer model (see “Methods”) for each individual box (histograms at the bottom panel) and the overall system (histogram at the top panel, overlain with the kernel distributions). The two crustal melting values in the upper crustal calculations shown with histograms correspond to the sole crustal assimilation in the upper crust between 2 and 4 Myr (left), and corrected crustal assimilation by summing up the lower and upper crustal assimilation (right) considering that the upper crustal intruded magmas initially contain 20% crustal contribution (%) on average from assimilation of lower crustal wall rocks between 0 and 2 Myr (right) (see text for details)

Full size image

In the upper crust (Fig. 5B), overlapping kernel density distributions of upper crustal samples reveal close similarity of all upper crustal units; the crustal contribution mainly ranges from 10 to 50%, with a peak around 20%. In the thermal calculations, we do not observe any significant difference between boxes A, B, C and the overall system. The numerical results (bar plot) show that the additional crustal contribution from upper crustal rocks is very low, around 1–2% for the long-term averaged system (compared to instantaneous calculations in Fig. 4D). However, it is important to consider that the upper crustal magmas were extracted from the lower crustal system where crustal contamination happened to a greater extent, as the average, space- and time-integrated crustal melt contribution in the lower crust is 20% (Fig. 5A purple bar). We therefore considered this by adding this average crustal contribution from the lower crustal system (20%) to that of the upper crust (1–2%) (see Eq. 2 above) to reach a total crustal contribution of 21–22% in the granitoids. Again, this reasonably fits the median values of the Monte Carlo calculations based on Hf isotopes.

Monte Carlo mixing simulations based on Sr and Nd isotopic data

We have implemented the Sr, Nd and bulk rock isotope geochemical dataset of Voshage et al. (1990) from lower crustal rocks into a Monte Carlo mixing model similar to that conducted based on our Hf isotopic record (Fig. 5C; Supplementary Tables S3, S4). Kernel density functions of resulting mixing calculations follow an overall trend similar to the findings based on the Hf isotope data. Pyroxenites show a bimodal distribution (ca. 8–45% for Sr; ca. 10–50% for Nd) with a main peak representing the lowermost crustal contribution, followed by gabbros and “diorites”, with overlapping crustal melt contribution (ca. 15–50% for Sr; ca. 20–50% for Nd). The concentration ratios for Sr and Nd by Voshage et al. (1990) show a large spread of a potential primary mantle-derived melt mixing with a crustal contaminant (Fig. 5C).

Discussion

Isotopic architecture of a crustal-scale magmatic system

Hafnium

When compared to the extreme εHf variability observed between the local crust (εHf of inherited zircons between 0 and –55 and average values between –6 and –20; Fig. 3 and Ewing et al. 2014) and putative mantle components (εHf from ca. +4 to +2; see “Methods” section), zircon Hf isotopic compositions across different igneous bodies of the SMS show a relatively tight range (Fig. 3). Specifically, 90% of the dataset shows εHf between −0.5 and −6.7 and 50% between −2.0 and −4.0 (with typical 2σ uncertainties of ± 1.1 εHf-unit). Ranges of εHf observed in all individual samples overlap with each other, suggesting that the crustal scale magmatic system is not only temporally (Karakas et al. 2019), but also genetically linked. Our results based on Hf isotopes therefore confirm the vertical stratigraphic link between all igneous units in the system, in agreement with the observations and interpretations of Voshage et al. (1990) based on Nd isotopes, and Sinigoi et al. (2016) based on Sr isotopes and trace elements. Together with previous work (Voshage et al. 1990; Peressini et al. 2007; Quick et al. 2009; Sinigoi et al. 2016; Karakas et al. 2019), our current εHf measurements suggest a polybaric differentiation column that evolved over a few million-year timescales, where the first mafic magma intrusions were emplaced in the lower crust, forming the incrementally growing lower crustal mush system. Over time, the temperature of the system increased, and magma started to be emplaced in the upper crust forming silicic bodies and caldera-forming eruptions at the peak of the magmatic activity.

While crustal-scale observations confirm the coeval and co-genetic nature of all igneous units in the entire SMS, the heterogeneity of mantle and crustal end-members and the complexity of their interactions inevitably leads to the observation of local heterogeneities between individual samples. The range of Hf isotopic compositions in zircon grains of a given sample indeed record the complex history of magma formation and evolution over time and space. Depending on their location in the evolving magmatic body, different samples are therefore expected to show deviating ranges in zircon Hf isotopic composition, which can in turn allow to place constraints on the history of the system as a whole.

In the lower crust, the only outlier in our Hf isotopic dataset is the pyroxenite body located in the southeastern portion of the Mafic Complex (near Castello di Gavala). This kind of pyroxenites have been studied by Garuti et al. (2001) and interpreted as ultramafic pipes (εNd: +1.6 to +3.69). Voshage et al. (1990) measured highest Nd isotopic values in cumulus peridotites (εNd: +7) and pyroxenites whereas all other units showed negative values. Likewise, Sinigoi et al. (1991, 1994, 2011) report very homogeneous Sr and O isotopic compositions across most of the Mafic Complex stratigraphy, but shows more scatter (notably values trending towards “mantle-like” compositions) in the lowermost part (richer in ultramafic cumulates).

In post-collisional settings an isotopically heterogeneous mantle can not be excluded (e.g., Soder and Romer 2018). Considering all available data, these isotopic differences between pyroxenite vs. gabbro/diorite becomes more complex and suggests several possibilities: (1) pyroxenites and diorites/gabbros may be originated from different batches of magma derived from isotopically different mantle sources (respectively more depleted and more enriched); (2) if the magma source is relatively homogeneous, the pyroxenites recorded a more primitive isotopic composition while the gabbro/diorites are more strongly, and uniformly, affected by crustal assimilation/contamination (Voshage et al. 1990; Sinigoi et al. 1994, 2011, 2016). In addition to the coherent set of Hf isotopes in zircon presented here, we argue that the second possibility is more likely because of (1) the close overlap of U–Pb ages of the pyroxenite body compared to other lower crustal units (Karakas et al. 2019); and (2) the occurrence of distinctively more radiogenic isotopic composition of possible Permian mantle components in the area (εHf from ca. +9 to +12 inferred from Nd isotopic data of Stille and Schaltegger 1996; Pinarelli et al. 2002 and Monjoie 2004) compared to the pyroxenite zircons (εHf of 0 to +4). Voshage et al. (1990) further argue that dramatically different mantle sources as well as complex mixing between them would be necessary to explain their observations, i.e., “stratigraphic variations” in both Sr and Nd isotopic compositions. Altogether, contamination of magmas derived from a common mantle source rather than multiple (or isotopically heterogeneous) mantle sources is therefore considered to be a likely scenario for the SMS.

Based on our thermal model, the crustal melt contribution in the lower crustal magmatic system can vary temporally from 0 to 40% (Fig. 4C upper panel, red dashed line), and spatially between 0 and 80% (Fig. 4C lower panel). It is therefore possible, that pyroxenites intruded either at a time when, or in a location where thermal conditions were not favourable to crustal melting, such that they crystallized rapidly without assimilating significant amount of crustal material.

In the upper crust, zircons from the granites have more diverse, and more negative εHf values (average −3.9 ± 3.9, 2 S.D.) than in the lower crustal rocks (average −2.5 ± 1.5, 2 S.D. excluding the pyroxenite sample). In addition, εHf values of granites are on average higher than those of most inherited zircons (for which 85% of the values are below −6; Fig. 3). The shift towards more negative values relative to the lower crustal zircons may be explained by additional assimilation of crustal melts during magma emplacement in the upper crust. However, this hypothesis is precluded by both field evidence, as anatectic, two-mica granites in the main upper crustal granites are restricted to volumetrically very minor dykes and pods intruded late in the igneous story (Sinigoi et al. 2016); and the present thermal model (Fig. 4) strongly suggesting that melting the upper crust during granite emplacement is rather inefficient, representing only a few % of the upper crustal magmas. We cannot exclude that larger extents of assimilation took place locally in the upper crust and/or the newly crystallized zircons inherited variable εHf from dissolved xenocrysts due to inefficient homogenization of Hf isotopic compositions in the melt, as described for S-type granites (Villaros et al. 2012). However, we emphasize that these processes would only have local effects and cannot explain the general trend observed at the scale of the entire section.

Alternatively and speculatively, the shift towards lower and more diverse εHf in the granites without significant upper crustal contamination can be explained by the contrasted thermal conditions and physical mechanisms of hybridization between the lower and upper crustal systems. The lower crustal system is relatively mafic in composition and kept at temperatures hotter than required for zircon saturation for most of its lifetime (Karakas et al. 2019), so that significant zircon saturation would only occur when the mush has sufficiently cooled down and the melts chemically evolved towards more silicic compositions, which happened within the last 50–150 °C before the solidus (Karakas et al. 2019). In addition, the higher lower crustal temperatures and melt fractions promote efficient mixing and hybridization between residual melts fractionated from the Mafic Complex and crustal melts delivered from rocks progressively incorporated in the growing Mafic Complex (Sinigoi et al. 2016), explaining the tight range of zircon εHf values. In contrast, extraction of melts from the lower crust and their emplacement in the colder upper crust (mostly < 850 °C throughout the lifetime of the entire system; Karakas et al. 2019) will promote rapid zircon saturation after emplacement. The extracted melts may have sampled different instantaneous and/or local proportions of melts derived from basalt fractionation vs. crustal melting, which vary temporally from 0 to 40%, and spatially between 0 and 80% (locally), in the lower crustal system (Fig. 4C). The extracted melts were potentially poorly mixed and their isotopic composition not necessarily representative of the space- and time-integrated proportion at the scale of the entire lower crustal system. This would translate in a larger scatter of zircon εHf as observed in the natural samples. Based on these results, we argue that (1) protracted extraction of isotopically heterogeneous hybrid melts from the evolving lower crustal mush played an important role in producing the more evolved, upper crustal bodies; and (2) crustal melting alone is not the dominant mechanism to form these magmas as their zircon εHf values still largely overlap with those of the Mafic Complex.

Strontium and neodymium vs. hafnium

Monte Carlo mixing calculations of Sr and Nd isotopic data on lower crustal lithologies in the Ivrea system show moderate crustal contribution (20–50%) when the lowest and highest concentrations measured in the Voshage et al. (1990) dataset were used in the mixing model. Although the absolute shift in crust/mantle ratios of different isotope systems (e.g., Hf vs. Sr and Nd) depends on the endmembers chosen for the mixing model, the overall trend (increasing crustal contribution from pyroxenites to gabbros and diorites) remains similar for all of them.

A possible explanation to the wide range in modelled crustal contributions is that the Hf, Sr and Nd concentrations and isotopic signatures are controlled by different mineral phases (zircon for Hf; micas, feldspars and apatite for Sr; apatite, monazite and titanite for Nd). Due to different phase stabilities and implications in melting reactions, the relative timing of crystallization or mineral break down reactions likely varies throughout the interactions between mantle- and crust-derived magmas, resulting in some extent of isotopic decoupling. For instance, the main melting reactions in the pelitic systems involve the breakdown of micas with high Rb/Sr ratios, imparting very radiogenic and out-of-equilibrium Sr isotopic compositions to the melts (Farina and Stevens 2011). Likewise, although Nd theoretically shows a behaviour similar to Hf (Vervoort et al. 1999), the different solubility of monazite and zircon in melt at a given temperature (Kelsey et al. 2008) would readily decouple Nd and Hf concentrations and isotopic compositions in crustal melts. Additionally, we note that Sr and Nd isotopes are from bulk-rock measurements in the same samples (Voshage et al. 1990), while Hf isotopes were conducted on single zircon grains in different samples (current study). This sampling issue may further increase the apparent mismatch between different isotopic systems, considering the large spatial and temporal variations of crust vs. mantle proportions as observed in the thermal model (Fig. 4).

Considering the abovementioned sources of isotopic decoupling, the discrepancy between the different models is minimal, and we conclude that the results obtained with different isotopes likely overlap within the resulting uncertainties. In addition, despite having different concentration ratios, Nd and Hf agree well in their isotopic signature and eventually probe the same magmatic processes precluding the need of two isotopically distinct mantle sources. Furthermore, models based on all isotopic systems tend to show a subordinate but significant proportion of crustal contribution relative to the fractionation of mantle-derived basalt, corroborating our conclusions based on Hf only.

Spatial and temporal heterogeneity of magmatic record

Our thermal results show that, although crustal melting can show locally significant values in the lower crust, mantle-derived magmas are volumetrically dominant in the overall crustal column (Figs. 4, 5). In the lower crust, melting of pre-existing lithologies typically reaches a maximum at the top of the growing magma body. However, the amount of melting in the crustal rocks strongly varies in time and space; parts of the systems can reach >55% melting locally (Fig. 4C, box III early timescales; box V late timescales) while others may remain mostly below their solidi (Fig. 4C, box I; Fig. 4D box A, B).

Previous thermal studies have quantified the crustal melting efficiency (or amount of assimilation) using different scenarios and assumptions in the lower and upper crust (e.g., Annen and Sparks 2002; Dufek and Bergantz 2005; Barboza and Bergantz 1996, 1997). Other studies used geochemical trends together with mixing models to quantify the amount of crustal melting and magma evolution in a wide variety of settings (e.g., Taylor 1980; DePaolo 1981). All these models have argued a wide range of possible degree of assimilation in their systems, either low (< 5%) or high (> 60%), suggesting very different implications considering the heat and mass balance of the system. Crustal lithology and volatile budget (fertility) have been shown to affect the resulting melting process and therefore dictate the magnitude of crustal melt/contamination of the magmatic system (e.g., Storck et al. 2020).

In our model, we used a refined systematic quantification, where we quantified possible amount of melting of pre-existing crustal rocks from both the lower crust and upper crust, while calculating the spatial and temporal changes over the lifetime of the system. We suggest that the amount of assimilation in the same system vary spatially and temporally by more than an order of magnitude (Fig. 4C), depending on the local conditions. However, we stress that when considering the entire system, the volume of crustal magma is subordinate, but necessary to form the upper crustal granites of the SMS (see also Sinigoi et al. 2016). Based on both modelling and isotopic constraints, the SMS column has a contribution from the crust likely ranging from ~20% (average value for the thermal model, Figs. 3, 4, 5) to ~40% (based on some of the highest values provided by the isotopic mixing models tailing to 70–80%). The former value is probably more representative, because it integrates variations of the crustal contributions in space and time over the lifetime of the system. In contrast, as discussed previously, isotopic data may sample spatially or temporally discrete variations of the crust/mantle ratio, can be already contaminated in mantle-derived material (crustal contribution from the lithospheric mantle is not considered in the thermal model) and depend on the behaviour of the mineral phase hosting the isotope of interest (i.e., zircon in the case of Hf). We also show that most of the assimilation occurs in the lower crust (with only negligible amounts occurring in the upper crust).

Finally, we stress out that the Kinzigite Formation was at least partly melt-depleted by the time the Mafic Complex has intruded (Barboza et al. 1999; Barboza and Bergantz 2000). Therefore, the crustal melting proportion estimated from our model may be regarded as a maximum, as modelling of phase relations in the crust was performed in our model using a fertile metapelite composition in the mid-upper crust (Vielzeuf and Holloway 1988) and more restitic composition (from Sinigoi et al. 2011) only for the lowermost 7 km of the crust.

Constraints and limitations of thermal and geochemical models

Overall, the crustal contribution estimated by the bulk rock Sr and Nd isotopic data in the lower crust agree well between each other. This dataset also shows a similar range with the estimations from the zircon Hf isotope measurements, the thermal model, and the Hf mixing model using our range of Hf concentration ratios. As discussed above, we conclude that this discrepancy between crustal contribution estimations are not dramatic considering the limitations of both geochemical and thermal models, complexity of the system, and resulting assumptions and simplifications. Following the detailed discussions in the preceding sections, we summarize these limitations and assumptions below.

The limitations in the geochemical models include (1) selection of the end-members in concentration ratios and isotopic compositions, which depend on crust and mantle source compositions, melting reactions and their spatial–temporal evolution, (2) bulk rock analysis of Sr and Nd isotopic measurements versus single-grain analysis of Hf measurements, (3) the use of the same samples for Sr and Nd analyses while Hf analyses were done on a different set of samples, and (4) possible sampling biases in and between different studies. The ranges of possible Hf concentrations and isotopic compositions in the geochemical model presented here were selected to be very conservative on purpose, in an attempt to mitigate these sources of uncertainty. Despite this, on the first order the results are consistent with modelling based on other systems. This gives us confidence on the robustness of our main interpretations, i.e., that the contribution of the crust is higher in the upper vs. lower crust of the SMS, yet remains subordinate to that of mantle-derived magmas at the scale of the entire system.

On the other hand, we acknowledge that the thermal model includes several simplifications and assumptions. For the upper- to mid-crustal lithology, we used an unmetamorphosed, fertile end-member composition parameterized from Vielzeuf and Holloway (1988) experiments. Fieldwork and previous thermal model studies showed that the Kinzigite Formation went through high-grade metamorphism, partial melting and to some extent large-scale melt depletion before the intrusion of the Mafic Complex (Barboza and Bergantz 1996, 2000; Barboza et al. 1999). Therefore, our thermal modelling results probably place a maximum bound for the amount of crustal contribution in the SMS. We note that we did not account for water transfer from crystallizing basalt to the crust, which would result in higher amounts of crustal melting (Collins et al. 2016). However, the primary magmas in extensional settings such as the SMS may be slightly drier (<2 wt.% H2O) than in typical arc systems (e.g., Plank et al. 2013) where the magmatic water contents are ~4 wt.% on average. We also note that in our crustal assimilation calculations, we consider all the crustal material that can be melted (anything with some fraction of melt in it) while in reality extracting melt at low melt fraction is extremely difficult, and would not allow those melts to participate in constructing upper crustal reservoirs and volcanic units. Also, the thermal model does not consider the possible efficiency of reactive bulk rock assimilation (Beard et al. 2005), which possibly acted in the Mafic Complex (Sinigoi et al. 2016).

Nevertheless, the combination of thermal modelling with (isotope-)geochemical datasets as done in this study, (1) delivers new insights on the magmatic evolution of the SMS as a whole, (2) places constraints (outcrop scale) spatial heterogeneities that form and vanish through space as a function of time.

Conclusions

Igneous rocks of the SMS show a remarkable coherence of zircon Hf isotopic compositions when compared to local mantle and crustal components, consistent with the consanguineous relationships of the different units proposed based on single grain ID-TIMS age determinations (Karakas et al. 2019). This suggests a continuous evolution, involving differentiation coupled to crustal assimilation, in the entire crustal-scale igneous system. Although the evidence shows that the system is coeval, the εHf values span a range between the deep and upper crustal levels. Zircons from gabbroic rocks of the upper Mafic Complex show a tight εHf range (−4.2 ± 1.0 to +0.2 ± 1.0; average −2.5 ± 1.5, 2 S.D.) in all analyzed samples, except one pyroxenite containing zircons with positive εHf (0 to +4). Geochemical modelling indicates that this most probably corresponds to between 10 and 40% crustal assimilation from mantle-derived magmas. Upper crustal silicic rocks cover a wider range in zircon εHf values (−9.5 ± 1.3 to −0.2 ± 1.1), and thus display on average a slight but significant shift to more negative values. Considering that the ages of both lower and upper crustal igneous bodies are virtually the same (Karakas et al. 2019) and that gabbroic rocks are cumulates, which evidently lost a significant amount of residual melt, it appears that those zircons in granite that share the same isotopic composition of gabbro cumulates crystallized from residual melts delivered from the Mafic Complex. On the other hand, zircons shifted towards more “crustal” compositions crystallized from silicic melts with a more crustal signature. Such melts cannot derive from melting in upper crust, which is minor as suggested by both field evidence and thermal model, and likely were produced in the deep crust, in agreement with the thermal model. Therefore, the main granitoids of the SMS are hybrids of residual and anatectic melts (in agreement with Sinigoi et al. 2016).

In addition to the general large-scale differences in Hf isotopic composition between the lower and upper crust, we observe that there is significant scatter from a sample to another. We explain this scatter by the thermal and compositional differences over time and space within and around the magma column. The thermal model shows that the crustal contributions vary extensively in space and time during the evolution of the lower crustal mush system (from 0 to 80%), but the average, space- and time-integrated value is about 20% of the total mantle-derived magma input. While zircon Hf isotopic variations between hand samples likely reflect spatially and temporally discrete events of melt extraction/crystallization from this evolving lower crustal mush, the bulk crustal contribution obtained based on the thermal model (~20%) averages out spatial and temporal differences and should be fairly representative at the crustal scale. Melt extraction from different parts of the magma body as well as timing of extraction as shown in this study can therefore lead to different crustal melt contributions. As the SMS did not achieve crustal-scale homogenization throughout its lifetime, we emphasize that one can only quantify/model “average crustal contributions” by integrating all spatial scales up to the entire crust.

Comparison of the analytical measurements together with detailed numerical modelling of the thermal evolution of the igneous bodies gives a first-order understanding of the long-term, large-scale evolution of the system. We suggest that hot, mafic magmas incrementally emplaced in the lower crustal depths, where they interacted with deep crustal refractory lithologies and overlying slightly more fertile ones. Over 2 million years of lower crustal magmatism, the thermal conditions evolved such that long-term melt presence occurred, leading to magma differentiation. This differentiation was fostered by assimilation of silicic partial melt originating from the pre-existing crustal lithologies and by melt extraction and percolation from the highly crystalline lower crustal mush zone. Such process of formation of hybrid silicic magmas cannot work in every geotectonic environment: the SMS developed in a slowly extending, warm and potentially still fertile crust due the previous Variscan orogenic event, an ideal condition to facilitate the involvement of crustal components in igneous systems (Moyen 2020). Similar conditions are more difficult to be attained, for instance, in rift environments, where the crust is thinner and extension faster, facilitating the rise of mafic melts at higher crustal levels or at the surface, without significantly impacting the crust. Hence, understanding the evolution of the SMS may help better constraining how local tectonic conditions influence the proportions of crustal recycling by mantle-derived melts to produce evolved magmas.