# Mantle versus crustal contributions in crustal-scale magmatic systems (Sesia Magmatic System, northern Italy) from coupling Hf isotopes and numerical modelling

Julian-Christopher Storck, Oscar Laurent, Ozge Karakas, Jörn-Frederik Wotzlaw, Andrea Galli, Silvano Sinigoi, Olivier Bachmann, Cyril Chelle-Michou

► To cite this version:

Julian-Christopher Storck, Oscar Laurent, Ozge Karakas, Jörn-Frederik Wotzlaw, Andrea Galli, et al. Mantle versus crustal contributions in crustal-scale magmatic systems (Sesia Magmatic System, northern Italy) from coupling Hf isotopes and numerical modelling. Contributions to Mineralogy and Petrology, 2021, 176 (11), 10.1007/s00410-021-01847-9. http://hal-03436644

HAL Id: hal-03436644

https://hal.science/hal-03436644

Submitted on 9 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mantle versus crustal contributions in crustal-scale magmatic systems (Sesia Magmatic System, northern Italy) from coupling Hf isotopes and numerical modelling

Julian-Christopher Storck<fnref n="1" />,<fnref n="2" />✉ · Oscar Laurent<fnref n="1" />,<fnref n="3" /> · Ozge Karakas<fnref n="1" /> · Jörn-Frederik Wotzlaw<fnref n="1" /> · Andrea Galli<fnref n="1" /> · Silvano Sinigoi<fnref n="4" /> · Olivier Bachmann<fnref n="1" /> · Cyril Chelle-Michou<fnref n="1" />

Received: 9 October 2020 / Accepted: 21 September 2021 / Published online: 23 October 2021
© The Author(s) 2021

## Abstract

The growth and evolution of crustal-scale magmatic systems play a key role in the generation of the continental crust, the largest eruptions on Earth, and the formation of metal resources vital to our society. However, such systems are rarely exposed on the Earth's surface, limiting our knowledge about the magmatic processes occurring throughout the crust to indirect geochemical and petrographic data obtained from the shallowest part of the system. The Hf isotopic composition of accessory zircon is widely used to quantify crust-mantle evolution and mass transfers to and within the crust. Here we combine single-grain zircon Hf isotopic analysis by LA-MC-ICP-MS with thermal modelling to one of the best-studied crustal-scale igneous systems (Sesia Magmatic System, northern Italy), to quantify the relative contribution of crustal- and mantle-derived magmas in the entire system. Zircons from the deep gabbroic units define a tight range of $\epsilon$Hf (−2.5 ± 1.5). Granites and rhyolites overlap with this range but tail towards significantly more negative values (down to −9.5). This confirms that the entire system consists of hybrid magmas that stem from both differentiation of mantle-derived magmas and melting of the crust. Thermal modelling suggests that crustal melting and assimilation predominantly occurs during emplacement and evolution of magmas in the lower crust, although melt production is heterogeneous within the bodies both spatially and temporally. The spatial and temporal heterogeneity resolved by the thermal model is consistent with the observed Hf isotope variations within and between samples, and in agreement with published bulk-rock Sr–Nd isotopic data. On average, the crustal contribution to the entire system determined by mixing calculations based on Hf isotopic data range between 10 and 40%, even with conservative assumptions, whereas the thermal model suggests that this space- and time-averaged contribution does not exceed 20%. However, spatial and temporal variations in the crustal melt proportion (from 0 up to 80% as observed in the thermal model) may impart significant isotopic variability to different batches of magma observed on the outcrop scale, emphasizing the need to consider a magmatic system as a whole, i.e., by integrating all spatial and temporal scales, to more precisely quantify crustal growth vs. reworking.

**Keywords** Crustal-scale magmatic systems · Zircon Hf isotopes · Thermal modelling · Sesia Magmatic System · Granite formation · Crustal evolution

Communicated by Othmar Müntener.

<figure><img src="image_2.png" /></figure>

## 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

<footnote>Julian-Christopher Storck
julian.storck@geo.unibe.ch</footnote>

<footnote n="1">¹ Department of Earth Sciences, ETH Zurich, Clausiusstrasse 25, 8092 Zurich, Switzerland</footnote>

<footnote n="2">² Institute of Geological Sciences, University of Bern, Baltzerstrasse 1+3, 3012 Bern, Switzerland</footnote>

