From long-lived batholith construction to giant porphyry copper deposit formation: petrological and zircon chemical evolution of the Quellaveco District, Southern Peru

Porphyry Cu ore deposits are a rare product of arc magmatism that often form spatiotemporal clusters in magmatic arcs. The petrogenetic evolution of igneous rocks that cover the temporal window prior to and during porphyry Cu deposit formation may provide critical insights into magmatic processes that are key in generating these systems. This study documents the magmatic evolution of the Palaeocene–Eocene Yarabamba Batholith, Southern Peru, that was incrementally assembled between ~ 67 and ~ 59 Ma and hosts three, nearly contemporaneous, giant porphyry Cu–Mo deposits that formed at 57–54 Ma (Quellaveco, Toquepala and Cuajone). Whole-rock geochemistry, U–Pb geochronology and zircon trace element chemistry are reported from Yarabamba rocks that span the duration of plutonic activity, and from six porphyry intrusions at Quellaveco that bracket mineralisation. A change in whole-rock chemistry in Yarabamba intrusive rocks to high Sr/Y, high La/Yb and high Eu/Eu* is observed at ~ 60 Ma which is broadly coincident with a change in vector of the converging Nazca plate and the onset of regional compression and crustal thickening during the first stage of the Incaic orogeny. The geochemical changes are interpreted to reflect a deepening of the locus of lower crustal magma evolution in which amphibole ± garnet are stabilised as early and abundant fractionating phases and plagioclase is suppressed. Zircons in these rocks show a marked change towards higher Eu/Eu* (> 0.3) and lower Ti (< 9 ppm) compositions after ~ 60 Ma. Numerical modelling of melt Eu systematics and zircon-melt partitioning indicates that the time series of zircon Eu/Eu* in these rocks can be explained by a transition from shallower, plagioclase-dominated fractionation to high-pressure amphibole-dominated fractionation at deep crustal levels from ~ 60 Ma. Our modelling suggests that any redox effects on zircon Eu/Eu* are subordinate compared to changes in melt composition controlled by the fractionating mineral assemblage. We suggest that growth and intermittent recharge of the lower crustal magma reservoir from ~ 60 Ma produced a significant volume of hydrous and metallogenically fertile residual melt which ascended to the upper crust and eventually generated the three giant porphyry Cu–Mo deposits at Quellaveco, Toquepala and Cuajone from ~ 57 Ma. Our study highlights the importance of high-pressure magma differentiation fostered by strongly compressive tectonic regimes in generating world-class porphyry Cu deposits.


Introduction
Porphyry Cu ore deposits are the products of large, longlived, trans-crustal arc magma systems that release and focus metal-charged fluids during the intrusion of porphyritic stocks and dykes (Dilles 1987;Seedorf et al. 2005;Sillitoe 2010). The occurrence of these economically valuable deposits is typically restricted to narrow arc segments-in lineaments or clusters that are formed during specific temporal windows within arc evolution (Sillitoe 2010). The formation of porphyry Cu deposits has been linked to distinctive changes in the trace element chemistry of igneous rocks that crystallised immediately prior to ore genesis, including increased whole-rock Sr/Y (> 50), La/Yb (> 20) and Eu/Eu* (> 0.3) (Lang and Titley 1998;Rohrlach and Loucks 2005;Richards and Kerrich 2007;Richards 2011;Richards et al. 2012;Loucks 2014).
The origin of the characteristic trace element signatures of magmas parental to porphyry ore deposits has been attributed to a protracted evolution in the lower crust (e.g., Rohrlach and Loucks 2005;Chiaradia et al. 2009;Wilkinson 2013;Chelle-Michou et al. 2015), linked to an inhibition of magma ascent. This is thought to be related to perturbations in regional geodynamics that generate intense but transient compression and crustal thickening (Cooke et al. 2005;Rosenbaum et al. 2005). Longterm evolution in the lower crust causes mantle-derived magmas to undergo extensive fractional crystallisation and assimilation of pre-existing crust, thereby attaining diverse, evolved compositions (DePaolo 1981;Hildreth and Moorbath 1988;Lang and Titley 1998;Annen et al. 2006). Here, high pressure (> 0.7 GPa) combined with high melt water contents (> 5 wt.%) cause a change in the fractionating mineral assemblage and its crystallisation sequence (Loucks 2014;Müntener and Ulmer 2018). This increased water content and pressure promotes the stability of amphibole ± garnet (in which Y and MREE/ HREEs are compatible) as early and abundant crystallising and fractionating phases, whereas crystallisation of plagioclase (in which Sr and Eu are compatible) is suppressed (Müntener et al. 2001;Alonso-Perez et al. 2008;Melekhova et al. 2015).
The reason for the association of giant porphyry Cu ore deposits with high Sr/Y and high La/Yb ratios in coeval igneous whole rocks remains equivocal. However, extensive evolution of magmas in the lower crust is thought to enrich them in water (e.g., Melekhova et al. 2015), in part because water is more soluble in melts at high pressure (e.g., Newman and Lowenstern 2002). Modelling suggests that following ascent to the shallow crust, such hydrous (> 5 wt.% H 2 O) magmas are predisposed to exsolve enough metalliferous fluid during decompression and crystallisation to generate and maintain large magmatichydrothermal ore systems (Chiaradia and Carrichi 2017). Furthermore, a higher melt fO 2 obtained during lower crustal fractionation (Lee et al. 2005;Tang et al. 2018;Ulmer et al. 2018) has been suggested to increase sulphur and metal solubility, resulting in delayed sulphide saturation that could otherwise lead to loss of chalcophile elements from the melt (Jenner et al. 2010;Lee et al. 2012). Recent studies have shown that sulphide saturation in arc magmas is inevitable even under oxidised (> NNO + 2) and lower crustal conditions (Matjuschkin et al. 2016;Du and Audétat 2020), indicating that exceptionally high melt fO 2 may not be fundamental in forming porphyry Cu deposits.
Further studies of the petrogenesis of magmas parental to porphyry Cu deposits have focused on zircon trace element chemistry, because zircon is resistant to re-equilibration during the extensive hydrothermal alteration associated with these systems. A key observation has been the presence of higher Eu/Eu* in zircons from igneous rocks spatiotemporally associated with porphyry Cu systems, compared to older intrusions not associated with mineralisation in the same igneous complex. This has been used to argue for an elevated oxidation state (fO 2 ) of magmas parental to porphyry Cu deposits (e.g., Ballard et al. 2002;Liang et al. 2006;Dilles et al. 2015), because Eu is more compatible in zircon under oxidised conditions, where Eu 3+ is dominant over Eu 2+ (Burnham and Berry 2012). However, this remains ambiguous, since zircon Eu/Eu* is sensitive to several other parameters such as the co-crystallising mineral assemblage, melt composition and temperature (Burnham et al. 2015;Loader et al. 2017).
In this study, we aim to better understand the magmatic processes that determine the formation of giant porphyry Cu ore deposits by tracking the long-lived magmatic evolution of the Yarabamba Batholith and its associated porphyry Cu-Mo deposits, Southern Peru. The compositional evolution of the batholith and the igneous rocks associated with porphyry mineralisation provide constraints on the tempo and nature of the magmatic evolution that preceded and ultimately generated three giant porphyry Cu-Mo deposits in the district: Quellaveco, Cuajone and Toquepala. Our approach integrates whole-rock chemistry with zircon U-Pb geochronology and trace-element chemistry to show that both whole-rock and zircon data record a transition in magma chemistry at ~ 60 Ma that reflects a deepening of the locus of crustal magma evolution, prior to district-wide mineralisation.

