Establishing a protocol for the selection of zircon inclusions in garnet for Raman thermobarometry

Zircon is one of the most common accessory minerals in various igneous, sedimentary, and metamorphic rocks. Furthermore, due to its large stability field and its physical robustness, zircon often hosts ultrahigh pressure (UHP) metamorphic minerals such as coesite and diamond (Parkinson and Katayama 1999; Ye et al. 2000b; Hermann et al. 2001; Rubatto and Hermann 2007; Katayama and Maruyama 2009). Zircon can be used for a wide range of applications, including U-Pb geochronology, and potential use as a host phase for the disposal of excess weapons-grade Pu because of its capacity to incorporate radioactive elements, (e.g., Ewing et al. 1995). Further, radiation-damaged zircon has been successfully applied as a model system to elucidate the behavior of partially ordered structures at high pressure and temperature conditions (e.g., Colombo et al. 1999; Binvignat et al. 2018).

A recent development in zircon petrology is its exploitation as a mineral inclusion in the frame of Raman elastic thermobarometry (Campomenosi et al. 2018; Stangarone et al. 2019; Zhong et al. 2019). The main idea is that the geochemical and isotopic information that can be obtained from zircon is complemented by elastic thermobarometry of zircon and its host to retrieve P-T-time-fluid-deformation paths of metamorphic and igneous rocks. For this new application, however, zircon composition and structural metamictization due to the decay of radioactive elements such as U and Th within its crystal structure has to be taken into account, because the radiation-induced structural changes can lead to large variations in the Raman spectrum. Indeed, the Raman peaks of radiation-damaged crystals shift toward lower wavenumbers, while the peak width increases with respect to pristine crystals taken as reference (Colombo et al. 1999; Ríos et al. 2000; Nasdala et al. 2001; Geisler and Pidgeon 2002; Binvignat et al. 2018). In this regard, useful diagrams such as peak width (Γ), in terms of full-width at half maximum, vs. peak position of the major Raman peak near 1008 cm⁻¹ (B_1g phonon mode related to antisymmetric SiO₄ stretching) can be exploited to recognize the zircon crystals whose Raman spectra is unaffected by metamictization processes and that can be safely used for elastic thermobarometry purposes (e.g., Zhong et al. 2019). However, even though the degree of metamictization is the main factor affecting the Raman spectra of zircon, it is not the only one. Cathodoluminescence (CL) studies of zircon internal zoning (e.g., Rubatto and Gebauer 2000; Corfu et al. 2003) reveal that single crystals of zircon, particularly in metamorphic rocks, commonly show a complex internal chemical and structural heterogeneity at the micrometer-scale. Furthermore, variations in the Raman peak positions can be caused also by chemical substitution: a typical example is Hf replacing Zr at the dodecahedral crystallographic site (e.g., Hoskin and Rodgers 1996).

To address these problems, we have chosen a particularly well-known sample suite from the Dora Maira UHP unit (Western Alps, Italy), where garnet megablasts host abundant zircon inclusions. We report micrometer-scale structural and chemical information of partially exposed zircon crystals studied by complementary Raman spectroscopy, charge contrast (CC) imaging, and laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS). The purpose is (1) to determine quantitatively how the structural and chemical heterogeneities of zircons can affect their Raman spectrum and (2) to propose a best-practice protocol for the selection of zircon inclusions suitable for elastic thermobarometric estimations.

We selected exposed and buried zircon inclusions in garnet megablasts from the UHP whiteschists of the Dora Maira Massif. The investigated samples come from the two major outcrops of the Brossasco-Isasca Unit (e.g., Chopin 1984): Vallone Gilba and Vallone Martiniana close to the Case Fapina and Case Parigi localities, respectively. Several polished sections of about 350 μm in thickness have been prepared from three different garnet megablasts (up to 15–20 cm in size). Zircon crystals are included in garnet and kyanite and range in size from ~10 to 200 μm; they are either single, isolated crystals or clusters of two or more inclusions. The shape of the inclusions varies from rounded to elongated and almost idiomorphic with sharp corners and edges. Some of the zircon crystals also host mineral inclusions such as rutile, coesite, sheet silicates, and apatite, in agreement with previous observations (e.g., Gautiez-Putallaz et al. 2016).