<footnote n="3">³ CNRS, Géosciences Environnement Toulouse, Observatoire Midi-Pyrénées, 14 av. E. Belin, 31400 Toulouse, France</footnote>

<footnote n="4">⁴ DST-Universita di Trieste, via Weiss 8, 34127 Trieste, Italy</footnote>
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).

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

<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 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</figcaption></figure>

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)
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⁻². The carrier gas consisted of high-purity helium (ca. 0.35 L·min⁻¹) and argon sample gas from the MC-ICP-MS (ca. 1.0 L·min⁻¹). 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): ¹⁷¹Yb (L4), ¹⁷³Yb (L2), ¹⁷⁵Lu (Ax),

<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 2** Plots of εHf (calculated at crystallization age, i.e., 283 Ma for the samples) as a function of ¹⁷⁶Yb/¹⁷⁷Hf 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 ¹⁷⁶Hf/¹⁷⁷Hf ratios and crystallization ages from Fisher et al. (2014). Error bars correspond to 2 S.E. uncertainties</figcaption></figure>

¹⁷⁶(Yb + Lu + Hf) (H1), ¹⁷⁷Hf (H2), ¹⁷⁸Hf (H3), ¹⁷⁹Hf (H4), ¹⁸⁰Hf (H5), ¹⁸¹Ta (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 ¹⁷³Yb/¹⁷¹Yb and ¹⁷⁹Hf/¹⁷⁷Hf and the isobaric interferences of ¹⁷⁶Yb and ¹⁷⁶Lu on ¹⁷⁶Hf were corrected using the natural abundance ratios of ¹⁷⁶Yb/¹⁷³Yb and ¹⁷⁶Lu/¹⁷⁵Lu, 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 ¹⁷⁶Hf/¹⁷⁷Hf ratios obtained on all reference materials are within uncertainties identical to the recommended values (Fig. 2A, see also Supplementary Table S1).

The initial ¹⁷⁶Hf/¹⁷⁷Hf ratios were calculated for each analysis using a ¹⁷⁶Lu decay constant of 1.865·10⁻¹¹ (Scherer

et al. 2001), the measured ¹⁷⁶Lu/¹⁷⁷Hf 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/m² 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 preheated 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 (L–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).

<figure><img src="image_7.png" /></figure>
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}}(\%) = \frac{\text{volume}_{\text{crustal melt}}^{t}}{\text{volume}_{\text{total melt}}^{t}} * 100 \quad (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}_{\text{t(end)}}^{\text{box, upper crust}} + \text{average crustal melting}_{\text{t(end)}}^{\text{overall system, lower crust}} \quad (2) $$

where, average crustal melting$_{t(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 $10^3$ 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 $\epsilon Hf = +4$, corresponding to the most primitive zircon $\epsilon Hf$ measured in this study; and $\epsilon 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 $\epsilon Nd$ value of +5 (Stille et al. 1996; Pinarelli et al. 2002) corresponding to $\epsilon 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 $\epsilon Nd$ of +7 (Monjoie 2004) corresponding to $\epsilon 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 $\epsilon 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 $\epsilon 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 $\epsilon Hf$ of the assimilated zircons. Therefore, we have used as bounds the range of average $\epsilon 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 $\epsilon Hf$ observed in our inherited zircons, whereby the main population ranges in $\epsilon Hf$ between -3 and -15 with a few grains with much lower $\epsilon Hf$ (-35 to -55). The range also compares well with average $\epsilon 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.

# 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 ¹⁷⁶Yb/¹⁷⁷Hf ratio for reference materials and samples (Fig. 2). This suggests that isobaric interferences of ¹⁷⁶Yb 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

<figure><img src="image_8.png" /><figcaption>**Fig. 3** Crustal-scale sketch of the Ivrea-Sesia magmatic system with expected sample locations in the reconstructed stratigraphy. Distribution histograms of single grain εHf₂₈₃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 εHf₂₈₃Ma values in 24 inherited zircon grains (LA-ICP-MS age > 300 Ma)</figcaption></figure>
(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).

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.
<figure><img src="image_9.png" /><figcaption>**Fig. 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)</figcaption></figure>

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).

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
A

LOWER CRUST

B

UPPER CRUST

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

C

<figure><img src="image_12.png" /></figure>