Regional geology
The Yarabamba Batholith is part of the Coastal Batholith which extends 1100 km along much of the length of the Peruvian Andes (Myers 1975;Cobbing and Pitcher 1972;Clark et al. 1990). The Yarabamba Batholith intrudes the Toquepala Group, which predominantly comprises intermediate-felsic volcanics of mid to late Cretaceous age (91-69 Ma;Simmons 2013). The Toquepala volcanics formed during a period of major tectonic perturbation that led to the closure of back-arc basins, amalgamation of arc terranes and deformation of sedimentary and volcanic rocks (Mpodozis and Ramos 1989). Segmentation of the arc occurred due to lateral heterogeneity in the subducting Nazca plate, creating changes in the angle of subduction and, therefore, the style and intensity of compression and volcanism (Ramos 2009;Capitanio et al. 2011;Hu et al. 2016). The youngest rocks of the Toquepala Group (~ 67 Ma), typically found at high altitudes between Quellaveco and Toquepala, overlap in age with the oldest intrusive activity of the Yarabamba Batholith (67-69 Ma) (Simmons 2013).
The Yarabamba Batholith is a composite plutonic complex ( Fig. 1) assembled by the emplacement of discrete batches of dioritic-granodioritic magma between ~ 67 and 59 Ma (Mukasa 1986;Demouy et al. 2012;Simmons et al. 2013) during a period of relatively oblique and low rates of convergence of the Nazca plate (Pardo-Casas and Molnar 1987;Mpdozis and Cornejo 2012). Lead isotopic data indicate that these magmas experienced little interaction with the Precambrian granulitic basement (Barreiro and Clark 1984). The final phase of this magmatic activity was the emplacement of three, near contemporaneous, major porphyry Cu-Mo systems in the Eocene (Quellaveco 58.4-54.3 Ma;Toquepala 57.0-54.0 Ma;and Cuajone 56.5-53.0 Ma;Sillitoe and Mortensen 2010;Simmons et al. 2013). The emplacement of these porphyry Cu-Mo systems has been associated with high rates of tectonic compression (e.g., Benavides-Cáceres 1999). Post-mineralisation, Fig. 1 Regional geological map of the Yarabamba Batholith and giant porphyry Cu-Mo deposits. Colourless areas are Palaeocene volcanics that pre-date the batholith, younger Miocene volcaniclastic and sedimentary rocks, and Plio-Pleistocene sediments. The Quellaveco deposit is hosted within the Yarabamba Batholith. Shapes show the location of regional samples of the Yarabamba Batholith with ages shown for samples that were dated, all other samples were collected from drill core at Quellaveco. Symbol shapes and colours are those used in subsequent figures (see key Fig. 2). Simplified after Bellido (1979). Inset map of Peru shows the location of the study area (star).
†Ages published in Simmons et al. (2013) and whole-rock geochemical data for these samples are reported herein 12 Page 4 of 21 the region underwent up to 40° counter-clockwise rotation during later crustal shortening and oroclinal bending during the Eocene-Oligocene (Arriagada et al. 2008).

The Quellaveco porphyry Cu-Mo deposit
The Quellaveco porphyry Cu-Mo deposit (Fig. 2) is one of the world's largest unexploited copper resources with ore reserves of 1.3 billion tonnes at 0.57% Cu (Anglo American, Annual Report 2019). At least six intrusive porphyry phases have been identified at the deposit which form stock and dyke complexes (Fig. 2) emplaced over ~ 4 Myr (Sillitoe and Mortensen 2010). The Quellaveco porphyries are hosted within a pluton that is a late member of the Yarabamba Batholith (LA-ICP-MS U-Pb zircon date of 59.46 ± 0.24 Ma; Sillitoe and Mortensen 2010), referred to as the Quellaveco Granodiorite. This host rock is an equigranular granodiorite, with alteration grading from weak propylitic (chlorite + epidote) at the periphery of the system, to strong potassic (biotite + K-feldspar) at the core.
Following the lithological classification of Simmons (2013), drill core observations of cross-cutting relationships allow the relative chronology of the porphyry phases to be reconstructed, especially where intrusive contacts truncate vein generations and alteration assemblages (Fig. 3). All units at Quellaveco exhibit a porphyritic texture dominated by feldspar phenocrysts; the emplacement of each of these was accompanied by separate phases of sulphide-mineralisation and hydrothermal alteration that generally decrease in intensity through the lifespan of the system (Simmons 2013). Each phase of hydrothermal alteration consisted of an initial potassic alteration that replaced mafic minerals and felspars with secondary biotite and K-feldspar, respectively. This was followed by texturally destructive phyllic alteration associated with quartz, sericite and pyrite deposition.
The earliest porphyry unit observed at Quellaveco is termed the "Granodiorite Porphyry". Previously only found as xenoliths within the later porphyry units, it has been intersected at depth in recent drilling. The Granodiorite Porphyry is characterised by a more equigranular texture compared to the other porphyry units (euhedral > 2 mm amphibole phenocrysts) and exhibits intense potassic alteration and extensive chalcopyrite mineralisation. Cross-cutting relationships indicate that its emplacement was associated with weak sulphide-mineralisation and preceded the bulk of Cumineralisation at Quellaveco (Simmons 2013). The second unit emplaced at Quellaveco, the "Early Porphyry", forms a stock ~ 1 × 2 km in surface area and hosts most of the Cu-ore in the deposit. This unit is granodioritic and is characterised by plagioclase (25%), K-feldspar (15%), euhedral quartz (15%) and biotite phenocrysts (10%). The Early Porphyry is cut by a sequence of porphyry dykes and a stock complex in the centre of the deposit (Fig. 2): (1) The first inter-mineral porphyry, the Monzonite Porphyry, occurs as a series of dykes and is characterised by its high phenocryst abundance (60%) including feldspars (50%), anhedral quartz (10%) and biotite (5%).