Previous results (Gebauer et al. 1997; Gautiez-Putallaz et al. 2016) suggested the existence of two main zircon generations that can coexist within the same crystal in these rocks. The first one records the crystallization of the Permian protolith (~275 Ma), while the second domain is Eocene in age (35.1 ± 0.8 Ma) and related to metamorphism during the Alpine subduction. Chopin (1984) was the first to constrain this metamorphic event at UHP conditions by finding coesite inclusions within garnet. Optically single crystals of coesite inclusions have been found also in this study, especially within the garnet neoblasts (SiO₂ saturated whiteschists). Within the garnet megablasts, rare coesite inclusions have been found only at the external rim. More recent studies have constrained the metamorphic peak conditions in the diamond stability field at about 4–4.3 GPa and 730–750 °C (e.g., Hermann 2003; Ferrando et al. 2009; Gautiez-Putallaz et al. 2016).

The selection of completely buried zircon inclusions follows the protocol given by Campomenosi et al. (2018) for elastic thermobarometry applications (i.e., small, similar in size, non-fractured, and isolated inclusions within the center of a non-fractured host). On the other hand, the selection of the exposed grains (about 60% of the total number of zircon inclusions present in the samples) was based on criteria such as variable size of the inclusion and of the surrounding intact host (the effective host), exposure degree of the inclusions (optically estimated in terms of the ratio between the exposed and buried inclusion surface) and presence of fractures or any other kind of discontinuity in both the host and the inclusion.

Exposed zircon inclusions were investigated by CC imaging first and Raman spectroscopy later to obtain textural and structural information, and finally by LA-ICP-MS for the chemical composition.

CC images were acquired with a ZEISS EV050 scanning electron microscope at the Institute of Geological Sciences, University of Bern, at low vacuum conditions (18 Pa), 12 kV, a beam current of 100 mA and a working distance of 9.5 mm. CC images were obtained from well-polished, uncoated samples. It has been demonstrated that CC images correlate exactly to cathodoluminescence images and result from the complex interaction between the electron beam, the positive ions generated by electron-gas interactions in the chamber, a biased detector, and the sample (Griffin 2000; Watt et al. 2000). Internal check in the Bern laboratory confirmed that CC images are identical to panchromatic cathodoluminescence images.

Raman measurements were collected with two different Horiba Jobin-Yvon spectrometers: a T64000 triple-monochromator system operating in a subtractive mode and an Explora_Plus single-monochromator spectrometer at the University of Hamburg and Genova, respectively. The first one was equipped with a symphony LN₂-cooled CCD detector, 1800 gr/mm holographic gratings, an Olympus BH-41 optical microscope and a Coherent Ar⁺ laser, whereas the second with a Peltier-cooled CCD detector, a 2400 gr/mm grating, an Olympus BH-41 optical microscope, and Nd:YAG solid state laser. In both cases, using a slit width of 100 μm, the spectral resolution was approximately 2 cm⁻¹ (determined using the photoluminescence line of Sm: SrB₄O₇ for the T64000 and a Ne lamp for the Explora_Plus spectrometers), and therefore apparatus corrections of the peak widths were not required (Nasdala et al. 2001). Sufficient acquisition times and a surface density of the laser power were used to avoid sample overheating during the Raman measurements (Zhong et al. 2019). The instrumental accuracy in the peak position determination was about 0.35 cm⁻¹ for the T64000 spectrometer and about 0.55 cm⁻¹ for the Explora_Plus. These values refer to the half of the spectral pixel-pixel distance, which results from several instrumental factors, including the CCD detector array, groove density of gratings, and excitation-light wavelength.

However, both the instruments were manually calibrated to the silicon Raman peak at 520.5 cm⁻¹ at each measurement session, and the stability of the instrumental setup was double-checked measuring reference crystals every 4 h. The resulting standard deviation in the determination of the silicon peak along different days of measurements was approximately 0.04 cm⁻¹. Furthermore, a crosscheck of the two data set measured with the two different spectrometers has been performed on several zircon grains to verify their consistency. Finally, the lateral spatial resolution of the two instruments is approximately 1 μm and the probed volume at the sample surface is of about 1 μm cube.

