Stochastic modelling of deep magmatic controls on porphyry copper deposit endowment

Porphyry deposits, our main source of copper and of significant amounts of Mo, Re and Au, form at convergent margins in association with intermediate-felsic magmas. Although it is accepted that copper is transported and precipitated by fluids released by these magmas, the magmatic processes leading to the formation of economic deposits remain elusive. Here we perform Monte Carlo petrological and geochemical modelling to quantitatively link crustal magmatic processes and the geochemical signatures of magmas (i.e., Sr/Y) to the formation of porphyry Cu deposits of different sizes. Our analysis shows that economic deposits (particularly the largest ones) may only form in association with magma accumulated in the lower-middle crust (P > ~0.5 GPa) during ≥2–3 Ma, and subsequently transferred to and degassed in the upper crust over periods of up to ~2.0 Ma. Magma accumulation and evolution at shallower depths (<~0.4 GPa) dramatically reduces the potential of magmatic systems to produce economic deposits. Our modelling also predicts the association of the largest porphyry deposits with a specific Sr/Y interval (~100 ± 50) of the associated magmatic rocks, which is virtually identical to the range measured in giant porphyry copper deposits.

Scientific RePoRTS | 7:44523 | DOI: 10.1038/srep44523 depths and timescales of formation of magmatic systems best suited to produce economic porphyry Cu deposits of different sizes within the time intervals constrained by radiometric dating (Table S1.1 in Supplementary Information 1).