Sample selection
To represent the time period of magmatic evolution prior to and during porphyry Cu ore deposit formation, samples were collected from Yarabamba intrusions in the district that predate mineralisation and from intrusions at the Quellaveco porphyry Cu-Mo deposit. Thirteen samples of Yarabamba intrusive rock were collected from several outcrops found in the district (Figs. 1, 2) for whole-rock geochemical analysis, including the Yarabamba host rocks for each of the three major porphyry Cu-Mo deposits: (1) six from a granodiorite intrusion that hosts and surrounds the Toquepala deposit, including one sample from the Toquepala open pit; (2) two from the diorite-monzonite intrusions that host, and crop out to the west of, the Cuajone deposit; (3) four from an extensive diorite-monzonite intrusion at the southern limit of the district; (4) one from the Quellaveco Granodiorite, sampled to the east of the Quellaveco deposit. Samples of least-altered Quellaveco porphyry intrusions (n = 51) and a further seven samples of the Quellaveco Granodiorite were collected from logged drill core at the Quellaveco deposit for whole-rock geochemical analysis. From this sample-set, a representative sub-set of fifteen duplicate samples were chosen for zircon separation. Six samples from the Yarabamba intrusive suite were selected plus nine drill core samples from Quellaveco: one sample of mineralised host rock (Quellaveco Granodiorite), one sample of the earliest Granodiorite Porphyry, one sample of the Early Porphyry, two samples of inter-mineral Monzonite Porphyry, one sample of late inter-mineral Monzodiorite, two samples of the Late Porphyry and one sample of a postmineralisation dacite dyke.

Whole-rock major and trace element analysis
Crushing and whole-rock geochemical analysis were carried out at Bureau Veritas (Vancouver, Canada) and the codes for  each package are provided below. Further details are available at http://acmel ab.com/servi ces/metho d-descr iptio ns/. Samples were crushed to ≥ 70% passing 2 mm, and a 250 g split of this material was then pulverised to ≥ 85% passing 75 μm. The remainder of the initial coarse product was used for mineral separation. Pulps were then dried at 105 ºC prior to analysis, with fractions taken from the pulp for different analytical procedures.
A 0.2 g aliquot of each sample was fused using lithium borate flux to make a glass bead. This was subsequently dissolved in 5% nitric acid. Major elements (SiO 2 , Al 2 O 3 , Fe 2 O 3 , MgO, CaO, Na 2 O, K 2 O, TiO 2 , P 2 O 5 , MnO, Cr 2 O 5 , Ba and Sc) were determined by inductively-coupled plasma atomic emission spectrometry (ICP-AES) under procedure LF300. Refractory and rare earth elements were determined using the same solution using inductively-coupled plasma mass spectrometry (ICP-MS) under procedure LF100. Elements measured were Ba, Be, Ce, Co, Cs, Dy, Er, Eu, Ga, Gd, Hf, Ho, La, Lu, Nb, Nd, Pr, Rb, Sm, Sn, Sr, Ta, Tb, Th, Tm, U, V, W, Y, Yb and Zr. Loss on ignition (procedure TG001) was determined by heating 1 g of sample to 1000 ºC and measuring the mass difference. Carbon and sulphur (procedure TC000) were determined by combusting 0.1 g of sample using LECO (infrared absorption and thermal conductivity). Information on the standards analysed and the accuracy and precision of whole-rock analysis is available in Supplementary Material 2.

Zircon geochronology and trace element analysis
Coarse rejects for selected samples were sieved to 500 μm, panned for heavy mineral separates, and passed under a neodymium magnet to remove magnetic grains. Zircons from each heavy mineral separate were picked, mounted in epoxy and polished to expose crystal cores in preparation for microanalysis.
To characterise zircon textures, each population was imaged by scanning electron microscope-cathodoluminescence (SEM-CL) using a Zeiss EVO SEM located in the Imaging and Analysis Centre (IAC), Natural History Museum, London. The instrument was operated with a 3.0 nA beam current and 10 kV accelerating voltage. Subsequent zircon trace element analysis was conducted in the LODE Laboratory, IAC, using an Agilent 7700x quadrupole ICP-MS coupled to an ESI New Wave Research NWR193 excimer laser. The laser was operated using a 5 Hz repetition rate, a fluence of ~ 3.5 J cm −2 and a spot size of 30 μm. Spots were selected using SEM-CL images to allow targeting of individual crystal domains for intra-crystal, core-rim analyses. The following isotopes were analysed: 27 Al, 28 Pb,207 Pb and 238 U (30 ms); and 27 Al, 43 Ca, 57 Fe and 91 Zr (5 ms). NIST-610 was the primary reference material used for trace element quantification and NIST612 and zircon 91,500 were used to monitor accuracy and precision. A stoichiometric value for Si (15.3 wt.%) was used as the internal standard. GJ-1 (ca. 601 Ma; Jackson et al. 2004) was used as the primary geochronology standard and Fish Canyon Tuff and zircon 91,500 were analysed between each set of unknowns to ensure consistency and accuracy in the U-Pb ages obtained (Wiedenbeck et al. 2004;Wotzlaw et al. 2013). Data reduction was performed using Iolite v. 3.61 (Paton et al. 2011), where the 206 Pb/ 238 U and 207 Pb/ 235 U spectra were monitored for isotopic heterogeneity. Analyses that were found to be discordant or with anomalously high Al, P, Ca, Fe and La (indicating the presence of mineral and melt inclusions) are not included in the reported data. Further information on the analytical methods, standards, typical analytical precision and limits of detection is available in Supplementary Material 2.
Three samples were dated by secondary ion mass spectrometry (SIMS) using the Sensitive High Resolution Ion Microprobe-Reverse Geometry facility (SHRIMP-RG) at Stanford University, in co-operation with the US Geological Survey. Zircons were sputtered with an O 2-primary ion beam varying from 4-5 nA with a spot diameter of 20-40 μm and a depth of 1-2 μm. The acquisition routine began with the high mass normalizing species ( 90 Zr 2 16 O + ), followed by 204 Pb + , a background measured at 0.050 mass units above 204 Pb + , 206 Pb + , 207 Pb + , 208 Pb + , 238 U + , 232 Th 16 O + and 238 U 16 O + and 232 Th + . The obtained dates were standardised against the R33 zircon standard (419 Ma; Black et al. 2004). A full description of the analytical technique and data reduction is provided in Simmons et al. (2013) and in Supplementary Material 1.