<table><thead><tr><th colspan="2">primary magma</th></tr><tr><th></th><th>[Sr]</th><th>Sri</th></tr></thead><tbody><tr><td>Min</td><td>91</td><td>0.7026</td></tr><tr><td>Max</td><td>188</td><td>0.7030</td></tr></tbody></table><table><thead><tr><th colspan="2">bulk crustal contaminant</th></tr><tr><th></th><th>[Sr]</th><th>Sri</th></tr></thead><tbody><tr><td>Min</td><td>137</td><td>0.7128</td></tr><tr><td>Max</td><td>296</td><td>0.7230</td></tr></tbody></table>

<figure><img src="image_13.png" /></figure>

<table><thead><tr><th colspan="2">primary magma</th></tr><tr><th></th><th>[Nd]</th><th>eNd</th></tr></thead><tbody><tr><td>Min</td><td>7.11</td><td>6.79</td></tr><tr><td>Max</td><td>9.65</td><td>7.18</td></tr></tbody></table><table><thead><tr><th colspan="2">bulk crustal contaminant</th></tr><tr><th></th><th>[Nd]</th><th>eNd</th></tr></thead><tbody><tr><td>Min</td><td>31.2</td><td>-11.22</td></tr><tr><td>Max</td><td>59.1</td><td>-7.04</td></tr></tbody></table>

**Fig. 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, overall with the kernel distributions). The two crustal melting values in

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

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)

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 Montjoie 2004) compared to the pyroxenite zircons (εHf of 0 to +4). Voshage et al. (1990) further argue that dramatically

<figure><img src="image_14.png" /></figure>
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 temporarily 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 over-
all 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 melt-
ing 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%), suggest-
ing 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-exist-
ing 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 assimi-
lation 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 subordi-
nate, 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 mod-
els tailing to 70–80%). The former value is probably more
representative, because it integrates variations of the crus-
tal 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 composi-
tion (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 estima-
tions 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 summa-
rize 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 mant-
tle source compositions, melting reactions and their spa-
tial-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 pre-
sented here were selected to be very conservative on pur-
pose, 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 unmetamor-
phosed, 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.% H₂O) 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 cal-
culations, 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 dif-
ficult, and would not allow those melts to participate in con-
structing 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 con-
sanguineous relationships of the different units proposed
based on single grain ID-TIMS age determinations (Karakas
et al. 2019). This suggests a continuous evolution, involv-
ing 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 gab-
broic 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 deliv-
ered 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 composi-
tional 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.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00410-021-01847-9.

**Acknowledgements** We thank J. Dufek for thermal modelling discussions, and P. Brack for support during fieldwork. In depth discussions with P. Ulmer, P. Brack, L. Tavazzani, and R. Economos on thermal modelling, crustal melting, and Ivrea magmatism are greatly appreciated. We acknowledge O. Müntener for editorial handling and constructive reviews by two anonymous reviewers that helped to sharpen the presentation of the main points of this manuscript. Funding of this project was provided to J.-C.S., O.K. and C.C.-M. by ETH Zürich. O.B. was supported by SNF Grant #200021_178928.

**Funding** Open Access funding provided by Universität Bern.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long

as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Annen C, Sparks RSJ (2002) Effects of repetitive emplacement of basaltic intrusions on thermal evolution and melt generation in the crust. Earth Planet Sci Lett 203:937–955. https://doi.org/10.1016/s0012-821x(02)00929-9

Annen C, Blundy JD, Sparks RSJ (2006) The genesis of intermediate and silicic magmas in deep crustal hot zones. J Petrol 47:505–539

Bachmann O, Miller CF, de Silva SL (2007) The volcanic-plutonic connection as a stage for understanding crustal magmatism. J Volcanol Geoth Res 167:1–23. https://doi.org/10.1016/j.jvolg.eores.2007.08.002

Barboza SA, Bergantz GW (1996) Dynamic model of dehydration melting motivated by a natural analogue: applications to the Ivrea-Verbano zone, northern Italy. Earth Environ Sci Trans R Soc Edinb 87(1–2):23–31

Barboza SA, Bergantz GW (1997) Melt productivity and rheology: complementary influences on the progress of melting. Numerical Heat Transfer, Part A Appl 31(4):375–392

Barboza SA, Bergantz GW (1998) Rheological transitions and the progress of melting of crustal rocks. Earth Planet Sci Lett 158(1–2):19–29