Volumes of intermediate-felsic melts associated with porphyry Cu deposits. We focus on
Andean-type porphyry Cu deposits associated with syn-subduction arc magmatism. This magmatism has a rather restricted range of Cu concentrations 18 , thus, simple mass balance considerations set a lower limit to the minimum amount of felsic-intermediate melt required to provide the Cu precipitated in a porphyry deposit (see also refs 6 and 19). This implies that the larger the amount of Cu deposited, the larger must be the volume of magma that delivered the mineralising fluids. In order to calculate volumes of intermediate to felsic melts we have used the thermal models developed in refs 16 and 20. Specifically, the petrologic model of ref. 16 is used to quantify the volumes of intermediate to felsic magmas generated through time by concomitant fractional crystallization of parent hydrous basalts and partial melting of host rocks (amphibolites at lower crustal levels and graywackes at upper crustal ones) at different crustal levels. We considered a typical long-term arc magma flux of 0.0009 km 3 /a into magmatic systems located at depths ranging from 5 to 30 km (corresponding to pressures between 0.15 and 0.9 GPa: Fig. 1). The model results presented in ref. 20 are used to quantify volumes of intermediate to felsic melts generated at shallower depths (6-10 km corresponding to 0.18-0.3 GPa) by a higher magma flux of 0.016 km 3 /a (Fig. 1). These two end-member situations encompass the conditions and rates of magma transfer in the arc crust reasonably well.
Following the model of ref. 16 we simulate injection of mantle-derived basaltic melt as circular sills of 50 m thickness and 7.5 km radius every 10 ka at depths ranging between 5 and 30 km over timescales of up to several Ma. The injected basaltic melt cools and crystallizes at a rate that is inversely proportional to the depth of emplacement (assuming a geothermal gradient of 20 °C/km in the model). Progressive magma injection leads to a temperature increase of the surrounding rocks until a residual melt (from fractionation of the injected basalt) starts to accumulate. Eventually, the temperature of the magmatic system may reach values appropriate for partial  16). In the arc environment intermediate/felsic melts form as a consequence of injection of basaltic sills at different crustal depths. They result from residual melts issued from the fractionation of the injected basalt that eventually mix with partial melts of the host rocks, forming a hybrid melt. Because of thermal constraints hybrid melts generated at greater depths can grow in size much more than hybrid melts formed at shallower depths after a given time since the onset of basaltic sill injection. Hybrid melts formed at deeper crustal levels are H 2 O-undersaturated and must rise to shallower crustal levels to be able to exsolve fluid and form porphyry deposits. Intermediate-felsic magmatic systems generated by basaltic melts injected at shallow crustal levels under average magma fluxes may exsolve fluid readily but are generally limited in volume. Intermediate-felsic magmas injected at shallower levels under high magma fluxes result in the transient construction of large magmatic reservoirs that may readily exsolve fluids and could form porphyry deposits. However, the size of these potential deposits is limited by the duration of the high magmatic flux (< ~200 ka) 20,26 . melting of the surrounding rocks and a crustal melt starts to accumulate and mix with the residual melt (from basalt fractionation), forming a hybrid melt.
The most important results of the modelling are: 1) the mass of residual melt increases with the rate and duration of magma injection; 2) injected basaltic melts cool increasingly rapidly at shallower levels in the Earth's crust and partial melting of host rocks becomes increasingly difficult and less efficient, because the temperature of the host rocks steadily decreases with decreasing depth according to the geothermal gradient (20 °C/km). Therefore, melt productivity (the amount of hybrid melt accumulated after a certain time since the onset of the injection process with respect to the amount of total injected basaltic melt: Supplementary Information 2) decreases with decreasing depth; 3) over time, hybrid melts have compositions increasingly closer to that of the recharging basaltic magma, which is the opposite of what happens during the fractionation of a single magma batch 21 .
The model of ref. 20 is parameterised to evaluate melt productivity in shallow magmatic systems under short-lived (up to 200 ka), high magma fluxes (melt injection rate of 50 mm/a equivalent to 0.016 km 3 /a for a disk-shaped pluton with 10 km radius). The modelling approach is the same as the scenario for mantle-derived basaltic melts, with the exception that the injected magma is dacitic in composition and no partial melting of the host rocks occurs.
H 2 O content of the hybrid melts and Cu contents of their fluids. We calculated the H 2 O concentration of hybrid melts for both scenarios assuming that H 2 O in the parent basaltic melt varies between 2 and 4 wt.% and in the host crustal rocks between 0.2 and 1 wt.% (Table 1). While we considered a fluid phase consisting solely of H 2 O, we took into account the variation of H 2 O solubility in silicate melts with pressure and chemical composition 17 (Supplementary Information 2).
The amount of fluid needed to transport and precipitate Cu depends on the concentration of Cu in such a fluid. To calculate the concentration of Cu in the fluid phase, we first considered a range of Cu contents of magmas typical for thick volcanic arcs 17 (Supplementary Information 2), and then we performed Monte Carlo simulations for a range of fluid-melt partition coefficients for Cu between 2 and 100 (e.g., ref. 22). We assume that there is no recycling of pre-existing sulphide cumulates, or, if there is through partial melting and incorporation into the hybrid melt, this is computed in the Cu budget of arc magmas with normal Cu concentrations 18 . The Cu concentrations in the exsolvable fluid are on the order of a few hundreds of ppm (up to 700 ppm). This compares well with the Cu concentrations measured in intermediate density fluid inclusions from porphyry systems, also taking into account the high diffusivity of Cu in quartz, which might be responsible for a post-entrapment enrichment of measured Cu in intermediate density S-bearing fluid inclusions 23 .
Chemical composition of the hybrid melts. We have used an empirical relationship between melt productivity and SiO 2 to determine the SiO 2 content of the hybrid melt (see Supplementary Information 2 for further details). To calculate the Sr/Y ratios of the hybrid melts we have used available experimental petrology data 24,25 on hydrous basaltic parental melts undergoing fractional crystallizion at pressures ranging from < 0.1 to 1.2 GPa. The mineralogy as well as the proportions of minerals, and thus the bulk partition coefficients of Sr and Y, change both with pressure and residual melt fractions. For the fractionating minerals of the above experiments (olivine, spinel, plagioclase, orthopyroxene, clinopyroxene, amphibole, and garnet) we pooled all the available mineral/ melt partition coefficients for Sr and Y (GERM database: https://earthref.org/KDD/) subdividing them according to melt composition (equivalent to residual melt fractions in a process of fractional crystallization). The median values of the partition coefficients for both elements and for every mineral display systematic changes with the composition (or residual fraction) of the melt (Supplementary Information 2 and Dataset 1). We fit the partition coefficients as a function of residual melt composition and model their contribution to the evolution of Sr/Y in the residual melt as a function of pressure and amount of melt (Supplementary Information 2 and Dataset 1).

Monte Carlo simulations.
In the Monte Carlo simulations initial H 2 O and Cu content of the melt, fluid-melt partition coefficients for Cu, duration of sill injection at depth (t), and depth of magma accumulation (P) were contemporaneously and randomly varied within specified ranges for the two selected magma fluxes (Table 1; Supplementary Information 3, Dataset 2). We performed more than 100,000 simulations.

Results and Discussion
Melt injection at shallow crustal levels. High magma fluxes (0.016 km 3 /a) at shallow crustal depths Our results show that, under these conditions, magmas are H 2 O-saturated throughout the period of magma accumulation with more than 30 vol.% of excess fluid in > 90% of the simulations, which is sufficient for a continuous extraction of fluids 27 , potentially Cu-bearing (Fig. 2a). Assuming that these extracted fluids generate mineralization 28 , the median value of Cu accumulation after the maximum lifetime (~200 ka) permissible for such a sustained flux 20,26 would be < 9 Mt (< 4.5 Mt Cu for 50% efficiency: Fig. 2b). This could be a viable mechanism to form at best deposits with < ~5 Mt Cu (50% efficiency), but not larger ones (Fig. 2c). Interestingly, the Cu exsolution rates and Sr/Y values calculated for these short-lived magmatic systems are comparable to Cu emission rates and Sr/Y compositions of persistently degassing volcanoes (Etna 29 , Masaya 30 , White Island 31 , Stromboli 32 ; Fig. 2d). Mafic magma injection under average arc magma flux (0.0009 km 3 /a) at shallow crustal levels (P < ~0.4 GPa) prevents the accumulation of significant amounts of hybrid melt and dissolved Cu-bearing fluids (Cu < 1 Mt even after 3 Ma: Fig. 3a) because of the large amount of heat lost to the relatively cool host rock at shallow crustal levels 16 . The rate of potential Cu-release by these magmatic systems is significantly lower than the rates at which Cu is precipitated in porphyry systems (all Cu endowments of porphyry deposits are much higher for a given time than our model results in Fig. 3b). This indicates that shallow level accumulation and fractionation of mafic magmas under average arc fluxes is not a viable mechanism to generate giant porphyry deposits in the timeframes measured in natural deposits.
Basaltic melt injection at mid-to deep crustal levels. Basaltic melt injection at pressures > 0.4 GPa (> ~13 km) and with an average arc magma flux (0.0009 km 3 /a) results in a higher hybrid melt productivity because of the higher temperature of the surrounding rocks 16 (see km 3 of hybrid melt in Fig. 3a). Additionally, the H 2 O content of the hybrid melt is higher (Fig. 3a) because of the strong dependence of fluid solubility on pressure 17 . Therefore, basaltic melt injection at pressures > 0.4 GPa over times > 3 Ma, leads to the accumulation of hybrid melt containing sufficient H 2 O to deliver up to 10 Mt Cu (Fig. 3a). At accumulation pressures > 0.5 GPa and accumulation times > 3.5 Ma, the hybrid melt contains enough H 2 O to deliver > 30 and up to 240 Mt Cu (Fig. 3a). Our simulations indicate that magmatic systems capable of forming porphyry deposits with > 10 Mt Cu must develop at mid-to lower crustal depths (> ~17 km, P > 0.5 GPa) and require accumulation times > 2.5 Ma with an average arc magma flux (Fig. 3a,c,d). The modelling results are in agreement with the evidence that large porphyry systems develop at the end of magmatic cycles lasting several Ma 33-36 (2-5 Ma; Supplementary Information 1) with average magma fluxes 37 .
Transfer of magma from the deep accumulation reservoir to shallower crustal levels. Our results suggest that magmas containing the largest amounts of H 2 O and Cu (> 30 Mt) are systematically water-undersaturated because they formed at high pressures (> 0.5 GPa; Fig. 4a). Therefore, fluid-and Cu-bearing magmas accumulated in the lower to middle crust may rise towards shallower levels carrying the bulk of their water and Cu, before starting to release them upon reaching fluid saturation pressures. Calculations show that if these magmas rose adiabatically, the majority of the simulated fertile magmas would become H 2 O-saturated at pressures between 0.25 and 0.55 GPa (Fig. 4b). This pressure interval corresponds to depths of ~8-18 km, which are in the same range of the main magmatic reservoirs inferred for several porphyry-type deposits from geophysical, geochemical and mineralogical data 1,6,19,38 .
Thus, the formation of large accumulations of hybrid melts at mid-to lower crustal levels is an essential step in the formation of the porphyry system which sets the upper size limit of the potential deposit. However, it is the next step (the transfer of this large accumulation of hybrid melt to shallower levels from where fluids can be exsolved) that determines the formation and the ultimate size of the deposit (see also ref. 28). This is suggested by data on Cu endowment and overall duration of the ore deposition period for different deposits ( Supplementary Information 1, Table S1.1). Here, the overall duration of the ore deposition period was bracketed using only Re-Os dating of molybdenite from different hydrothermal pulses and/or U-Pb zircon dating of syn-and late/post-mineral porphyries (Table S1.1). Although the measurements of both Cu endowment and ore deposition period are subject to large uncertainties related to methodological approaches, uncertainties in cross-calibration of dating techniques and sampling bias, a broadly positive correlation exists between Cu endowment and the overall duration of the ore deposition period (Fig. 2c). Therefore, it is reasonable to consider that Cu in most porphyry deposits is precipitated at a long-term average rate of a few tens (~40) of tons of Cu per year. This does not mean that Cu is continuously precipitated at this rate. Rather, as pointed out by several  (Table S1.1). The onset of the mineralization has been placed in correspondence with the beginning of magma accumulation to allow comparison of the two timescales. For abbreviations of PCD see Table S1.1; (c) Cu (Mt) in the exsolvable fluid of hybrid melts after injection times of 0 to 5 Ma versus the overall duration of the ore deposition period. Color classes indicate different pressures of magma accumulation; (d) same as (c) but with colour codes indicating different times of magma accumulation. geochronological studies 9,10,39 , this suggests that larger deposits are formed by a higher number of short-lived (few tens of ka) pulses distributed over a longer timescale than in smaller deposits, with the long-term average rate of Cu deposition being similar in both cases. The overall longer lifetime of magmatic-hydrothermal activity in the larger deposits could be due either to a larger availability of hybrid melt in the deep reservoir of the larger deposits or to earlier tectonic interruption of the hybrid melt transfer from depth for the smaller deposits. The latter seem to be characterized by durations of precursor magmatic activity as long as those in the largest deposits ( Supplementary Information 1). Hence, it is possible that tectonics play an important role in modulating the transfer of magma from the deep to the shallow reservoir.

Sr/Y composition of magmatic rocks associated with porphyry Cu deposits. Our modelling
results quantitatively account for the association between supergiant porphyry Cu deposits (> 10 Mt of Cu at 50% efficiency) and magmas with Sr/Y values between ~50 and ~150 12 (Figs 5 and 6). In fact, in our model this range of natural Sr/Y values associated with major porphyry copper deposits corresponds to timescales (t) and crustal levels (P) of magma accumulation that are the most favourable for the development of the largest and H 2 O-richest magma volumes (Figs 3a and 5a,b). The availability of such a large amount of H 2 O-rich magma (with its corresponding "diagnostic" Sr/Y values) allows a long-lived transfer to shallower levels from where Cu can be released to form mineralization.
Our model suggests, in agreement with natural porphyry data, that for both higher (> 150) and lower (< 50) Sr/Y values, the probability to form deposits with > 10 Mt Cu is extremely low (Fig. 6), because total amounts of magmas and associated H 2 O formed at the corresponding P-t conditions are not sufficiently large (Figs 3a and  5a,b). Such limited H 2 O and magma quantities are not sufficient to sustain the longer-lived ore deposition period required to form supergiant deposits (Fig. 5c). In contrast, for magmatic systems with Sr/Y between 50 and 150 the probabilities to deliver 30-50 and 70-90 Mt Cu (50% efficiency) are ~9 and 4% of the simulations, respectively (Figs 5c and 6). The probability to form a deposit with > 90 Mt Cu (El Teniente-type) drops to about 0.2% (Figs 5c  and 6). Additionally, the largest number of simulations (66.5%) results in the production of porphyry deposits with < 10 Mt Cu, which is also observed in natural deposits (Fig. 6). The remarkable fit between model predictions and natural porphyry data (Figs 5c and 6) suggests that our model (based on measurable parameters like the overall duration of ore deposition period and average Sr/Y values of associated magmatic rocks; Fig. 5c) could be used for preliminary evaluation of mineral resources during early stages of mineral exploration.

Conclusions
Our results suggest that the formation of arc magmatic systems associated with porphyry Cu deposits occurs in two steps: 1) long-lasting (> ~2.5 Ma) injection of hydrous basalts in the mid-to lower crust (> ~17 km) leading to the formation of large amounts (> 800 km 3   lasting up to ~2 Ma. Whereas the first step (formation of large accumulation of hybrid melts at mid-to lower crustal levels) sets the upper size limit of the potential deposit, it is the next step (duration of the transfer of this large accumulation of hybrid melt to shallower levels from where fluids can be exsolved) that seems to determine the ultimate size of the deposit. The occurrence of these two steps also provides an explanation for the observed transition of tectonic regimes associated with the formation of giant porphyry systems in magmatic arcs, sometimes accompanied by a magmatic lull, from a mainly compressional stage (favouring multi-Ma accumulation of magma at depth) to a near neutral stress or slightly extensional stage 38,40 (promoting the migration of magma towards shallower depths and fluid exsolution within the timescales of porphyry deposit formation). The termination of the magma transfer during this second step could be due either to exhaustion of the mid-crustal magma reservoir or, more likely, to external tectonic processes. Magmatic systems that do not grow enough in size at depth cannot sustain the long-lived magmatic-hydrothermal activity that is recorded by the largest porphyry systems.
Overall, while not discounting the potential role of specialized magmas 18,41 , our model suggests that typical arc processes and melt compositions are able to explain the wide range of Cu endowments observed in porphyry deposits.  (Table S1.1). The colour codes correspond to different amounts of Cu potentially delivered by the magmatic systems; (c) Simulations of overall durations of the ore deposition periods versus Sr/Y values of associated magmas subdivided by classes of exsolvable Cu (Mt) with 50% efficiency (half of total exsolvable Cu). The model data (except Los Pelambres) show a good fit with natural porphyry system data (coloured circles: the colors of the circles indicate which Cu tonnage class reported in the legend each deposit belongs to; Rio Blanco is bicolor because, taking into account the whole district, Cu endowment may rise to > 200 Mt Cu).

Methods
Our model uses a Monte Carlo approach to simulate (> 100000 simulations) the amounts of hybrid melt produced in the crust (melt productivity = amount of hybrid melt accumulated/amount of total intruded basaltic melt), their water contents (both in solution and exsolvable), the Cu contents in the exsolvable water and the SiO 2 as well as Sr/Y compositions of the hybrid melts produced within the crust.
Melt productivity. We quantify melt productivity at different crustal depths under (i) an average arc magma flux [5 mm/a of basaltic melt injection rate through a circular section of 7500 m of radius, equivalent to 0.0009 km 3 /a] and (ii) an episodically high magma flux [50 mm/a injection rate of dacitic melt through a circular section with 10000 m radius, equivalent to 0.016 km 3 /a], using the models of refs 16 and 20.
For the average arc magma flux (0.0009 km 3 /a) we parameterized melt productivity for a time interval between 0 and 5 Ma and for pressures between 0.15 and 0.9 GPa (corresponding to crustal depths of ~5 to ~30 km) 16 .
We allowed pressure (P) and time (t) of maturation of the magmatic systems to vary randomly within the above mentioned fixed limits (Table 1) to obtain, for any random value of P and t, the corresponding value of melt productivity using a Monte Carlo method.
We have parameterized the curves of melt productivity from ref. 16   = + + + + + + x/y/z a' " "'OT b' " "'OT c' " "'OT d' " "'OT e' " "'OT f' " "'OT g' " "' where OT is the time since the onset of the injection of basaltic magma and a '/"/"' , b '/"/"' , c '/"/"' , d '/"/"' , e '/"/"' , f '/"/"' , g '/"/"' are constant values different for each one of the x, y and z variables. For the high magma flux (0.016 km 3 /a) scenario we parameterized melt productivity for a time interval up to 0.2 Ma (a maximum estimate for a continuously high magma flux 20 ) and for depths of magma emplacement shifting steadily through time from 6 to 10 km due to continuous injection of the sills underneath the previous ones 20 ( Supplementary Information 2, Figure S2.2). Volume of melt is calculated from volume of mobile magma assuming that the melt fraction in the mobile magma ranges randomly between 60 and 80% of the magma volume 20 . Simulation of the overall duration of the ore deposition period. The overall duration of the ore deposition period in our simulations was obtained considering the average Cu flux rate in natural porphyry Cu deposits which is provided by the broad linear correlation (Supplementary Information 2, Figure S2.7) between Cu endowment and overall duration of the ore deposition periods for selected porphyry Cu deposits (i.e., those with the best constrained radiometric dating available: see Table S1.1 in Supplementary Information 1). Although measurements of both these parameters may be subject to large uncertainties related to methodological approaches, uncertainties in cross-calibration of dating techniques and sampling bias, the positive correlation between them (Overall duration of ore deposition = 0.0202 * Cu (Mt), R 2 = 0.74, p < 0.00001) provides a reasonable first order approximation of the average rate at which Cu is precipitated in porphyry Cu deposits of variable sizes (~40 tons of Cu per year). This rate is used to quantify the time that is needed to flux the Cu contained in the exsolvable fluid associated with the variable amounts of hybrid melts produced by our simulations, using the equation = .
* Overalldurationoforedepositionperiod 0 0202 (50%exsolvable Cu) (6) assuming that 50% of the exsolved Cu is actually deposited. and at a temperature of 950 °C (typical for andesitic melts) was calculated using the online fugacity calculator at "https://www.esci.umn.edu/people/researchers/withe012/fugacity.htm" and parameterizing the variation of the molar volume of H 2 O with pressure by the following best fit equation where P = pressure (in GPa). Using temperatures of ± 100 °C has a minimal effect on the molar volume of H 2 O.