U-Pb geochronology
A total of 762 zircon analyses are reported from twelve samples with both U-Pb dates and trace element chemistry, acquired using LA-ICP-MS (Supplementary Material 2). A further 42 zircon analyses are reported with only U-Pb dates, determined by the SHRIMP-RG method. For the individual samples, weighted mean dates were calculated for concordant 206 Pb/ 238 U dates together with mean square weighted deviation (MSWD) values, which ranged from 0.21 to 5.31 ( Fig. 4   Pb/ 238 U dates of ± 2%. Annotated dates are given with two uncertainties: (1) the 2σ uncertainty from mean square weighted deviation (MSWD) and (2) with propagated external uncertainty of 2%. MSWD and number of analyses are also given. Samples are plotted left-toright according to age ranking for the Yarabamba rocks and plotted left-right according to cross-cutting relationships for the Quellaveco porphyry rocks. Symbols besides sample units show the symbology used in zircon trace-element plots; for the three samples with no symbol shown, these samples were dated using the SHRIMP-RG method and zircon trace elements are not reported here. Black vertical bars indicate the range of U-Pb zircon LA-ICP-MS dates reported by Simmons et al. (2013) for the neighbouring Cuajone and Toquepala porphyry Cu systems. †One Quellaveco Granodiorite sample yielded few concordant ages, so a discordia date was calculated and reported instead of a weighted mean and Table 1). For both analytical techniques, weighted mean dates are reported here with two uncertainties, the first calculated as the two-sigma error calculated from the weighted mean and the second calculated with a propagated external uncertainty of 2% (Horstwood et al. 2016). A 2% external uncertainty was selected based on long-term reproducibility of reference materials across different positions in the sample cell for the LA-ICP-MS system used. The propagation of this added external uncertainty allows systematic bias to be accounted for and allows comparison of dates determined by the two analytical techniques. The single sample from the monzonite pluton southwest of Toquepala (sample QVC075) gave a weighted mean date of 67.2 ± 0.79/1.56 Ma (n = 14, MSWD = 0.79). The two samples from the granodiorite pluton west of Toquepala (samples QVC082, QVC083 and Qu380) yielded weighted means of 60.96 ± 0.66/1.39 Ma (n = 23, MSWD = 1.82), 61.64 ± 0.26/1.26 Ma (n = 95, MSWD = 2.63) and 60.28 ± 1.00/1.57 (n = 16, MSWD = 5.31), overlapping within the reported standard error. An aplite dyke (Qu381) cutting the Toquepala Granodiorite gave a similar date to its host rock at 60.18 ± 0.67/1.38 Ma (n = 14, MSWD = 2.37). The Quellaveco Granodiorite (samples QVC028 and Qu008) yielded a weighted mean date of 58.99 ± 0.23/1.20 Ma (n = 94, MSWD = 5.41) and a discordia date of 59.99 ± 0.89/1.49 (n = 12, MSWD = 0.81).

Zircon trace element chemistry
Zircon trace element data show systematic trends with positive correlations between Hf and Yb/Gd (Fig. 8a), and between Th/U and Ti (Fig. 8b). However, in terms of several parameters (particularly Ti and Eu/Eu*), the zircons can be separated into two, chemically-distinct groups (Fig. 8c, d). High Ti (5-20 ppm) and low Eu/Eu* (0.10-0.35) characterise those from the regional Yarabamba rocks, and low Ti (< 9 ppm) and high Eu/Eu* (0.3-0.7) typify zircons from the Quellaveco Granodiorite and porphyry units. Zircons from the eight Quellaveco porphyry intrusions cannot be distinguished; those from the Quellaveco Granodiorite overlap with the porphyries but extend to lower and more restricted Hf (9000-11,000 ppm) and Ti contents (3-8 ppm), and higher and more restricted Th/U ratios (0.5-1.3).
Core and rim analyses were performed on a sub-set of zircon crystals that showed variable intra-grain trends (Fig. 9). In the regional Yarabamba rocks, zircon Eu/Eu* (Eu/Eu* zircon ) remains low (< 0.4) and there is a general increase in Hf from core to rim (Fig. 9a). In the Quellaveco Granodiorite, there is a lot of variability in core-rim Eu/ Eu* and Hf trends (Fig. 9b). Conversely, Eu/Eu* zircon in the Quellaveco porphyry intrusions more consistently decreases from core to rim, coupled with an increase in Hf (Fig. 9c).

Timescale of assembly of the Yarabamba Batholith and emplacement of the Quellaveco porphyries
Many of the reported weighted mean 206 Pb/ 238 U zircon dates return MSWDs that are greater than the acceptable MSWD expected for the number of data points (Wendt and Carl 1991). This either indicates that the overdispersion is the result of real geological scatter or an underestimation of the uncertainties of each spot analysis. While it is difficult to exclude real geological scatter as the main cause of the overdispersion, the similar MSWDs and the normal distribution of the data for each sample suggests that under-propagated uncertainties on each spot may play a more significant role, rather than continuous zircon crystallisation over extended durations (e.g., Schoene 2014; Large et al. 2020). We provide recalculated uncertainties of each grain to fit the criteria of a single zircon population (see Supplementary Material 1). However, recalculating the weighted means and standard errors with theoretically higher uncertainties on each date does not change the main findings of the study.
The reported weighted mean dates for the Yarabamba intrusive phases and the Quellaveco porphyry units represent the average age of zircon crystallisation but also constrain the age of emplacement, since the uncertainties on the weighted mean of each intrusive phase are greater than the typical hundred kyr duration of zircon crystallisation in mid-upper crustal magma reservoirs (e.g., Schoene et al. 2012;Wotzlaw et al. 2013;Samperton et al. 2015;Large et al. 2020). We herein refer to the weighted mean dates as the age of emplacement.
The emplacement ages determined for the thirteen samples from the Yarabamba Batholith (Fig. 4), in combination with their spatial separation, cross-cutting relationships and previously determined U-Pb ages (Simmons 2013), The emplacement ages of the Quellaveco porphyry units indicate that all the porphyries were emplaced within 2.39 ± 1.61 Myr. Considering the ages of the Granodiorite Porphyry and the post-mineralisation dacite dykes that bracket mineralisation, ore formation at Quellaveco is constrained to between 57.91 ± 1.33 Ma and 55.52 ± 1.16 Ma. However, the bulk of sulphide-mineralisation is thought to have been associated with the emplacement of the early and intermineral porphyry units (Simmons 2013), suggesting that mineralisation occurred over a timescale smaller than the analytical error for the U-Pb methodology used here. The emplacement ages reported here are consistent within uncertainty with previous work on the deposit (Sillitoe and Mortensen 2010;Simmons 2013), with the exception of the "Early Porphyry" we sampled which is considerably younger (by ~ 4 Myr) than the LA-ICP-MS U-Pb date of 58.41 ± 0.53 Ma for an "Early Porphyry" reported by Sillitoe and Mortensen (2010). However, our results are consistent with data reported for the "Early Porphyry" by Simmons (2013), suggesting that Sillitoe and Mortensen (2010) sampled an older intrusive phase.
The emplacement ages determined for the Quellaveco porphyries overlap within uncertainty with those for the neighbouring deposits of Cuajone (56.5-53.0 Ma) and Toquepala (57.0-54.0 Ma) ( Fig. 4; Simmons et al. 2013). The resolution of the current dates is insufficient to establish the degree of contemporaneity of the three mineralisation centres (e.g., Large et al. 2020); however, the new ages suggest that the magmatic processes that were conducive to porphyry Cu-Mo deposit formation occurred broadly synchronously on a district-scale, in at least these three  mineralisation centres. Considering the distance between these deposits (10-20 km; Fig. 1), and the size of exposed, inferred source intrusions for porphyry deposits elsewhere (e.g., Yerington and Bingham Canyon: Dilles et al. 2000;Steinberger et al. 2013), it is speculated that the magma and fluids that formed these three systems were sourced from a common, upper crustal magma reservoir.