The Origin Lab-Pro 2018 software package was used for data evaluation. The collected spectra were baseline corrected for the continuum luminescence background when necessary, temperature-reduced to account for the Bose-Einstein occupation factor (Kuzmany 2009), and normalized to the acquisition time (7 s). Peak positions, FWHMs, and integrated intensities were determined from fits with pseudo-Voigt functions [PV = (1 – q)·Lorentz + q·Gauss, q is the weight coefficient]. The criterion for the maximum number of fitted peaks was ΔI < I/2, where I and ΔI are the calculated magnitude and uncertainty of each peak intensity, respectively.

Trace element analyses were obtained with an ASI-RESOlution 193 nm laser system coupled to an Agilent 7900 quadrupole Inductively Coupled Plasma mass spectrometer (ICPMS) at the Institute of Geological Sciences, University of Bern. The laser was tuned to a repetition rate of 5 Hz and an energy output of 4 mJ (corresponding to an HV of about 26–27 kV). The ICPMS was tuned for maximum sensitivity and minimum production of molecular species, maintaining ThO⁺/Th⁺ at <0.2%. The spot size on the zircon crystal was 20 μm. Analyses were standardized to glass NIST 612, and zircon 91500 was run as a secondary standard to monitor accuracy. Stoichiometric Si was employed as an internal standard for zircon (SiO₂: 31.6 wt%). Reproducibility and accuracy were within 10% or less across all analyzed elements. The data were reduced with the free-ware Iolite (Paton et al. 2011) and its data reduction scheme for trace elements (Woodhead et al. 2007).

CC imaging of 64 inclusions from samples DM17-13 and DM17-49 (Martiniana locality) and 83 grains from samples DMG4-6 and DM17-35 (Gilba locality) shows large variability in zircon internal texture. By combining the data on various crystals, it was possible to recognize at least four major domains (Fig. 1): (1) a dark (i.e., low CC-emission) domain, commonly corresponding to the crystal cores, (2) an oscillatory-zoned domain (medium-high CC emission), which usually corresponds to entire crystals or the core of elongated crystals and are further subdivided into oscillatory-dark and oscillatory-bright domains based on CC emission, (3) a transition (or undefined) domain, whose appearance is between the darker and a brighter domains and without a defined internal texture, and (4) a bright domain (i.e., high CC-emission), which usually belongs to the external rim of the crystals. About 70% of zircon inclusions from the Gilba locality show a more elongated shape with a typical oscillatory-zoned core surrounded by 1 or 2 thin brighter rims. On the other hand, about 80% zircon inclusions from the Martiniana locality present more sub-idiomorphic shapes and more homogeneous internal texture. Furthermore, as already reported by Gautiez-Putallaz et al. (2016) for zircon hosted in Si undersaturated whiteschists, some of the analyzed inclusions in Martiniana, where garnet is more abundant, display oscillatory dark domains.

Zircon has tetragonal symmetry with space group I4₁/amd. According to group-theory analysis, the optical phonons at the Brillouin-zone center of zircon are (Kroumova et al. 2010):

The A_1g, B_1g, B_2g, and E_g modes are Raman-active and therefore a total of 12 Raman peaks can be observed in the spectrum of a randomly oriented zircon. According to previous studies (Knittle and Williams 1993; Nasdala et al. 2001; Binvignat et al. 2018), both pressure and metamictization affect most strongly the B_1g mode near 1008 cm⁻¹, originating from anti-symmetrical SiO₄ stretching. For this reason, our discussion is mainly focused on this phonon mode.

The Raman spectroscopic analysis is based on 250 spectra collected from 60 zircon inclusions in different CC domains. As an example, Figure 1 shows the variation of the internal SiO₄ antisymmetric stretching mode (B_1g) near to 1008 cm⁻¹ of representative zircon inclusions exhibiting heterogeneous CC emission. Generally, dark CC domains correspond to Raman peak broadening and shift toward lower wavenumbers when compared to those collected across brighter or oscillatory-zoned domains. Furthermore, transition CC domains show transitional average features also for the Raman spectra.

A more quantitative evaluation of the structural state of zircon crystals can be done plotting the relationship between the phonon wavenumber ω and the corresponding peak widht Γ (Geisler and Pidgeon 2002; Zhong et al. 2019). The data measured on partially exposed zircon inclusions (Fig. 2) show two different trends that correlate with the four different zircon domains (bright, oscillatory-zoned, undefined, and dark).

The first trend is characterized by ω < 1008 cm⁻¹, 4.5 cm⁻¹ ≤ Γ < 12 cm⁻¹, and ω(Γ) with a negative slope (inverse correlation). The data points on this trend mostly correspond to zircon domains with dark and oscillatory-dark CC emission.