Barboza SA, Bergantz GW (2000) Metamorphism and anatexis in the mafic complex contact aureole, Ivrea Zone, Northern Italy. J Petrol 41(8):1307–1327

Barboza SA, Bergantz GW, Brown M (1999) Regional granulite facies metamorphism in the Ivrea zone: is the Mafic Complex the smoking gun or a red herring? Geology 27(5):447–450

Beard JS, Ragland PC, Crawford ML (2005) Reactive bulk assimilation: a model for crust-mantle mixing in silicic magmas. Geology 33:681–684. https://doi.org/10.1130/g21470.1

Bertotti G, Picotti V, Bernoulli D, Castellarin A (1993) From rifting to drifting—tectonic evolution of the south-alpine upper crust from the triassic to the early cretaceous. Sed Geol 86:53–76. https://doi.org/10.1016/0037-0738(93)90133-p

Bindeman IN, Fu B, Kita NT, Valley JW (2008) Origin and evolution of silicic magmatism at Yellowstone based on ion microprobe analysis of isotopically zoned zircons. J Petrol 49:163–193. https://doi.org/10.1093/petrology/egm075

Boriani A, Burlini L, Caironi V, Giobbi Origoni E, Sassi A, Sesana E (1988) Geological and petrological studies on the hercynian plutonism of Serie dei Laghi–geological map of its occurrence between Valsesia and Lago Maggiore (N-Italy). Rend Soc It Mineral Petrol 43(367):84

Bouilhol P, Schaltegger U, Chiaradia M, Ovtcharova M, Stracke A, Burg J-P, Dawood H (2011) Timing of juvenile arc crust formation and evolution in the Sapat Complex (Kohistan–Pakistan). Chem Geol 280:243–256

Bouvier A, Vervoort JD, Patchett PJ (2008) The Lu-Hf and Sm-Nd isotopic composition of CHUR: Constraints from unequilibriated chondrites and implications for the bulk composition
of terrestrial planets. Earth Planet Sci Lett 273:48–57. https://doi.org/10.1016/j.epsl.2008.06.010

Bowen NL (1928) The evolution of igneous rocks. Dover publications, New York

Brack P, Ulmer P, Schmid SM (2010) A crustal-scale magmatic system from the Earth's Mantle to the Permian surface. Swiss Bull Für Angewandte Geologie 15:3–21

Caricchi L, Simpson G, Schaltegger U (2014) Zircons reveal magma fluxes in the Earth's crust. Nature 511:457. https://doi.org/10.1038/nature13532

Carrigan CR (1988) Biot-number and thermos bottle effect- implications for magma-chamber convection. Geology 16:771–774. https://doi.org/10.1130/0091-7613(1988)016%3c0771:bnatbe%3e2.3.co;2

Cashman KV, Sparks RSJ, Blundy JD (2017) Vertically extensive and unstable magmatic systems: a unified view of igneous processes. Science. https://doi.org/10.1126/science.aag3055

Chelle-Michou C, Laurent O, Moyen JF, Block S, Paquette JL, Couzinié S, Gardien V, Vanderhaeghe O, Villaros A, Zeh A (2017) Pre-Cadoman to late-Variscan odyssey of the eastern Massif Central, France: formation of the West European crust in a nutshell. Gondwana Res 46:170–190

Chu NC et al (2002) Hf isotope ratio analysis using multi-collector inductively coupled plasma mass spectrometry: an evaluation of isobaric interference corrections. J Anal at Spectrom 17:1567–1574. https://doi.org/10.1039/b206707b

Clemens JD, Stevens G, Mayne MJ (2021) Do arc silicic magmas form by fluid-fluxed melting of older arc crust or fractionation of basaltic magmas? Contrib Miner Petrol 176(6):44

Collins WJ, Huang HQ, Jiang XY (2016) Water-fluxed crustal melting produces Cordilleran Batholiths. Geology 44:143–146. https://doi.org/10.1130/g37398.1

Couzinié S, Laurent O, Bouilhol P, Chelle-Michou C, Paquette JL, Moyen JF (2019) Detrital zircon U-Pb–Hf systematics of Ediacaran metasediments from the French Massif Central: consequences for the crustal evolution of the north Gondwana margin. Precambr Res 324:269–284