Long-lived magmatic evolution culminating in porphyry Cu-Mo mineralisation
The distribution of ages obtained for the samples from the Yarabamba Batholith indicate at least four phases of intrusive activity prior to mineralisation between ~ 67 and 57 Ma (Fig. 4). The identification of these four phases does not preclude the existence of further magmatic activity that was not sampled; however, these samples are considered to be representative of several stages of the extended time period of batholith construction. The chemical compositions of these samples provide critical petrogenetic information on the magmatic evolution that led up to the formation of the three giant porphyry Cu-Mo deposits at ~ 55 Ma. In addition, the porphyry intrusive phases at Quellaveco provide a series of aliquots of the underlying magma system that released magmatic-hydrothermal fluid resulting in mineralisation. The observed differences in whole-rock compositions between the Yarabamba and Quellaveco porphyry suites reflect a distinct change in magma chemistry at ~ 60 Ma, characterised by a transition to higher Sr/Y, La/Yb, Eu/Eu* and Gd/Yb (Fig. 7), as well as generally more evolved major element compositions (e.g., SiO 2 > 65 wt.%). These differences could reflect changes in intra-crustal magma evolution processes and/or changes in the mantle source composition. The relatively evolved and limited range in SiO 2 and MgO concentrations of the studied samples does not permit a rigorous assessment of whether the recorded geochemical changes are linked to variable primary melt compositions and mantle source processes. However, the concomitant change in SiO 2 content with the shifts in Sr/Y, La/Yb and Eu/Eu* provides evidence that such signatures more likely reflect intra-crustal evolution. Furthermore, in modern arcs, there is a broad correlation between Sr/Y and La/Yb with  Thus, while a change in mantle source cannot be excluded as a control on the shift in magma geochemistry observed in this study, we suggest intra-crustal processes primarily drove this fundamental change.
Experimental studies (e.g., Müntener and Ulmer 2018;Nandedkar et al. 2014) have demonstrated that in hydrous melts (> 4 wt.% H 2 O), at pressures found at depth in thickened arc crust (> 0.7 GPa), amphibole (in which Y, MREEs and HREEs are compatible: Nandedkar et al. 2016) is stabilised as a liquidus phase at the expense of plagioclase (in which Sr and Eu are compatible; Aigner-Torres et al. 2007) and clinopyroxene (in which Y, MREEs and HREEs are weakly compatible; Luhr et al. 1984a, b). This fractionation of an amphibole-dominant assemblage from a primitive arc magma can account for the high Sr/Y and La/Yb values observed in global porphyry Cu datasets (Richards and Kerrich 2007), and suppression of plagioclase would further increase Sr/Y and increase Eu/Eu*. However, the increases in Sr/Y and La/Yb observed here are accompanied by an increase in Gd/Yb (Fig. 7), which is not characteristic of an amphibole-dominated fractionation pathway, because amphibole preferentially partitions Gd over Yb (e.g., Nandedkar et al. 2016). This may be accounted for by the additional presence of garnet in the fractionating assemblage (or as a residual phase in equilibrium with the magma), which has partition coefficients an order of magnitude higher for HREEs than for MREEs (e.g., Green et al. 2000). Garnet has been stabilised in experimental studies at either high pressures (> 0.8 GPa) or at high melt water contents (> 8 wt.%) (Müntener et al. 2001;Alonso Perez et al. 2009). The traceelement geochemical evolution, therefore, suggests a change in the fractionating assemblage in the lower crust, as a result of increased pressure and/or melt water contents, stabilising amphibole ± garnet at the expense of plagioclase.
In porphyry-related rocks, trace element signatures such as high Sr/Y have been shown to become more pronounced during magmatic differentiation (Loucks 2014). Thus, it is possible to suggest that the absence of these signatures in the older, Yarabamba rocks is because they had undergone less magmatic differentiation (< 65 wt.% SiO 2 ) than the younger, more evolved (> 65 wt.% SiO 2 ), high Sr/Y Yarabamba and Quellaveco rocks. However, this amphibole fractionation signature would be imparted early in magma evolution, because amphibole generally crystallises early in experimental studies at lower crustal pressures (e.g., Nandedkar et al. 2014;Melekhova et al. 2017;Ulmer et al. 2018). A change in melt silica content from 60 to 65 wt.% SiO 2 would not change the fractionating assemblage from plagioclase-to amphibole-dominated, and, therefore, cannot solely explain the change in trace element chemistry. Instead, the stabilisation of amphibole in lower crustal cumulates at the expense of clinopyroxene would drive the derivative melts to granodioritic compositions (e.g., Klaver et al. 2018) which could explain the more evolved major element composition of the batholith after ~ 60 Ma.
We, therefore, interpret the inferred change in the fractionating assemblage as most likely related to an increase in Fig. 9 Core-rim Eu/Eu* vs. Hf composition of zircon crystals for (a) Toquepala Granodiorite (QVC083), (b) Quellaveco Granodiorite (QVC028), (c) Early Porphyry (QVC034). Red circles indicate core analyses, blue circles-rims, yellow-midway between core and rim. Lines join analyses from the same crystal and arrows point from core to rim pressure and melt water content, reflecting an increase in the depth of the locus of magma evolution in the lower crust. This could be the result of a stronger compressional tectonic regime, causing the timescale of magma storage and differentiation in the lower crust to increase and allowing residual melts to attain higher concentrations of incompatible trace elements (e.g., Annen et al. 2006) as well as water. The marked change in the geochemistry of magmas is broadly coincident with a clockwise rotation of the convergence of the subducting Nazca plate towards more orthogonal subduction at ~ 59 Ma which caused a significant increase in the rate of convergence (Pardo-Casas and Molnar 1987;Jaillard and Soler 1996). This initiated an episode of major compression in the Central Andes-referred to as the earliest phase of the Incaic orogeny, termed 'Incaic I' (59-55 Ma;Noble et al. 1985;Jaillard and Soler 1996;Benavides-Cáceres 1999). An equivalent tectonic phase in Northern Chile, the 'K-T' compressive event, has been attributed to the formation of a porphyry Cu belt consisting of Relincho (61 Ma), Spence (57 Ma) and Cerro Colorado (52 Ma) (Sillitoe and Perelló 2005;Mpodozis and Cornejo 2012).
The transition in melt compositions towards high Sr/Y and La/Yb predates the onset of the bulk of mineralisation in the district (~ 55 Ma; Sillitoe and Mortensen 2010) by ~ 3.5 ± 1.7 Myr. The delay suggests that despite the emplacement of magmas typically associated with porphyry Cu mineralisation in the shallow crust by ~ 59 Ma (the Quellaveco Granodiorite), this was insufficient to generate major porphyry systems at this time. Perhaps, the Quellaveco Granodiorite was sourced from an upper crustal magma reservoir that was not yet capable of supplying the heat and volume of fluid necessary for effective porphyry Cu mineralisation (Chelle-Michou et al. 2017). We envisage this extended time period as a necessary priming phase, where growth and recharge of the lower crustal magma reservoir allowed a sufficient volume of hydrous melt (> 800 km 3 ) to accumulate at depth (Rohrlach and Loucks 2005;Chiaradia and Caricchi 2017). Such magmas would be rapidly and episodically transferred to upper crustal levels (< 500 Kyr;Chelle Michou et al. 2017), where cooling and degassing would allow efficient migration and accumulation of a magmatic volatile phase at apices of the batholith (Huber et al. 2012;Parmigiani et al. 2016;Korges et al. 2020). Laterally focused fluid outbursts from the roof of the batholith would provide the source of melt and fluid to initiate and sustain the magmatic-hydrothermal system at Quellaveco (Lamy-Chappais et al. 2020).