The second trend, at higher wavenumbers, is characterized by 1008 cm⁻¹ < ω < 1011 cm⁻¹, 3.5 cm⁻¹ < Γ < 5.5 cm⁻¹, and ω(Γ) ≈ constant. This trend is mostly defined by zircon domains with bright and oscillatory(-bright) CC emission.

Note that, on average, oscillatory-dark domains display the largest Γ and lowest wavenumber with respect to oscillatory-bright domains. On the other hand, undefined domains are difficult to classify since they scatter over the entire data.

The different types of zircon domains also show differences in trace-element composition. Tables 1a–1b report representative major, minor, and trace elements composition of analyzed zircon crystal inclusions as determined by LA-ICP-MS with the associated CC-emission domains. Hf concentration is usually around 11 000 μg/g and only few analyses show higher values up to 14 000 μg/g. U content varies significantly between oscillatory-bright zoned, oscillatory-dark, and dark CC domains, with the last showing, usually, the highest concentration (see for example DM17-35-3a-zrc10-p1 and p2 or even DM17-35-3a-zrc14-p3 and p2 in Table 1a). U contents vary from a few hundred μg/g in the bright domains to 5000 μg/g for the darker domains and generally show an inverse relationship with CC emission, indicating that U suppresses luminescence (e.g., Rubatto and Gebauer 2000).

In line with previous REE data sets (e.g., Gautiez-Putallaz et al. 2016), zircon cores with oscillatory or convolute zoning have a steep HREE-enriched pattern with a pronounced negative Eu anomaly (e.g., Fig. 3c p1). This REE pattern is characteristic of conditions where plagioclase was present and garnet was absent, i.e., during the crystallization of the granitic protholith in Permian times (e.g., Gautiez-Putallaz et al. 2016).

The rims and many crystals with oscillatory zoning have lower REE contents and a flat HREE pattern with no or weak Eu-anomaly (Fig. 3). This REE pattern is diagnostic for zircons formed during prograde to peak metamorphism, where feldspar is no longer stable (e.g., Gautiez-Putallaz et al. 2016). The variable slope of the HREE patterns ranging from slightly positive to even negative (Fig. 3c) indicates zircon growth together with garnet, as both minerals compete for the HREE. Therefore a first distinction between pre-Alpine inherited igneous zircon and Alpine metamorphic zircon can be obtained from the trace element analysis. Notably, several large crystals that display a fine oscillatory to sector zoning more or less surrounded by a dark or light rim present a depletion in HREE across the entire grain (e.g., Figs. 3d and 3e). This demonstrates that oscillatory zoning occurs also in Alpine metamorphic zircon and is not restricted to inherited igneous zircon, demonstrating that both CC and trace elements are required for the correct classification of the zircon.

Previous analyses on the same locality demonstrated that inherited zircon cores are ca. 275 Ma old and metamorphic rims are 35.1 ± 0.8 Ma old (Gebauer et al. 1997; Gauthiez-Putallaz et al. 2016). In our case, a rough estimate of the age can also be obtained from the LA-ICP-MS data, even if proper age standardization and common Pb correction have not been performed. The ²⁰⁶Pb/²³⁸U values estimated in such a way are in general agreement with those obtained before on zircon from the same locality and correlate with the distinction made from the REE patterns (Gauthiez-Putallaz et al. 2016). This rough age estimate allows distinguishing between Alpine domains (usually lower U content and bright or oscillatory-bright CC domains) and pre-Alpine domains (higher U content and dark or oscillatory-zoned CC domains at the zircon cores) especially in cases where REE patterns might be ambiguous (see Table 2). Notably the large crystals with oscillatory-sector zoning and low HREE content are of Alpine age.

Assuming an age of 275 Ma for the pre-Alpine and 35 Ma for the Alpine domains (Gebauer et al. 1997; Gauthiez-Putallaz et al. 2016) and based on the measured U-Th composition of the zircon domains, their accumulation radiation doses D_a was calculated as follows (e.g., Nasdala et al. 2001):