Davies JH (2013) Global map of solid Earth surface heat flow. Chem Geophys, Geosyst 14:4608–4622

Demarchi G, Quick JE, Sinigoi S, Mayer A (1998) Pressure gradient and original orientation of a lower-crustal intrusion in the Ivrea-Verbano Zone, northern Italy. J Geol 106:609–622

DePaolo DJ (1981) Trace-element and isotopic effects of combined wallrock assimilation and fractional crystallization. Earth Planet Sci Lett 53:189–202. https://doi.org/10.1016/0012-821x(81)90153-9

Devoir A, Bloch E, Müntener O (2021) Residence time of igneous garnet in Si-rich magmatic systems: Insights from diffusion modeling of major and trace elements. Earth Planet Sci Lett 560:116771

Dhuime B, Hawkesworth CJ, Storey CD, Cawood PA (2011) From sediments to their source rocks: Hf and Nd isotopes in recent river sediments. Geology 39(4):407–410

Dufek J, Bergantz GW (2005) Lower crustal magma genesis and preservation: a stochastic framework for the evaluation of basalt-crust interaction. J Petrol 46:2167–2195. https://doi.org/10.1093/petrology/egi049

Ewing TA, Rubatto D, Hermann J (2014) Hafnium isotopes and Zr/Hf of rutile and zircon from lower crustal metapelites (Ivrea-Verbano Zone, Italy): implications for chemical differentiation of the crust. Earth Planet Sci Lett 389:106–118

Farina F, Stevens G (2011) Source controlled Sr-87/Sr-86 isotope variability in granitic magmas: the inevitable consequence of mineral-scale isotopic disequilibrium in the protolith. Lithos 122:189–200. https://doi.org/10.1016/j.lithos.2011.01.001

Fountain DM (1976) Ivrea-Verbano and Strona-Ceneri zones, northern Italy–cross section of continental crust–new evidence from seismic velocities of rock samples. Tectonophysics 33:145–165

Garuti G, Bea F, Zaccarini F, Montero P (2001) Age, geochemistry and petrogenesis of the ultramafic pipes in the Ivrea Zone, NW Italy. J Petrol 42(2):433–457

Geist D, Howard KA, Larson P (1995) The generation of oceanic rhyolites by crystal fractionation—the basalt-rhyolite association at Volcan-Alcedo, Galapagos-Archipelago. J Petrol 36:965–982. https://doi.org/10.1093/petrology/36.4.965

Guergouz C, Martin L, Vanderhaeghe O, Thébaud N, Fiorentini M (2018) Zircon and monazite petrochronologic record of prolonged amphibolite to granulite facies metamorphism in the Ivrea-Verbano and Strona-Ceneri Zones, NW Italy. Lithos 308–309:1–18

Handy MR (1987) The structure, age and kinematics of the pogallo fault zone—southern alps, northwestern Italy. Eclogae Geol Helv 80:593

Handy MR, Franz L, Heller F, Janott B, Zurbriggen R (1999) Multistage accretion and exhumation of the continental crust (Ivrea crustal section, Italy and Switzerland). Tectonics 18(6):1154–1177

Hellstrom J, Paton C, Woodhead J, Hergt J (2008) Iolite: software for spatially resolved LA-(quad and MC) ICPMS analysis. Mineral Assoc Can Short Course Series 40:343–348

Henk A, Franz L, Teufel S, Oncken O (1997) Magmatic underplating, extension, and crustal reequilibration: insights from a cross-section through the Ivrea Zone and Strona-Ceneri Zone, northern Italy. J Geol 105(3):367–378

Hildreth W (2004) Volcanological perspectives on Long Valley, Mammoth Mountain, and Mono Craters: several contiguous but discrete systems. J Volcanol Geoth Res 136:169–198. https://doi.org/10.1016/j.volgrores.2004.05.019

Hildreth W, Moorbath S (1988) Crustal contributions to arc magmatism in the Andes of central Chile. Contrib Miner Petrol 98:455–489. https://doi.org/10.1007/bf00372365

Jacob JB, Moyen JF, Fiannacca P, Laurent O, Bachmann O, Janousek V, Farina F, Villaros A (2021) Crustal melting vs fractionation of basaltic magmas: part 2, attempting to quantify mantle and crustal contributions in granitoids. Lithos. https://doi.org/10.1016/j.lithos.2021.106292