The zircon record of the evolving Yarabamba Batholith and Quellaveco porphyry system
The trace element compositions of zircons not only record the trace element budget of the melts building the Yarabamba Batholith and Quellvaco porphyry Cu deposit, but also provide insights into the evolution of each intrusive phase from zircon saturation to the solidus. High precision zircon U-Pb studies have documented that zircons within a single sample of an intrusive or volcanic rock typically record several hundred kyr of magmatic history (e.g., Schoene et al. 2012;Wotzlaw et al. 2013;Samperton et al. 2015;Buret et al. 2016;Large et al. 2018). Some of these studies have resolved systematic trends in trace elements in zircon with time, indicating that zircon can record melt compositional evolution over such timescales. Furthermore, core to rim analyses on single grains provides a tool to track a relative liquid line of descent during zircon growth.
Here, zircon data between samples of different age (Yarabamba units vs. Quellaveco porphyry units, Fig. 8), and within individual grains (Fig. 9), show systematically increasing Yb/Gd, decreasing Th/U, increasing Hf and decreasing Ti. These trends are consistent with a previous zircon chemistry study of the Quellaveco intrusive complex (Simmons 2013) and are interpreted to reflect the changing trace element composition of the residual melt during zircon growth. Both increasing Yb/Gd and decreasing Th/U suggest the co-crystallisation of titanite and/or apatite which preferentially partition Gd over Yb and Th over U (Wotzlaw et al. 2013;Samperton et al. 2015;Buret et al. 2016;Loader et al. 2017;Rezeau et al. 2019). These ratios involve isovalent cations, which should minimise the effect of other potential influences on partition coefficients such as temperature and melt composition. Increasing Hf is consistent with progressive magma evolution, because Hf is incompatible in all common rock-forming minerals except zircon (Claiborne et al. 2006) and the partition coefficient of Hf would be expected to increase significantly with decreasing temperature (Blundy and Wood 1994). Titanium in zircon is typically considered as a proxy for the temperature of zircon crystallisation (Ferry and Watson 2007), thus the trend of decreasing Ti (which correlates with Th/U; Fig. 8b), could indicate progressive cooling of the melt after zircon saturation. Overall, these systematic trends are consistent with melt evolution towards more evolved compositions during cooling and crystallisation of the magmatic system after zircon saturation. Although systematic bivariate trends in zircon trace element chemistry are observed in some parameters, a more bimodal distribution is observed for Eu/Eu*, and less strongly for Ti (Fig. 8c, d), which increase and decrease, respectively after ~ 60 Ma. These two zircon populations correspond to samples that define the change in whole-rock chemistry after ~ 60 Ma (Fig. 7), with higher Eu/Eu* zircon in the younger rocks correlating with higher whole-rock Eu/ Eu* (Fig. 10). Based on our interpretation of the wholerock signatures, this coupling suggests that the change in Eu/ Eu* zircon after ~ 60 Ma is the result of a change in early, deep 12 Page 14 of 21 crustal processes, prior to magma emplacement, cooling and zircon saturation in the upper crust.
Lower Ti concentrations in zircons associated with porphyry Cu deposits compared with precursor granitoids has been observed elsewhere-at Yanacocha (Peru), Yerington (USA) and Highland Valley (Canada) (Dilles et al. 2015;Lee et al. 2020)-and may, therefore, be an indicator of petrogenetic processes that are linked to the formation of giant porphyry Cu systems. As noted above, the decrease in zircon Ti may reflect a decrease in the temperature of the melts from which zircon crystallised (Ferry and Watson 2007). It could be speculated that the inferred higher melt water contents of the magmas that formed the Quellaveco Granodiorite and the porphyry units, may have delayed zircon saturation to lower temperatures (Harrison and Watson 1983), however more recently it has been suggested that melt water content does not exert a major control on zircon saturation (Boehnke et al. 2013). Calculated temperatures using the Ti-in-zircon thermometer, assuming activities of TiO 2 and SiO 2 of 0.5 and 1.0, yield temperatures that range from 650 to 800 ºC for the rocks younger than ~ 60 Ma (Quellaveco rocks) and from 750 to 890 ºC for the Yarabamba rocks older than ~ 60 Ma (see Supplementary Material 2). The occurrence of many below/near-solidus temperatures suggests that other variables (e.g., variable activities of SiO 2 and TiO 2 , additional substitution mechanisms for Ti) may influence Ti concentrations in zircon (Fu et al. 2008). The melt TiO 2 activity can vary between different granitic melts and at different stages of melt differentiation (e.g., between 0.1 and 0.6; Schiller and Finger 2019) and can, therefore, influence calculated Ti-in-zircon temperatures. Lower Ti activity in the melt could be the result of stabilising amphibole and/ or biotite (in which Ti is compatible; e.g., Nandedkar et al. 2016;Azadbakht et al. 2020) as earlier and more abundant crystallising phases in the inferred higher water content melts after ~ 60 Ma. Therefore, these estimated temperatures should be treated with caution, and the origin of the marked decrease in zircon Ti concentrations at ~ 60 Ma could be related to temperature, a change in the TiO 2 activity of the melt or other unknown factors. Regardless of the origin, lower Ti concentrations in zircon appear to be a potential empirical discriminator of granitoids associated with giant porphyry Cu deposits.