where λ_n and m_n are, respectively, the nuclear decay constant and mass of the corresponding isotope (²³⁸U, ²³⁵U, and ²³²Th) (Steiger and Jaeger 1977), the coefficient n represents the number of α decays per nucleus (n₂₃₈ = 8, n₂₃₅ = 7, and n₂₃₂ = 6) and N_A is the Avogadro's number c_U and c_Th are the measured concentrations of U and Th, respectively, and t the time. Note that Equation 1 presumes a U isotopic composition of 99.28% ²³⁸U and 0.72% ²³⁵U. The resulting D_a values are reported in Table 2.

The Raman spectra collected on the completely buried zircon crystals have higher wavenumbers (from 1010 cm⁻¹ up to 1013.5 cm⁻¹) with respect to the partially exposed grains, with a Γ mostly between 3.5 and 5.5 cm⁻¹. Only a few analyses reach values up to 6.5 cm⁻¹ (Fig. 4). The overlap between the two sets of inclusions is limited to about 10% of the analyses.

The data from the exposed and fully enclosed zircon inclusions can also be analyzed in terms of the wavenumber difference (Δω) between the B_1g mode near 1008 cm⁻¹ and the A_1g mode near 440 cm⁻¹ rather than the absolute wavenumber. In this way, since the data were collected along different sessions of measurements, we can avoid eventual effects due to instrumental drift. As shown previously by Knittle and Williams (1993) and more recently by Binvignat et al. (2018) both phonon modes show an increase in the phonon wavenumber as the hydro-static pressure increases. However, as reported from the same authors, the A_1g mode shows a weaker pressure dependence with respect to the B_1g mode and, therefore, to an increase of the pressure acting on the crystal should correspond an increase in the value of Δω (Fig. 5).

A large number of non-formula elements can be incorporated in the zircon crystal structure, however, most of them are usually far below 1 wt% (e.g., Hoskin and Schaltegger 2003), with the notable exception of Hf. Raman data measured on isomorphic series ZrSiO₄–HfSiO₄, suggest that even if as much as 25% of all Zr ions are replaced by Hf the frequency variation of the main peaks does not exceed 3 cm⁻¹ (Hoskin and Rodgers 1996). Furthermore, Nasdala et al. (2002) suggested that positions and Γ of the main peaks in annealed, metamict gemstone-quality zircon, containing up to 6000 µg/g of U and 16 300 µg/g of Hf, deviate less than 1 cm⁻¹ from the data of pure well-crystalline ZrSiO₄. Therefore, since our exposed zircon inclusions contain impurities and trace elements below these values (maximum U content is 5060 µg/g and maximum Hf content is 14 570 µg/g, see Table 1), the effect of chemical variations on zircon Raman spectra is negligible.

Zircon variation in HREE composition (i.e., depletion in the Alpine domains) has been previously interpreted as the result of growth zoning during metamorphism in a fractionating bulk composition where REE are largely incorporated in garnet (e.g., Gautiez-Putallaz et al. 2016). Indeed, the garnet host shows a similar REE pattern along a core-to-rim line profile (e.g., Gautiez-Putallaz et al. 2016), indicating equilibrium conditions with the associated zircon grains. The qualitative ²⁰⁶Pb/²³⁸U measurements confirm that zircon crystals showing a flat HREE pattern throughout the entire grain are completely metamorphic in origin. In this regard, it is important to note that, although REE chemical zonation does not influence the main Raman scattering features, such as peak broadening and position, it can give rise to heterogeneous photoluminescence, and hence to different background levels of the Raman spectra collected from zones. Besides, depending on the excitation laser wavelength, additional photoluminescence peaks may be observed next to the fundamental Raman peaks of zircon.

Charge contrast or cathodoluminesce imaging of metamorphic and inherited zircons usually show notable differences in terms of the corresponding internal texture (e.g., Rubatto and Gebauer 2000; Corfu et al. 2003). Oscillatory-zoned grains have been chemically interpreted as the result of alternating depletion and enrichment in trace elements (i.e., U and Y) during crystal growth at the crystal-melt interface. On the other hand, bright domains at the rim of the crystals, usually show higher and homogeneous CL-emission, presenting irregular shape that often overgrowths the pre-existing crystals whose texture can sometime be evident as a relict (Rubatto and Gebauer 2000). Our results indicate that metamorphic zircon can also include dark, bright, oscillatory, and undefined CC domains. This observation warns against using internal zircon zoning alone for distinguishing magmatic vs. metamorphic zircon.