Jagoutz O, Bouilhol P, Schaltegger U, Müntener O (2019) The isotopic evolution of the Kohistan Ladakh arc from subduction initiation to continent arc collision. Geol Soc, London, Spl Publ 483(1):165–182

Karakas O, Dufek J (2015) Melt evolution and residence in extending crust: thermal modeling of the crust and crustal magmas. Earth Planet Sci Lett 425:131–144

Karakas O, Degruyter W, Bachmann O, Dufek J (2017) Lifetime and size of shallow magma bodies controlled by crustal-scale magmatism. Nat Geosci 10:446. https://doi.org/10.1038/ngeo2959

Karakas O, Wotzlaw J-F, Guillong M, Ulmer P, Brack P, Economos R, Bergantz GW, Sinigoi S, Bachmann O (2019) The pace of crustal-scale magma accretion and differentiation beneath silicic caldera volcanoes. Geology 47(8):719–723

Kelsey DE, Clark C, Hand M (2008) Thermobarometric modelling of zircon and monazite growth in melt-bearing systems: examples using model metapelite and metapsammitic granulites. J Metamorph Geol 26:199–212. https://doi.org/10.1111/j.1525-1314.2007.00757.x

Kunz BE, Regis D, Engi M (2018) Zircon ages in granulite facies rocks: decoupling from geochemistry above 850 °C? Contrib Mineral Petrol 173:26

L'ecuyer P (1988) Efficient and portable combined random number generators. Commun ACM 31:742–751
Matsumoto A, Nakagawa M (2010) Formation and evolution of silicic magma plumbing system: petrology of the volcanic rocks of Usu volcano, Hokkaido. Jpn J Volcanol Geoth Res 196:185-207. https://doi.org/10.1016/j.jvolgeores.2010.07.014

Mintrone M, Galli A, Laurent O, Chelle-Michou C, Schmidt MW (2020) Quantifying frozen melt in crustal rocks: a new melt-o-meter based on zircon rim volumes. Chem Geol 551:119755

Monjoie P (2004) The Mont Collon mafic complex (Austroalpine Dent Blanche nappe). Permian evolution of the Western European mantle (Doctoral dissertation, Université de Lausanne)

Morel MLA, Nebel O, Nebel-Jacobsen YJ, Miller JS, Vroon PZ (2008) Hafnium isotope characterization of the GJ-1 zircon reference material by solution and laser-ablation MC-ICPMS. Chem Geol 255:231-235. https://doi.org/10.1016/j.chemgeo.2008.06.040

Moyen JF (2020) Granites and crustal heat budget. Geol Soc London, Spl Publ 491:77-100

Moyen JF, Janousek V, Laurent O, Bachmann O, Jacob JB, Farina F, Fiannacca P, Villaros A (2021) Crustal melting vs. fractionation of basaltic magmas: part 1, the bipolar disorder of granite petrogenetic models. Lithos. https://doi.org/10.1016/j.lithos.2021.106291

Peressini G, Quick JE, Sinigoi S, Hofmann AW, Fanning M (2007) Duration of a large mafic intrusion and heat transfer in the lower crust: a SHRIMP U-Pb zircon study in the Ivrea-Verbano Zone (Western Alps, Italy). J Petrol 48(6):1185-1218

Pin C (1990) Evolution of the lower crust in the ivrea zone—a model based on isotopic and geochemical data. Granulites Crustal Evol 311

Pinarelli L, Moro AD, Boriani A, Caironi V (2002) Sr, Nd isotope evidence for an enriched mantle component in the origins of the Hercynian gabbro-granite series of the “Serie dei Laghi” (Southern Alps, NW Italy). Eur J Mineral 14(2):403-415

Plank T, Kelley KA, Zimmer MM, Hauri EH, Wallace PJ (2013) Why do mafic arc magmas contain [sim] 4wt% water on average? Earth Planet Sci Lett 364:168-179

Price RC, Gamble JA, Smith IEM, Stewart RB, Eggins S, Wnight IC (2005) An integrated model for the temporal evolution of andesites and rhyolites and crustal development in New Zealand's North Island. J Volcanol Geoth Res 140:1-24. https://doi.org/10.1016/j.jvolgeores.2004.07.013

Quick JE, Sinigoi S, Negrini L, Demarchi G, Mayer A (1992) Syn-magmatic deformation in the underplated igneous complex of the Ivrea-Verbano zone. Geology 20(7):613-616