Modelling the effects of melt compositions and oxidation state on zircon Eu/Eu*
Higher Eu/Eu* zircon in igneous rocks associated with porphyry Cu deposits, compared to precursor or unrelated igneous rocks, has been observed in several studies (e.g., Ballard et al. 2002;Liang et al. 2006;Dilles et al. 2015;Shen et al. 2015;Lee et al. 2017;Loader et al. 2017) and has been attributed by many authors to an increase in melt fO 2 . However, it is clear that fO 2 alone cannot account for all of the variability in Eu/Eu* zircon, because values above 1 are reported (e.g., Loader et al. 2017); a limiting value for highly oxidised, primitive melts containing essentially only Eu 3+ is 1. An alternative factor that can strongly influence Eu/ Eu* zircon that we consider in more detail below is the intrinsic Eu/Eu* melt , itself a function of prior magmatic evolution (e.g., Trail et al. 2012;Buret et al. 2016;Loader et al. 2017;Rezeau et al. 2019).
Theoretically, changes in the Eu anomaly of zircon can reflect changes in either: (1) the Eu anomaly of the melt; or (2) the zircon-melt partition coefficient for Eu (D Eu ) relative to the partition coefficients of its neighbouring REEs (D Sm and D Gd ). Experimental studies have demonstrated that increasing melt fO 2 increases D Eu towards D Sm and D Gd because of a concurrent increase in the proportion of Eu 3+ in the melt which partitions more strongly into zircon than Eu 2+ (Burnham and Berry 2012), thus an fO 2 control on Eu/Eu* zircon is expected. To assess the relative influence of changes in Eu/Eu* melt and fO 2 on Eu/Eu* zircon , a model was developed that calculates theoretical Eu/Eu* zircon as a function of Eu/Eu* melt , fO 2 , temperature and optical basicity (a measurement of melt structure) (Fig. 11). The partition coefficient of Eu into zircon (D Eu(zircon) ) is calculated in the model as a function of the partition coefficients for Eu 2+ and Eu 3+ into zircon (D Eu2+ and D Eu3+ ) and the Eu 2+ / Fig. 10 Eu/Eu* in zircon (density bars-horizontal lines show median values) and whole-rocks (small symbols), ranked by relative age as determined by geochronology for the Yarabamba samples and by cross-cutting relationships for the Quellaveco porphyry intrusions. Both zircons and whole-rocks show a proportional increase in Eu/Eu* from the emplacement of the Quellaveco Granodiorite and onwards Eu 3+ ratio of the melt. This is formulated as (Wilke and Behrens 1999;Aigner-Torres et al. 2007): where D Eu 2+ ~ 0, because Eu 2+ is essentially incompatible in zircon's structure due to its low charge and large ionic radius. The partition coefficients D Eu 3+ , D Sm and D Gd were calculated using lattice strain theory by compiling and parameterising the lattice strain fit parameters for zircon (Blundy and Wood 1994;Claiborne et al. 2018). The valence state of Eu in the melt (Eu 2+ /Eu 3+ melt ) was calculated as a function of fO 2 using an experimental calibration (Burnham et al. 2015).
To interrogate the sensitivity of Eu/Eu* zircon to Eu/ Eu* melt , a model was run in which Eu/Eu* melt was varied across an appropriate range of Eu/Eu* expected for arc magmas (0.1-1.2), and Eu/Eu* zircon in equilibrium with these melt compositions was calculated. The model used a Monte Carlo approach in which the value of each variable (Eu/ Eu* melt , temperature, optical basicity and fO 2 ) in each iteration of the model is stochastically generated within defined ranges. The model was run 10,000 times and the calculated Eu/Eu* zircon and the conditions in each model iteration were compiled (Fig. 11). Further information on the modelling is available in Supplementary Materials 1 and 3. This modelling identified a strong dependence of Eu/ Eu* zircon on Eu/Eu* melt (r 2 = 0.87; Fig. 11). Thus, using this relationship, it is possible to estimate Eu/Eu* melt for a natural dataset using measured Eu/Eu* zircon . By this method, Eu/ Eu* melt (at the time of zircon crystallisation) was estimated from the Eu/Eu* zircon measured for the early Yarabamba rocks (Yarabamba Monzonite) and the Quellaveco porphyries (see kernel-density plots in Fig. 11). The estimated Eu/Eu* melt for the Yarabamba Monzonite (interquartile range = 0.18-0.24) and the Quellaveco porphyries (interquartile range = 0.52-0.69) form two, nearly non-contiguous Eu/Eu* melt populations. These two melt populations display significantly lower Eu/Eu* than whole-rock Eu/Eu* for the Yarabamba Monzonite (interquartile range = 0.55-0.67) and the Quellaveco porphyries (interquartile range = 1.00-1.08), respectively, which is unsurprising given that melt compositions strongly diverge from bulk magma compositions during crustal transit and storage (e.g., Reubi and Blundy 2009). We interpret these low Eu/Eu* melt values, compared to whole-rock Eu/Eu*, as being due to modification of melt compositions prior to zircon saturation by plagioclase crystallisation during magma ascent and storage in the shallow crust (Sisson and Grove 1993;Waters et al. 2015). Despite this modification of the putative deep crustal signature, the model indicates that zircons in the Yarabamba Monzonite and the Quellaveco porphyries still inherited distinct Eu/ Eu* signatures, as also recorded by the whole-rock Eu/Eu*.
The strong control of Eu/Eu* zircon by Eu/Eu* melt predicted by the model does not preclude an additional influence of melt fO 2 . The model outlined above was re-run using a smaller Eu/Eu* melt range (0.12-0.70), based on the modelled Eu/Eu* melt derived for the Yarabamba Monzonite and Quellaveco porphyries. Melt fO 2 was found to have a weak effect on Eu/Eu* zircon with considerable scatter (Fig. 12), indicating that other parameters (such as Eu/Eu* melt ) have a stronger effect on Eu/Eu* zircon . Significantly, an increase in melt fO 2 from FMQ to FMQ + 3 (the maximum range for arc magmas; Carmichael 1991), cannot solely account for the difference in Eu/Eu* zircon between the early Yarabamba and the Quellaveco porphyry rocks. Furthermore, re-running the model to isolate the individual effect of each intensive parameter on Eu/Eu* zircon (Eu/Eu* melt , temperature, fO 2 and optical basicity) indicates that Eu/Eu* melt has a substantially Fig. 11 Heat map of modelled Eu/Eu* zircon (n = 10,000) as a function of a range of Eu/Eu* melt compositions using a Monte Carlo simulation, where Eu/Eu* melt , temperature (T), optical basicity (Ʌ) and fO 2 are varied stochastically in each iteration of the model within the given range. Darker colours illustrate higher density of modelled data points. Modelled Eu/Eu* zircon shows strong dependence on Eu/ Eu* melt (r 2 = 0.87). Kernel density plot along y-axis shows the density of natural Eu/Eu* zircon for Yarabamba Monzonite (dotted line) and the Quellaveco porphyries (dashed line) fit with a gaussian kernel function. The Eu/Eu* melt that is in equilibrium with these two zircon populations is estimated from the heat map and is displayed as a kernel density plot along the x-axis. The grey bars along the x-axis indicate the interquartile range of the whole-rock Eu/Eu* for Yarabamba Monzonite and the Quellaveco porphyries. The calculated Eu/Eu* melt is considerably lower than the whole-rock Eu/Eu* for both populations, indicating that zircons crystallised from a melt that was modified by plagioclase crystallisation without segregation greater levering effect on Eu/Eu* zircon than the other three parameters (see Supplementary Fig. S7).
The importance of increasing Eu/Eu* melt in driving the change in Eu/Eu* at ~ 60 Ma is consistent with the higher whole-rock Eu/Eu* after ~ 60 Ma (Fig. 10). Increasing melt fO 2 would only affect Eu/Eu* zircon and not the evolution of whole-rock Eu/Eu*, because fO 2 changes would not strongly alter the bulk Eu content of the magma. In addition, because whole-rock Eu/Eu* increases in tandem with redox insensitive trace element ratios (e.g., Sr/Y), it is thus best explained by fractionation of amphibole (excludes Eu relative to other REEs) and suppression of plagioclase fractionation (preferentially incorporates Eu relative to other REEs). Therefore, the shift to higher Eu/Eu* zircon at ~ 60 Ma is interpreted to reflect the composition of magmas that evolved at deeper crustal levels than previously, a signature which zircon inherits upon its crystallisation later in the magma evolution, in the shallow crust.