As portrayed in Figure 2, bright and oscillatory-bright domains (yellow and blue spots, respectively), have the same Raman spectral features and define the “non-metamict” domain in the diagram (see details below). Therefore, despite their magmatic or metamorphic orgin, from a structural point of view, these two domains are effectively equivalent and there will be no difference in the calculated residual stress even if the two domains coexist within the same crystal.

Finally, if we plot our data in a Γ of the B_1g mode near 1008 cm⁻¹ vs. Dα diagram, most of the Alpine domains fall within the broad interpolation line given by Nasdala et al. 2001 (Fig. 6). Nasdala et al. (2001) interpreted the data plotting outside of such interpolation band as the possible effects of thermal annealing in the crystal. In our case, we interpret the few points outside the interpolation band as the result of possible partial annealing of inherited zircon cores (Permian and Caledonian in age) during the Alpine metamorphism.

Previous studies have established that partial metamictization can have a major influence on zircon Raman spectra (e.g., Zhang et al. 2000; Geisler et al. 2001; Nasdala et al. 2001; Binvignat et al. 2018). In general, depending on the degree of metamictization, zircon crystals can show a significant Raman peak broadening and frequency shift toward lower wavenumbers with respect to pristine samples taken as reference. Furthermore, Geisler et al. (2001) pointed out the possible heterogeneous (step-evolution) effect that annealing could have on peak broadening and position.

Based on these considerations, the relationship between the peak width and the Raman shift (Fig. 2) can be exploited as a discriminant between partially metamict/annealed crystals and pristine crystals, even when they are completely buried within their mineral host (see also Zhong et al. 2019). Indeed, differently to phonon wavenumber, variation in the Γ of a Raman peak is independent from minor amounts of stress (Binvignat et al. 2018) as those usually recorded in host-inclusion systems. It follows that, when considering zircon inclusions, a buried partially metamict crystal should present Γ comparable to partially metamict crystals exposed at the surface and the same is valid for well-crystalline zircon grains. Therefore, we conclude that measurements with Γ greater than 5 cm⁻¹ showing an inverse relationship with the Raman shift, represent partially metamict/annealed domains. On the other hand, domains with lower Γ (i.e., 3.5–5.0 cm⁻¹) and with no correlation with the Raman shift are associated to a non-metamict domain. In this regard, it is worth noting that also Alpine zircons with high U contents fall into the metamict domain.

Geisler and Pidgeon (2002) pointed out that possible annealing processes may influence both the Raman shift and the peak width, as well as the relationship between them during secondary geological processes. However, from their results, it is evident that such effects are critical for zircon with moderate to heavy levels of radiation damage (i.e., ω < 1004 cm⁻¹ and Γ > than 11–12 cm⁻¹). In our case, the zircon inclusion displaying a negative slope result all with Γ < 11 cm⁻¹ and ω > 1004 cm⁻¹. This correlation can be tentatively interpreted as an indication that, as previously stated (e.g., Fig. 6), annealing effects were negligible or absent in most of the samples over the relatively short geological evolution (280 or 35 Ma to present) and fast subduction metamorphism (Gebauer et al. 1997; Rubatto and Hermann 2001; Gauthiez-Putallaz et al. 2016).

However, it is difficult to make a rigorous prediction of the effect of metamictization on the determination of the residual pressure of a buried inclusion. For zircon with Γ larger than 5 cm⁻¹ there are at least two additional unknown variables:

• The reference wavenumber ω₀ of an equally metamict free crystal;

• The phonon mode compressibility β = [1/ω₀][dω/dP].

Both variables are very sensitive to the accumulated radiation dose D_α (e.g., Binvignat et al. 2018). Unfortunately, the dispersion of our data in the established Γ–D_α trend is too large to give a reliable value of D_α and consequently of ω₀ and β for the purposes of Raman thermobarometry.

For completely unstressed pristine crystals, the phonon wavenumber usually is expected not to exceed the values of 1008.5–1009 cm⁻¹ (see Geisler et al. 2001; Binvignat et al. 2018). Nevertheless, our partially exposed grains, even considering only the non-metamict domain in Figure 2, show a larger variation (up to 1011 cm⁻¹).