Quick JE, Sinigoi S, Mayer A (1994) Emplacement dynamics of a large mafic intrusion in the lower crust, Ivrea-Verbano Zone, northern Italy. J Geophys Res: Solid Earth 99(B11):21559-21573

Quick JE, Sinigoi S, Mayer A (1995) Emplacement of mantle peridotite in the lower continental crust, Ivrea-Verbano zone, northwest Italy. Geology 23(8):739-742

Quick J, Sinigoi S, Peressini G, Demarchi G, Wooden J, Sbisà A (2009) Magmatic plumbing of a large Permian caldera exposed to a depth of 25 km. Geology 37:603-606

Rioux M, Hacker B, Mattinson J, Kelemen P, Blusztajn J, Gehrels G (2007) Magmatic development of an intra-oceanic arc: high-precision U-Pb zircon and whole-rock isotopic analyses from the accreted Talkeetna arc, south-central Alaska. Geol Soc Am Bull 119:1168-1184. https://doi.org/10.1130/b25964.1

Rivalenti G (1975) The origin of the Ivrea-Verbano basic formation (western Italian alps)-whole rock geochemistry

Rivalenti G, Rossi A, Siena F, Sinigoi S (1984) The layered series of the Ivrea-Verbano igneous complex, we tern Alps, Italy. Tschermaks Mineralogische Und Petrographische Mitteilungen 33(2):77-99

Rudnick RL, Gao S (2014) Composition of the continental crust. Trea-tise Geochemi 4:1-51

Schaltegger U, Brack P (2007) Crustal-scale magmatic systems during intracontinental strike-slip tectonics: U, Pb and Hf isotopic constraints from Permian magmatic rocks of the Southern Alps. Int J Earth Sci 96:1131-1151. https://doi.org/10.1007/s00531-006-0165-8

Schaltegger U, Gebauer D (1999) Pre-Alpine geochronology of the Central, Western and Southern Alps Schweizerische. Mineral Petrogr Mitt 79:79-87

Scherer E, Munker C, Mezger K (2001) Calibration of the lutetium-hafnium clock. Science 293:683-687. https://doi.org/10.1126/science.1061372

Schmid R (1967) Zur Petrographie und Struktur der Zone Ivrea-Verbano zwischen Valle d'Ossola und Val Grande (Prov. Novara Italien). Schweiz Mineral Petrogr Mitt 47:935-1117

Schnetger B (1994) Partial melting during the evolution of the amphibolite- to granulite-facies gneisses of the Ivrea Zone, northern Italy. Chem Geol 113:71-101

Sills JD, Tarney J (1984) Petrogenesis and tectonic significance of amphibolites interlayered with metasedimentary gneisses in the Ivrea Zone, Southern Alps, Northwest Italy. Tectonophysics 107(3-4):187-206

Sinigoi S, Antonini P, Demarchi G, Longinelli A, Mazzucchelli M, Negrini L, Rivalenti G (1991) Intractions of mantle and crustal magmas in the southern part of the Ivrea Zone (Italy). Contrib Miner Petrol 108(4):385-395

Sinigoi S et al (1994) Chemical evolution of a large mafic intrusion in the lower crust, Ivrea-Verbano Zone, northern Italy. J Geophys Res: Solid Earth 1978-2012(99):21575-21590

Sinigoi S, Quick JE, Mayer A, Budahn J (1996) Influence of stretching and density contrasts on the chemical evolution of continental magmas: an example from the Ivrea-Verbano zone. Contrib Miner Petrol 123:238-250. https://doi.org/10.1007/s004100050153

Sinigoi S, Quick JE, Demarchi G, Klötzli U (2011) The role of crustal fertility in the gen- eration of large silicic magmatic systems trig- gered by intrusion of mantle magma in the deep crust. Contrib Mineral Petrol 162:691-707. https://doi.org/10.1007/s00410-011-0619-2

Sinigoi S, Quick JE, Demarchi G, Klotzli US (2016) Production of hybrid granitic magma at the advancing front of basaltic under-plating: inferences from the Sesia Magmatic System (south-western Alps, Italy). Lithos 252:109-122. https://doi.org/10.1016/j.lithos.2016.02.018