The decoupling of zircon chemistry from bulk magma compositions
The ability of Eu/Eu* zircon to record a deep-crustal fractionation history is dependent on it reflecting the bulk magma Eu/Eu* (cf., Nathwani et al. 2020). Although Eu/Eu* zircon appears to reflect changes in the inferred bulk magma Eu/ Eu* based on the correlation between Eu/Eu* zircon and whole-rock Eu/Eu*, modification of the melt composition (largely prior crystallisation of plagioclase in the uppercrust) prior to zircon saturation would be expected. The degree of this has implications for the applicability of Eu/ Eu* zircon as an indicator of long-lived, lower-crustal magma evolution.
Core-rim analyses of zircon provide a means to assess the trajectory and degree of late-stage melt evolution that may overprint the Eu/Eu* signature that zircons inherit from the prior magma evolution history. In the early Yarabamba rocks (older than ~ 60 Ma), Eu/Eu* zircon remains low from core to rim (Fig. 9a) which can be explained by Eu melt being initially low due to earlier, lower crustal plagioclase fractionation, and remaining low during zircon saturation and growth in the shallow crust (with or without plagioclase crystallisation). The variable core-rim vectors in the Quellaveco Granodiorite zircons (Fig. 9b) could reflect more open-system behaviour of the magma reservoir at this time, wherein periodic recharge of the system would reset Eu/ Eu* melt towards the higher Eu/Eu* of the magma batches being sourced from the deep crustal reservoir. In between recharge events, cooling of the reservoir and crystallisation of assemblages in which Eu is compatible (i.e., plagioclasedominated) would reduce Eu/Eu* melt and, therefore, also of Eu/Eu* zircon . In addition, in situ crystallisation of other accessory phases (titanite and apatite) which preferentially partition Sm and Gd relative to Eu, could have the opposite effect (e.g., Loader et al. 2017;Rezeau et al. 2019), leading to potentially complex trends. Finally, cycles of magma recharge, cooling and crystallisation could introduce greater melt compositional heterogeneity throughout the upper-crustal magma reservoir (Buret et al. 2016). In the Quellaveco porphyries, similar open-system behaviour is inferred; however, some samples demonstrate more consistently decreasing core-rim Eu/Eu* zircon values (Fig. 9c) which may represent cooling and crystallisation in the upper crustal magma reservoir between each phase of recharge that triggered fluid saturation and porphyry stock emplacement. The paucity of increasing core-rim Eu/Eu* zircon vectors in the Quellaveco porphyries (Fig. 9c) indicates that the process generating the high Eu/Eu* signature occurred prior to zircon saturation. This is consistent with magma evolution at depth producing high Eu/Eu* melt , which zircon inherits upon its late-stage crystallisation in the shallow crust.
Therefore, we conclude that the Eu/Eu* zircon record is a composite of both an early, lower crustal magma history, and a late, shallow crustal melt evolution that is dominated by plagioclase crystallisation. The Eu/Eu* zircon signature of deep crustal evolution will, therefore, be better preserved in earlier-crystallising zircons and/or in the cores of later crystallising zircons ( Fig. 8; lower Hf, higher Th/U) that more closely represent the bulk magma composition. Later crystallising zircons, would instead predominantly record the later in-situ melt evolution path during magma storage in Fig. 12 Heat map of modelled Eu/Eu* zircon (n = 10,000) as a function of a range of melt fO 2 (relative to the fayalite-magnetite-quartz redox buffer) using a Monte Carlo simulation, where Eu/Eu* melt , temperature (T), optical basicity (Ʌ) and fO 2 are varied stochastically in each iteration of the model within the given range. The range of Eu/Eu* melt selected is based on the modelled Eu/Eu* melt range for 5th and 95th percentile for Yarabamba and Quellaveco porphyries, respectively (Fig. 11). Darker colours illustrate higher density of modelled data points the upper crust, as seen by the lower Eu/Eu* in such zircons ( Fig. 8; higher Hf, lower Th/U). Although this effect does not appear to obscure the early, lower crustal magma history, it adds complexity to the development of Eu/Eu* zircon as a geobarometer or geohygrometer (e.g., Tang et al. 2020). The integration of Eu/Eu* zircon with indicators of melt differentiation in zircon (Hf and Th/U) in magma suites is, therefore, crucial in distinguishing and characterising early and late magmatic evolution processes.

Conclusions
Whole-rock and zircon chemistry and geochronology have been used to track the geochemical evolution of a long-lived, upper-crustal arc magma system in the ~ 12 Myr prior to, and during, the generation of giant porphyry Cu-Mo deposits in Southern Peru. Both whole-rock and zircon geochemical data document a distinctive shift in the geochemistry of the Yarabamba intrusive phases, ~ 4 Myr prior to the onset of significant mineralisation in the district. The rocks that represent magmas emplaced after this perturbation are characterised by high Sr/Y, Eu/Eu* and La/Yb, which we suggest were generated by increased water content and greater depth of magma evolution in the lower crust that promoted the stability of amphibole (± garnet) and suppressed plagioclase crystallisation. This geochemical change overlaps temporally with a change in the subduction vector of the Nazca plate from oblique to orthogonal and an increase in the rate of convergence, which are interpreted to have increased crustal compression and thickening and which could have impeded magma ascent from deeper levels.
Zircon records an increase in Eu/Eu* that occurred at the same time as the increase in Eu/Eu* observed in wholerock chemistry. We use mineral-melt partitioning theory to model Eu/Eu* zircon and show that the effect of melt fO 2 on Eu/Eu* zircon is subordinate to the effect of zircon crystallisation from a melt with elevated Eu/Eu*. Therefore, we propose that Eu/Eu* zircon is principally a monitor of Eu/Eu* melt that is predetermined by the lower crustal magma evolution recorded by whole-rock chemistry. Importantly, however, this fingerprint of lower crustal magma evolution in zircon can be partly overprinted by the prior or co-crystallisation of other minerals, predominantly by the onset of significant, upper crustal, plagioclase crystallisation.
Our study confirms previous studies that highlight the importance of high Sr/Y, high La/Yb magma compositions and lower crustal magma evolution in generating worldclass porphyry Cu deposits (e.g., Richards and Kerrich 2007;Loucks 2014;Chiaradia and Carrichi 2017). Our approach demonstrates the utility of zircon as a tracer of magmatic processes but emphasises that it is necessary to consider its crystallisation relative to other mineral phases in interpreting zircon trace element data for petrogenetic and provenance studies. support. We would particularly like to thank Victor Ramos for assistance during the field season in 2019. We thank Tobias Salge for help with scanning electron microscope work and Lauren Tuffield and William Brownscombe for laboratory support. This work benefited from thought provoking discussions with Matthew Loader. We are grateful to associate editor Othmar Müntener, Robert Loucks and an anonymous reviewer for their constructive reviews which greatly improved an earlier version of this manuscript. JW, SL and YB acknowledge funding under Natural Environment Research Council grant (NE/P017452/1) "From arc magmas to ores (FAMOS): A mineral systems approach".

Data availability
The reported datasets are provided in the electronic supplementary material.
Code availability Code for the Eu/Eu* modelling is provided as a Jupyter notebook in PDF format in Supplementary Material 3.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
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://creat iveco mmons .org/licen ses/by/4.0/.