As reported by Campomenosi et al. (2018), partially exposed inclusions can still preserve a notable stress state in terms of Raman peak shift as a function of the inclusion exposition degree and size. To better clarify this issue, Figure 7 shows a Γ vs. Δω diagram of partially exposed inclusions discriminating between grains with different sizes (<50 and >50 μm). Note that, in this case, we show only data from the exposed inclusions showing bright and or oscillatory-bright domains at CC because only for such domains can we safely exclude other effects described above. Data from inclusions of different size largely overlap. However, about 20% of the smaller inclusions have a higher Δω than the larger one (with one exclusion) and reach a Δω of 572 cm⁻¹ that are never observed in the larger inclusions. This is in agreement with theoretical predictions (Mazzucchelli et al. 2018) for which, within the same host, small inclusions, when partially exposed, tend to retain more stress, due to the larger surface-to-volume ratio, i.e., larger impact of the host upon relatively larger surface of the inclusion. On the other hand, small inclusions should exhibit less dispersion of the degree of exposure. Hence, the quite large spread of our data is most probably due to the superposition of both effects related to the inclusion size: the surface-to-volume ratio and the degree of exposure.

Based on the above considerations, we propose a simple protocol for selecting zircon inclusions to get reliable residual pressure estimates for elastic thermobarometry.

• Buried crystals must be isolated from other inclusions, section surfaces, host boundary, cracks, and from any other kind of boundaries. This implies that thicker sections (i.e., 300 μm) are better (see also Campomenosi et al. 2018).

• At least 2–3 measurement spots for each crystal moving from core to rim should be performed: this enables detection of eventual structural heterogeneity within the crystal. For this pourpose, it is better to use a spectrometer with a confocal optical microscope to reach the best spatial resolution.

• An instrumental sepctral resultion equal or better than 2 cm⁻¹ should be used and only those inclusions presenting a Γ of the B_1g mode near 1008 cm⁻¹ smaller than 5 cm⁻¹ should be considered. Note that this value represents the summation of the analytical uncertainty (i.e., instrumental spectral resolution) to the minimum value of about 3 cm⁻¹ that has been measured on several grains by the two different spectrometers (see Fig. 4).

• Only inclusions in which the phonon wavenumber and the Γ is constant across the entire crystal volume should be selected.

• If Γ is constant, but the phonon wavenumber changes across the crystal volume there could be effects due to the shape of the inclusion (Campomenosi et al. 2018). In this case it is reccomended to consider the measurement collected at the center of the crystal to calculate the residual pressure of the inclusion.

• Partially exposed crystals should not be used as reference to calculate residual stress in the buried inclusions: they can be still under a notable residual stress state. It is reccomended using a completely free crystal or large exposed inclusions for which the residual stress state and the metamictization effects are negligible (see for example Fig. 1d).

• Whenever possible, a statistically significant amount of partially exposed inclusions should be selected to double check their textural complexity by CC or CL imaging and their chemistry. As an alternative, imaging and chemical checks should be performed on inclusions that have been exposed after the Raman measurements.

• Chemical and age measurements should be considered as important corollary information to reconstruct the petro-genesis of the zircon inclusions and the garnet host.

Zircon inclusions are difficult to manage correctly for elastic thermobarometry applications, and a detailed characterization of the inclusions should be performed before extracting thermobarometric data. In our systematic study, the combination of structural and chemical information obtained by different analytical techniques on partially exposed zircon inclusions allows to reliably predict the structural state of buried crystals for Raman spectroscopic measurements.

Our results provide a solid basis for the selection of reliable zircon inclusions to use for elastic thermobarometry applications. This method, in combination with the already rich tool set that can be applied to zircon [e.g., Ti-in-zircon thermometry (Watson et al. 2006), U-Pb geochronology, oxygen isotopes], will provide an even more detailed characterization of P-T-t fluid and deformation hystory of metamorphic rocks.

This work was financially supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program grant agreements 714936 to M. Alvaro. N. Campomenosi acknowledges the University of Genova for funding. D. Rubatto acknowledges the financial support of the Swiss National Science Foundation (project 166280).

Thanks are due to F. Piccoli and to A. Berger (University of Bern) for assistance in the LA-ICP-MS laboratory and the SEM laboratory, respectively, at the Institute of Earth Science of Bern and to N. Waeselmann for assistance in the laboratory of Raman spectroscopy at the University of Hamburg. We are grateful to D. Gatta for manuscript editing and to M. Cisneros and an anonymous reviewer for constructive comments.