Slama J et al (2008) Plesovice zircon—a new natural reference material for U-Pb and Hf isotopic microanalysis. Chem Geol 249:1-35. https://doi.org/10.1016/j.chemgeo.2007.11.005

Soder CG, Romer RL (2018) Post-collisional potassic-ultrapotassic magmatism of the Variscan Orogen: implications for mantle metasomatism during continental subduction. J Petrol 59(6):1007-1034

Stille P, Schaltegger U, Basu A, Hart S (1996) The lithospheric mantle beneath Central Europe: Nd isotopic constraints for its Late Proterozoic enrichment and implications for early crustal evolution. Geophys Monograph-Am Geophys Union 95:269-276

Storck JC, Wotzlaw JF, Karakas O, Brack P, Gerdes A, Ulmer P (2020) Hafnium isotopic record of mantle-crust interaction in an evolving continental magmatic system. Earth Planetary Sci Lett. https://doi.org/10.1016/j.epsl.2020.116100

Tavazzani L, Peres S, Sinigoi S, Demarchi G, Musumeci G (2017) Structure and petrography of the Valle Mosso pluton, Sesia magmatic system, Southern Alps. J Maps 13(2):684-697

Tavazzani L, Peres S, Sinigoi S, Demarchi G, Economos RC, Quick JE (2020) Timescales and mechanisms of crystal-mush rejuvenation and melt extraction recorded in Permian plutonic and volcanic rocks of the Sesia magmatic system (Southern Alps, Italy). J Petrol 61(5):egaa049
Taylor HP (1980) The effects of assimilation of country rocks by mag-
mas on O-18 O-16 and Sr-87 Sr-86 systematics in igneous rocks.
Earth Planetary Sci Lett 47:243–254. https://doi.org/10.1016/
0012-821x(80)90040-0

Thompson AB, Matile L, Ulmer P (2002) Some thermal constraints on crustal assimilation during fractionation of hydrous, mantle-derived magmas with examples from central alpine batholiths. J Petrol 43:403–422. https://doi.org/10.1093/petrology/43.3.403

Turcotte DL, Schubert G (2002) Geodynamics. Cambridge University Press

Vervoort JD, Patchett PJ, Blichtert-Toft J, Albarede F (1999) Relation-
ships between Lu-Hf and Sm-Nd isotopic systems in the global
sedimentary system. Earth Planet Sci Lett 168:79–99. https://doi.
org/10.1016/s0012-821x(99)00047-3

Vielzeuf D, Holloway JR (1988) Experimental determination of the fluid-absent melting relations in the pelitic system. Contrib Miner Petrol 98:257–276

Villaros A, Buick IS, Stevens G (2012) Isotopic variations in S-type granites: an inheritance from a heterogeneous source? Contrib Mineral Petrol 163(2):243–257

Voller V, Swaminathan C (1991) General source-based method for solidification phase change. Numerical Heat Transfer 19B:175–189

Voshage H, Hofmann A, Mazzucchelli M, Rivalenti G, Sinigoi S, Raczek I, Demarchi G (1990) Isotopic evidence from the Ivrea

Zone for a hybrid lower crust formed by magmatic underplating.
Nature 347:731–736

Walker BA, Bergantz GW, Otamendi JE, Ducea MN, Cristofolini EA (2015) A MASH zone revealed: the mafic complex of the sierra valle fertil. J Petrol 56(9):1863–1896

Wiedenbeck MAPC et al (1995) Three natural zircon standards for U-Th-Pb, Lu-Hf, trace element and REE analyses. Geostand Newsl 19:1–23

Woodhead JD, Hergt JM (2005) A preliminary appraisal of seven natu-
ral zircon reference materials for in situ Hf isotope determination.
Geostand Geoanal Res 29:183–195

Zheng YF, Gao P (2021) The production of granitic magmas through crustal anatexis at convergent plate boundaries. Lithos 106232

Zingg A, Handy MR, Hunziker JC, Schmid SM (1990) Tectono-
metamorphic history of the Ivrea Zone and its relationship to
the crustal evolution of the Southern Alps. Tectonophysics
182(1-2):169–192

Zurbriggen R, Franz L, Handy MR (1997) Pre-Variscan deformation, metamorphism and magmatism in the Strona-Ceneri zone (southern Alps of northern Italy and southern Switzerland). Schweiz Miner Petrogr Mitt 77:361–380

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.