Resolving the timescales of magmatic and hydrothermal processes associated with porphyry deposit formation using zircon U–Pb petrochronology

Understanding the formation of economically important porphyry Cu–Au deposits requires knowledge of the magmatic-to-hydrothermal processes that act within the much larger magmatic system and the timescales on which they occur. We apply high-precision zircon geochronology (chemical abrasion–isotope dilution–thermal ionisation mass spectrometry; CA–ID–TIMS) and spatially resolved zircon geochemistry (laser ablation inductively coupled plasma mass spectrometry; LA-ICP-MS) to constrain the magmatic evolution of the underlying magma reservoir at the Pliocene Batu Hijau porphyry Cu–Au deposit. We then use this extensive dataset to assess the accuracy and precision of different U–Pb dating methods of the same zircon crystals.

Emplacement of the oldest pre- to syn-ore tonalite (3.736±0.023 Ma) and the youngest tonalite porphyry to cross-cut economic Cu–Au mineralisation (3.646±0.022 Ma) is determined by the youngest zircon grain from each sample, which constrains the duration of metal precipitation to fewer than 90±32 kyr. Overlapping spectra of single zircon crystallisation ages and their trace element distributions from the pre-, syn and post-ore tonalite porphyries reveal protracted zircon crystallisation together with apatite and plagioclase within the same magma reservoir over >300 kyr. The presented petrochronological data constrain a protracted early >200 kyr interval of melt differentiation and cooling within a large heterogeneous magma reservoir, followed by magma storage in a highly crystalline state and chemical and thermal stability over several tens of thousands of years during which fluid expulsion formed the ore deposit. Irregular trace element systematics suggest magma recharge or underplating during this final short time interval.

The comparison of high-precision CA–ID–TIMS results with in situ LA-ICP-MS and a sensitive high-resolution ion microprobe (SHRIMP) U–Pb geochronology data from the same zircon grains allows a comparison of the applicability of each technique as a tool to constrain dates and rates on different geological timescales. All techniques provide accurate dates but with different precision. Highly precise dates derived by the calculation of the weighted mean and standard error of the mean of the zircon dates obtained by in situ techniques can lead to ages of unclear geological significance that are older than the maximum ages of emplacement given by the CA–ID–TIMS ages of the youngest zircons in each sample. This lack of accuracy of the weighted means is due to the protracted nature of zircon crystallisation in upper crustal magma reservoirs, suggesting that standard errors should not be used as a means to describe the uncertainty in those circumstances. We conclude from this and similar published studies that the succession of magma and fluid pulses forming a single porphyry deposit and similarly rapid geological events are too fast to be reliably resolved by in situ U–Pb geochronology and that assessing the tempo of ore formation requires CA–ID–TIMS geochronology.

1 Introduction

Zircon geochronology is widely applied to date geological events and constrain timescales of geological processes. Combined with zircon geochemistry it has improved our understanding of crustal magmatic systems, such as those forming economically important magmatic–hydrothermal porphyry Cu–Au deposits. Advances in analytical techniques resulted in a shift from establishing the ages of magma emplacement or crystallisation to resolving the durations of magmatic and associated hydrothermal processes, such as magma accumulation or recharge, fractional crystallisation, or hydrothermal ore formation, and it has resulted in unprecedented information about the mechanisms and scales of magma ascent and storage in the earth's crust (e.g. Vazquez and Reid, 2004; Chamberlain et al., 2014; Barboni et al., 2016; Bucholz et al., 2017).

Porphyry copper deposits provide successively quenched samples of magma extracted from large crustal-scale hydrous magma systems. They are therefore a critical source of information about the processes and rates of magma ascent, magma storage and fluid generation, bridging those of volcanism and pluton formation. The identification of the processes that lead to porphyry deposit formation (e.g. Rohrlach et al., 2005; Audétat et al., 2008; Richards, 2013; Wilkinson, 2013) and the timescales on which they operate can provide us with valuable information about arc magmatic processes but could also potentially help in discriminating possibly fertile magmatic systems from ubiquitous infertile systems resulting in barren intrusions or volcanic eruptions.

Porphyry Cu–Au deposits commonly display clear field relationships of successive generations of porphyritic stocks or dikes, which were injected into sub-volcanic and other upper-crustal rock sequences (Sillitoe, 2010). The injected porphyry magmas thus provide snapshots of the underlying, vertically and laterally extensive magma reservoirs (e.g. Dilles, 1987; Steinberger et al., 2013). Cross-cutting relationships between veins and intrusive rocks suggest temporal overlap of hydrothermal alteration, ore mineralisation and porphyry emplacement (Proffett, 2003; Seedorff and Einaudi, 2004; Redmond and Einaudi, 2010). Strong hydrothermal alteration of the intrusive rocks associated with ore formation severely disturbs the geochemical information of most minerals and whole-rock compositions. While providing important insights into the hydrothermal history of a deposit (e.g. Roedder, 1971; Dilles and Einaudi, 1992; Landtwing et al., 2005; Cathles and Shannon, 2007; Seedorff et al., 2008; Large et al., 2016), it limits the investigation of the magma evolution, especially for the porphyries that are most intimately associated with ore formation. Zircon is a widespread mineral in intermediate to felsic rocks that is resistant to nearly all hydrothermal alteration and can thus provide coherent information about the evolution of a magmatic system. Recent advances in high-precision zircon geochronology by chemical abrasion–isotope dilution–thermal ionisation mass spectrometry (CA–ID–TIMS; e.g. Mattinson, 2005; Bowring et al., 2011; McLean et al., 2011a, 2015; Condon et al., 2015) now allow dating of the porphyritic intrusions associated with ore formation with unprecedented precision. The dramatically improved precision permits the constraint of rapid events, such as individual porphyry emplacement and hydrothermal mineralisation phases (<100 kyr; von Quadt et al., 2011; Buret et al., 2016; Tapster et al., 2016), that typically occur at the end of a longer-term period of volcanism and intrusive magma emplacement extending over several million years (e.g. Deino and Keith, 1997; Halter et al., 2004; Maksaev et al., 2004; Rohrlach et al., 2005; Lee et al., 2017). The integration of the temporal and chemical information gained from zircon is referred to as zircon petrochronology and can yield time-calibrated information about magma chemistry, thermal evolution and crystallinity during zircon crystallisation in magmatic systems (e.g Schoene et al., 2012; Chelle-Michou et al., 2014; Samperton et al., 2015; Buret et al., 2016; Szymanowski et al., 2017).

Timescales for magmatic and hydrothermal processes involved in porphyry ore formation have been suggested based on in situ U–Pb data (e.g. Garwin, 2000; Banik et al., 2017; Lee et al., 2017) and CA–ID–TIMS geochronology, which became increasingly precise (e.g. von Quadt et al., 2011; Chelle-Michou et al., 2014; Buret et al., 2016; Tapster et al., 2016; Gilmer et al., 2017; Large et al., 2018). However, several studies applying multiple techniques to the same sample sets have resulted in differing dates (von Quadt et al., 2011; Chiaradia et al., 2013, 2014; Chelle-Michou et al., 2014; Correa et al., 2016). The discrepancy demands a more detailed understanding of the precision and accuracy of the techniques and statistical data treatment that are applied to derive a geological age. This is not only fundamental for resolving dates and rates of geological processes in ore deposit research but provides more general insights into the geological meaning of magmatic dates and rates obtained by U–Pb geochronology.

For the present paper, we obtained a large dataset of zircon geochemistry and geochronology by laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) followed by high-precision geochronology of the same zircon crystals, utilising chemical abrasion–isotope dilution–thermal ionisation mass spectrometry (CA–ID–TIMS). These coupled data from the world-class Batu Hijau porphyry Cu–Au deposit allow us to resolve the chemical evolution and the changing physical state of the magma reservoir over time as well as the timescales of hydrothermal processes. In addition, previously published data on the same lithologies permit a critical comparison of two in situ microanalytical methods (sensitive high-resolution ion microprobe, SHRIMP, data by Garwin, 2000), LA-ICP-MS presented here) with high-precision U–Pb CA–ID–TIMS geochronology (this study). This allows us to critically compare the effects of variable degrees of precision and of the statistical treatment of data on the resulting interpreted ages, and it provides a means to test the accuracy of the different techniques.

2 Geological background

The Pliocene island-arc-hosted world-class porphyry deposit of Batu Hijau is located on the island of Sumbawa, Indonesia (Fig. 1), and it is one of the largest Cu and Au resources in the south-west Pacific region (7.23 Mt Cu and 572 t Au; Cooke et al., 2005). It is currently the only mined porphyry deposit in the Sunda–Banda volcanic arc, where Cu–Au porphyries are restricted to a narrow segment of the eastern Sunda–Banda arc from 115 and 120^∘ E (Fig. 1), where the Australian plate has been being subducted since the Eocene (Hall, 2002).

The exposed islands of the Sunda–Banda arc are characterised by late Oligocene to early Miocene calc-alkaline basaltic to andesitic arc rocks that are overlain or intruded by a late Miocene to Pleistocene calc-alkaline volcanic and plutonic rock suite ranging from basaltic to rhyolitic compositions (Hamilton, 1979; Hutchison, 1989). The magmatic arc hosts a variety of ore deposit types, including porphyry Cu–Au deposits; high-, intermediate- and low-sulfidation epithermal deposits; and a volcanic-hosted massive sulfide (VMS) deposit on Wetar (Fig. 1).

Figure 1 Tectonic map of South East Asia and the south-west Pacific. The Batu Hijau porphyry Cu–Au deposit (enlarged red square) is located on Sumbawa on the subduction-related magmatic Sunda–Banda arc within a small corridor between 110 and 120^∘ E that hosts several porphyry deposits. Arrows display plate motion relative to the Eurasian plate (Sundaland shield) but do not indicate velocity. HS: high sulfidation; IS: intermediate sulfidation; LS: intermediate sulfidation; VMS: volcanic hosted massive sulfide deposit. Most deposit locations are from Garwin et al. (2005).

The geology of the island of Sumbawa, hosting the Batu Hijau deposit, is dominated by early Miocene to Holocene volcanic arc successions deposited on oceanic crust that is 14–23 km thick (Hamilton, 1979; Barberi et al., 1987). Thickened continental crust observed in most other porphyry-mineralised magmatic arcs and commonly considered a prerequisite for porphyry Cu formation (Rohrlach et al., 2005; Chiaradia, 2015; Lee and Tang, 2020) is lacking beneath Sumbawa (Garwin et al., 2005). The distribution of volcano-sedimentary units, intrusions and the current coastline of Sumbawa are controlled by a major arc-transverse, left-lateral oblique-slip fault zone (Arif and Baker, 2004; Garwin et al., 2005). The fault zone strikes south-south-west to north-north-east about 30 km east of the Batu Hijau deposit, coinciding with the north-easterly projection of the Roo Rise oceanic plateau (Fig. 1).

Figure 2 Geological map (a) and north-to-south cross-sections with lithological information and grade contours (b–e) of the open pit at Batu Hijau. The intermediate and young tonalite intruded into a volcanic lithic breccia and the equigranular quartz diorite (a, b). Note that the extent of the old tonalite is not displayed but is included in the intermediate tonalite. The dashed line in panel (a) is the N–S section displayed in (b)–(e). Thin grey lines in panel (a) indicate mine benches. Cu and Au grades (c, d) are enveloped around the tonalites and a deep, central barren core. High-grade Cu and Au mineralisation is cut by the young tonalite. The ratios of Cu to Au (e) illustrate strong Au enrichment proximal to the intermediate tonalite and Cu-dominated distal mineralisation. Map, section and grades are based on company information from May 2016.

The hypabyssal stocks in the Batu Hijau district are intruded into an early to middle Miocene volcano-sedimentary rock sequence (<21 Ma based on biostratigraphy; Adams, 1984; Berggren et al., 1995) that reaches a thicknesses of up to 1500 m in south-western Sumbawa. The low-K₂O, calc-alkaline, sub-volcanic intrusive rocks in the Batu Hijau district have andesitic to quartz-dioritic and tonalitic compositions (Foden and Varne, 1980; Garwin, 2000) and were emplaced in several pulses during the late Miocene and Pliocene (Garwin, 2000). Over this multi-million-year magmatic history, a continuous geochemical evolution towards more fractionated lithologies is indicated by whole-rock chemistry and Fe-isotopic evidence of the magmatic rock suite in the Batu Hijau district (Garwin, 2000; Wawryk and Foden, 2017). Within the Batu Hijau deposit, andesite porphyries and different quartz diorite bodies are the earliest recognised stocks, whereas three tonalite porphyries are the youngest exposed intrusions (Clode, 1999). These tonalite porphyries, which are associated with economically important Cu–Au mineralisation and pervasive hydrothermal alteration at Batu Hijau, were emplaced as narrow, semi-cylindrical stocks into a broad east-north-east trending structural dome between ∼3.9 and 3.7 Ma (Fig. 2; Garwin, 2000). Based on petrography and cross-cutting relationships, they were termed old tonalite, intermediate tonalite and young tonalite (Fig. 3; Meldrum et al., 1994; Clode, 1999; Setynadhaka et al., 2008).

All three tonalite intrusions are petrographically similar and are geochemically described as low-K calc-alkaline tonalites (Idrus et al., 2007). The least altered specimens contain phenocrysts of plagioclase, hornblende, quartz, biotite and magnetite with or without ilmenite hosted in an aplitic groundmass of plagioclase and quartz (Fig. 3; Mitchell et al., 1998; Clode, 1999; Garwin, 2000; Idrus et al., 2007). Notably, all three porphyry intrusions lack potassium feldspar. Identified accessory minerals include apatite, zircon and rare titanites. Relicts of clinopyroxene can be identified within the tonalites. Vein density, ore grade and alteration intensity decrease from the old to young tonalite. The old tonalite is the volumetrically smallest, occurring mostly at the edges of the composite stock. It can clearly be identified in drill core, where its veins are truncated by later intrusions (Fig. 3), but it is currently not separated from the intermediate ionalite by the mine geology department at Batu Hijau because their phenocryst proportion is almost indistinguishable (Fig. 2). It locally contains the highest ore grades (>1 % Cu and >1 g t⁻¹ Au), and its matrix is characteristically the coarsest of the three tonalite intrusions. The intermediate tonalite is the volumetrically largest of the three porphyry intrusions and is strongly mineralised (Fig. 2). The intermediate tonalite is porphyritic with phenocrysts, including characteristic euhedral quartz phenocrysts, which are <8 mm in diameter (Fig. 3b). The young tonalite is the youngest intrusive rock, and it cuts most vein generations, ore mineralisation and alteration (Fig. 3c, e). It is strongly porphyritic with the largest-observed phenocrysts, including euhedral quartz phenocrysts, and contains elevated but subeconomic metal grades (<0.3 % Cu and <0.5 g t⁻¹ Au).

Figure 3 Rock specimens of the different tonalite porphyries and zircon CL images at Batu Hijau. Mineral assemblage in all tonalites is dominated by plagioclase, quartz and biotite. (a) Phenocrysts in the slightly propylitically altered, pre- to syn-Cu–Au-mineralisation, equigranular old tonalite are <3 mm. (b) Phenocrysts of the syn-Cu–Au-mineralisation, porphyritic intermediate tonalite are <5 mm. (c) The post-mineralisation porphyritic young tonalite contains the largest phenocrysts (<8 mm) and is characterised by a higher abundance of “quartz eyes”. (d) Abundant veins in the equigranular old tonalite are truncated by the later porphyritic intermediate tonalite. (e) Strongly veined intermediate tonalite is truncated by the barren and little-altered young tonalite. Dashed-green lines indicate intrusive contacts. (f) Representative zircons that display dominant oscillatory zoning and areas with little zoning. Circles indicate domains selected for LA-ICP-MS analyses (30 µm in diameter).

Copper and gold are systematically zoned within the deposit. High-Au zones are tightly enveloped around the tonalite stocks, whereas high copper grades extend further out into the volcanic lithic breccia and the equigranular quartz diorite (Fig. 2). The lowest Cu∕Au ratios occur towards and below the current pit floor, and higher Cu∕Au ratios are recorded peripheral to the central porphyry stock and towards the upper, already-mined part of the ore body. A positive correlation between vein density and Cu and Au contents was described at Batu Hijau (Mitchell et al., 1998; Clode, 1999; Arif and Baker, 2004). A veins were suggested to comprise ∼80 % of all quartz veins and contain a similar fraction of Cu (Mitchell et al., 1998). Most earlier authors suggested that the bulk of Cu and Au was precipitated as bornite during early A vein formation and converted to later chalcopyrite and gold associated with AB and B vein formation (Clode, 1999; Arif and Baker, 2004; Proffett, 2009). More recent studies on vein relationships and mineralogy using scanning electron microscopy cathodoluminescence (SEM-CL) petrography combined with fluid inclusion analyses suggest that Cu–Au ore mineralisation including bornite, chalcopyrite and gold all precipitated with a late quartz generation post-dating high-temperature A and AB vein quartz at a lower temperature together with the formation of brittle fractures associated with thin chlorite, white mica halos (∼ C veins and “paint veins”; Zwyer, 2011; Schirra et al., 2019). Irrespective of the debated timing of stockwork quartz veins and economic ore metal introduction in the old and intermediate tonalites, the young tonalite cuts through all high-grade Cu and Au zones, demonstrating its late, largely post-mineralisation emplacement (Figs. 2, 3e). Therefore, the maximum duration of economic mineralisation is bracketed by the emplacement ages of the Cu–Au-rich old tonalite predating it and the young tonalite post-dating it.

3 Materials and methods

Based on detailed core logging and outcrop mapping with company geologists in May 2016, one sample each from the old tonalite, intermediate tonalite, young tonalite and equigranular quartz diorite was selected from locations where the lithologies were in unequivocal relative time relationship (cross-cutting relationship; Fig. 3; see Supplement for sample locations). Rocks were crushed and zircons separated with conventional techniques, including Selfrag™ disintegration, panning and heavy liquid mineral separation (methylene iodide; 3.3 g cm⁻³). Selected zircons were annealed for 48 h at 900 ^∘C, mounted in epoxy resin and polished to reveal their crystal interior. Polished zircons were carbon-coated and imaged using SEM-CL (Tescan EOscan VEGA XLSeries 4 scanning electron microscope) prior to in situ LA-ICP-MS analysis for trace elements and U–Pb isotopes, employing a 193 nm ASI RESOlution (S155) ArF excimer laser with a 30 µm spot diameter, 5 Hz repetition rate and 2 J cm⁻² energy density coupled to an element SF-ICP-MS. A detailed description of the method including data reduction and results on secondary reference materials can be found in Guillong et al. (2014) and the Supplement. Generally, at least one spot was chosen in the interior (core) and one on the exterior (rim) part of the zircon, but up to four individual spots were analysed per zircon (Fig. 3f) to obtain in situ geochemical information and U–Pb dates. All ²⁰⁶Pb∕²³⁸U dates were corrected for initial ²³⁰Th–²³⁸U disequilibrium in the ²³⁸U–²⁰⁶Pb decay chain (e.g. Schärer, 1984). Ratios of Th∕U recorded by zircons cluster around 0.3–0.6, and the dates were therefore corrected assuming a constant Th∕U_melt of 2 based on partition coefficients (0.25) by Rubatto and Hermann (2007). Variation in the assumed Th∕U_melt by ±0.5 would result in changes in individual ²³⁸U–²⁰⁶Pb dates of <10 kyr, far below analytical uncertainty.

Titanium concentrations in zircon have been calibrated as a proxy for the crystallisation temperature of zircons (Watson and Harrison, 2005; Watson et al., 2006; Ferry and Watson, 2007) and have been widely used in igneous and ore deposit petrology (e.g. Claiborne et al., 2010b; Reid et al., 2011; Chelle-Michou et al., 2014; Dilles et al., 2015; Buret et al., 2016; Lee et al., 2017). The determination of accurate zircon crystallisation temperatures by Ti-in-zircon thermometry (Ferry and Watson, 2007) requires reliable estimates for the activity of SiO₂ and TiO₂ (aSiO₂ and aTiO₂) during zircon crystallisation. Based on previous studies on porphyry deposits we utilise an aSiO₂ of 1 and an aTiO₂ of 0.7 (Chelle-Michou et al., 2014; Buret et al., 2016; Tapster et al., 2016; Lee et al., 2017; Large et al., 2018), reflecting quartz and titanite saturation (Claiborne et al., 2006; Ferry and Watson, 2007). Titanite saturation during zircon crystallisation is ambiguous at Batu Hijau (see discussion), but changes in the assumed aTiO₂ result in systematic changes in all zircon crystallisation temperatures and will therefore not affect the interpretation of relative temperature changes: a change in the aTiO₂ of ±0.2 would result in a variation of about ±30 ^∘C.

Imaging by cathodoluminescence (CL), followed by low-precision but spatially resolved U–Pb dating and geochemical microanalysis by LA-ICP-MS, was used to evaluate potential inherited zircon populations and to select inclusion-free zircons for subsequent dissolution and analysis by high-precision U–Pb geochronology by CA–ID–TIMS. Selected crystals were removed from the epoxy mount and chemically abraded (CA) for 12–15 h at 180 ^∘C using techniques modified from Mattinson (2005). Zircons were spiked with 6–8 µg of the EARTHTIME ²⁰²Pb–²⁰⁵Pb–²³³U–²³⁵U tracer solution (ET2535; Condon et al., 2015; McLean et al., 2015) and dissolved in high-pressure Parr bombs at 210 ^∘C for >60 h. Dissolved samples were dried down and redissolved in 6N HCl at 180 ^∘C for 12 h. Sample dissolution, ion exchange chromatography modified from Krogh (1973) and loading onto zone-refined Re filaments were conducted at ETH Zürich and are described in detail by Large et al. (2018). High-precision U–Pb isotopic data were obtained by employing thermal ionisation mass spectrometry at ETH Zürich (Thermo Scientific TRITON Plus). Pb was measured sequentially on a dynamic MasCom secondary electron multiplier, and U was measured in static mode as U oxide using Faraday cups fitted with 10¹³ Ω resistor amplifiers (von Quadt et al., 2016; Wotzlaw et al., 2017). Data reduction and age calculation were performed using the algorithms and software described in McLean et al. (2011) and Bowring et al. (2011). All ²⁰⁶Pb/²³⁸U dates were corrected for initial ²³⁰Th–²³⁸U disequilibrium in the ²³⁸U–²⁰⁶Pb decay chain (e.g. Schärer, 1984) using a constant Th∕U partition coefficient ratio of 0.25 (Rubatto and Hermann, 2007), assuming that variations in Th∕U of the zircons result from different Th∕U of the crystallising melt and not from variations in relative zircon-melt partitioning of Th and U. High-precision U–Pb dates were obtained from 45 zircons, all of which were previously analysed by LA-ICP-MS.

4 Results

4.1 Optical zircon appearance and SEM-CL petrography

Zircon crystals from all three tonalite samples are colourless, euhedral to subhedral and variable in size, with c axis lengths of 100–500 µm and aspect ratios between 1:2 and 1:4 (Fig. 3f). Thin-section observations reveal zircons that are enclosed by phenocrysts and also occur within the fine-grained groundmass, suggesting protracted zircon crystallisation within the magma until emplacement of the tonalite porphyries. Investigation of mineral separates and mounts with a binocular microscope reveals that many zircons contain small (≪20 µm) mineral or melt inclusions. SEM-CL imaging reveals few non-zoned and sector-zoned zircon domains, but most zircons exhibit oscillatory zoning (Fig. 3f).

Only few broken zircon fragments could be identified from heavy mineral separates of the equigranular quartz diorite, but these indicate originally euhedral to subhedral shapes. Five of these broken grains, typically <200 µm long, could be identified and were mounted. Four zircons were non-zoned, and one was oscillatory-zoned.

4.2 Spatially resolved zircon trace element composition

At Batu Hijau, zircon geochemical analyses from the three tonalites display largely overlapping arrays and ranges for all analysed trace element concentrations and ratios (Fig. 4). Most zircons from the tonalites display systematically higher heavy rare earth elements (HREEs; e.g. Yb) over middle rare earth elements (MREEs; e.g. Dy) and light rare earth elements (LREEs; e.g. Nd) contents in their rims relative to their cores (Fig. 4a). This strongly correlates with core–rim systematics of other differentiation proxies, like increasing Hf or decreasing Th∕U (Fig. 4: Hoskin and Ireland, 2000; Claiborne et al., 2006, 2010b; Schaltegger et al., 2009; Samperton et al., 2015). However, some core–rim trends, especially from the young tonalite, display increasing Th∕U and decreasing Yb∕Dy ratios (Fig. 4b).

Figure 4 Covariation diagrams (a–e) and probability density plots (f–h) of in situ geochemical data obtained by LA-ICP-MS. Panels (a)–(c) are plotted against Th∕U as an indicator for fractionation, whereas panels (d) and (e) are plotted against Hf as the fractionation proxy. Arrows labelled “fractional crystallisation” indicate the approximate predicted direction zircon geochemistry would migrate given fractional crystallisation of zircon ± apatite ± titanite ± amphibole. Arrows labelled “plag” and “tit” point in the predicted direction of zircon geochemistry evolution during co-crystallisation with plagioclase or titanite. Zircons from the three tonalite porphyries are considered to have crystallised from the same magma reservoir, whereas zircons from the equigranular quartz diorite (purple) are unrelated (see text for discussion). Temperature lines in panel (d) are calculated with an assumed aSiO₂ =1 and aTiO₂ =0.7 based on Ferry and Watson (2007: see text for discussion). The crosses in the top right corners illustrate average analytical 2σ uncertainties. Probability density plots (after Vermeesch et al., 2013) illustrate differences between different samples and core and rim analyses within each sample. The axes of probability density plots in panels (f) and (g) are aligned with the axes of panels (d) and (e).

In most zircons, Ti concentrations decrease from core to rim (Fig. 4d, f). This decrease correlates well with increasing Hf and decreasing Th∕U. Maximum and minimum Ti contents for all intrusions are ∼10 and ∼2 ppm, resulting in model crystallisation temperatures of 770 to 650 ^∘C (see methods for details). The majority of zircons from the Batu Hijau deposit contain lower U concentrations (<75 ppm) compared to zircons from most other porphyry deposits (several 100 ppm), but individual zircons can contain up to 300 ppm (Fig. 4c). The zircons with high U concentrations do not correspond to the lower-Th∕U zircons but also contain high Th concentrations and cover the whole spectra of Th∕U ratios observed at Batu Hijau (Fig. 4c). The Eu anomaly (Eu ∕ Eu^*, which is a means to quantify the negative inflexure of the normalised rare earth element, REE, diagram), increases (Eu ∕ Eu^* decrease) with increasing Hf concentration (Fig. 4e). Zircon analyses from the equigranular quartz diorite plot towards the lowest Hf, Yb ∕ Dy, Yb ∕ Nd and Eu ∕ Eu^* but highest Th∕U and Ti end of the trends displayed by all tonalite zircons (Fig. 4).

Probability density functions (Vermeesch, 2012) are used to test for statistically significant differences between the overlapping zircon populations of the different tonalites and between core and rim analyses from the same tonalite porphyries (Fig. 4f, g, h). The Hf and Ti concentrations as well as the Eu anomaly of zircons display overlapping distributions for the intermediate and young tonalites. The old tonalite zircon population peaks at higher Ti concentrations and Eu ∕ Eu^* as well as lower Hf concentrations than the younger tonalites. Core and rim analyses from zircons of the old tonalite document decreasing Ti and Eu ∕ Eu^* together with increasing Hf concentrations from cores to rims. Hafnium contents of the rim analyses peak at higher concentrations than the core analyses within the intermediate and young tonalite, whereas Eu ∕ Eu^* displays the opposite effect. Populations illustrating titanium concentrations of the two younger tonalites, however, display no systematic changes between core and rim.

4.3 CA–ID–TIMS geochronology

We dated 16 zircons of each the old and intermediate tonalite and 13 zircons of the young tonalite by high-precision CA–ID–TIMS geochronology. The youngest zircon each of the old, intermediate and young tonalite yields ²³⁰Th⁻²³⁸U disequilibrium-corrected ²⁰⁶Pb∕²³⁸U zircon dates of 3.736±0.023, 3.697±0.018 and 3.646±0.022 Ma (individual grain $\pm \mathrm{2}\mathit{\sigma },$ Fig. 5). We interpret these dates as the time of respective porphyry emplacement based on the assumption that zircons grew at depth up to the point that the magma cooled rapidly upon injection into the sub-volcanic environment (cf. Oberli et al., 2004; Schaltegger et al., 2009; von Quadt et al., 2011; Samperton et al., 2015; Large et al., 2018). Reproducibility of the individual dates is corroborated by reproducible dates for all analysed CA–ID–TIMS standards analysed during the time of analyses of the presented samples. Standards include the Aus_Z7_5 (von Quadt et al., 2016), which contains similarly low amounts of radiogenic Pb (0.5–4 pg) and is of Pliocene age.

Despite the small differences in emplacement age, the age sequence is consistent with field observations documenting the emplacement sequence (Fig. 3). The time intervals between emplacement of the old and intermediate tonalite and between the intermediate and young tonalite can therefore be constrained to 39±29 and 51±28 kyr, respectively. The minimum duration of zircon crystallisation, as recorded by the oldest and youngest zircon of each sample, spreads over 246±28, 212±32 and 171±26 kyr for the old, intermediate and young tonalite (Fig. 5). The overall duration of recorded zircon crystallisation is 336±27 kyr. Using the youngest zircon population rather than the youngest individual zircon as the best approximation for porphyry emplacement (cf. Samperton et al., 2015; Buret et al., 2016; Tapster et al., 2016) would result in slightly older emplacement ages (∼20 kyr) but nearly identical durations of zircon crystallisation and time intervals between porphyry emplacement events (see Supplement). Our high-precision CA–ID–TIMS dates precisely constrain protracted zircon crystallisation over several hundred thousand years and successive emplacement of the three porphyritic tonalite bodies at Batu Hijau within 90±32 kyr.

Ratios of Th∕U obtained by CA–ID–TIMS analyses on the same sample volume illustrate no systematic variation with time. Values vary inconsistently between 0.4 and 0.6 over the whole recorded time interval (Fig. 6).

Figure 5 High-precision U–Pb CA–ID–TIMS zircon dates from the three tonalite porphyries. Vertical bars are individual analyses including analytical uncertainty (2σ). The youngest crystallisation age is used as the best approximation for porphyry emplacement. The extended range in zircon crystallisation ages in each sample indicates protracted crystallisation. The yellow box indicates maximum duration of ore formation as constrained by the emplacement age of the old tonalite and the young tonalite. The grey box illustrates the duration of zircon crystallisation recorded by CA–ID–TIMS geochronology. The vertical bars in the grey box are emplacement ages of the tonalites, demonstrating >200 kyr of zircon crystallisation before emplacement of the first porphyry intrusion and start of Cu–Au mineralisation.

4.4 In situ U–Pb geochronology

Trace element and U–Pb isotopic data were obtained for each LA-ICP-MS spot (Fig. 7) prior to CA-IC-TIMS dating. Low uranium concentrations and the young ages of the analysed zircons resulted in high individual uncertainties for individual in situ U–Pb dates (mean: 10 %; minimum: 3 %; maximum: 41 %). All individual spot analyses of the three tonalites that were not discarded due to common Pb or strong discordance yield Pliocene dates (2.98±1.06–4.95±0.54 Ma; Fig. 7) with no apparent inherited zircons. All in situ dates of individual samples illustrate continuous arrays and do not indicate more than one population of zircons per sample (Fig. 7). Weighted means of all zircon analyses from each tonalite are 3.879±0.027, 0.065, 0.32 (n=207; mean square weighted deviation, MSWD =2.1), 3.778±0.023, 0.061, 0.62 (n=189, MSWD =2.5) and 3.751±0.023, 0.060, 0.29 Ma (n=158, MSWD = 2.6) from oldest to youngest (Fig. 7), where the stated uncertainties are the standard error of the weighted mean, the standard error including an external uncertainty of 1.5 % as suggested by Horstwood et al. (2016) to incorporate excess variance, and the standard deviation of zircon dates from each sample. These weighted averages do not overlap within uncertainty with the emplacement ages constrained by CA–ID–TIMS but overlap with the mean of the respective population. The few in situ analyses (n=8) on zircons from the diorite result in overlapping late Miocene dates. The weighted mean of all LA-ICP-MS analyses of the equigranular diorite results in an apparent age of the equigranular diorite of 6.37±0.40, 0.41, 0.37 Ma (n=8, MSWD =0.46).

Figure 6 Th∕U ratios plotted against time. Both values obtained from CA–ID–TIMS analyses of the same sample volume.

Garwin (2000) previously published SHRIMP U–Pb data on zircons from the Batu Hijau tonalites. Similar to the LA-ICP-MS analyses in this study, individual uncertainties of the dates were elevated (0.12–0.30 Ma, ∼5 %–10 % uncertainty) due to low U concentrations and the young zircon crystallisation ages. As all dates of each sample statistically formed single populations (Supplement), weighted means were calculated of all zircons of a sample; these were interpreted as the intrusion ages of the tonalites by Garwin (2000). The reported zircon dates were not corrected for ²³⁰Th–²³⁸U disequilibrium. For comparability we will only consider zircon dates that are corrected for initial Th∕U disequilibrium (Schärer, 1984; for details consult the Supplement). Correction increases the individual zircon dates by ∼60–100 kyr, and recalculation of the weighted mean averages and standard errors results in dates of 3.74±0.14 Ma (MSWD =1.2, n=8), 3.843±0.094 (MSWD =1.2, n=18) Ma and 3.81±0.2 Ma (MSWD =2.35, n=7) for the old, intermediate and young tonalite, respectively.

Figure 7 In situ U–Pb geochronology by LA-ICP-MS of zircons from the equigranular quartz diorite (a), the old tonalite (b), the intermediate tonalite (c) and the young tonalite (d). Vertical lines illustrate individual U–Pb dates including analytical uncertainty (2σ). As no zircon populations can be separated, the weighted mean average of all analyses is calculated. Standard error (2SE, lightest grey), standard error +1.5 % (light grey) to incorporate excess variance (Horstwood et al., 2016) and standard deviation of all individual dates (dark grey) are calculated and plotted. Duration of zircon crystallisation as obtained by CA–ID–TIMS is illustrated as a box for comparison. Lower boundary of the coloured box indicates youngest CA–ID–TIMS date or the emplacement age. Note the different vertical scale in panel (a). Zircon dates from the equigranular quartz diorite are ∼2 Myr older than zircons from the tonalites.

5 Discussion

5.1 Timing and duration of magmatic and hydrothermal processes leading to porphyry Cu formation

The three tonalite intrusions each record protracted zircon crystallisation over ∼200 kyr, as resolved by high-precision ID-TIMS geochronology. The older zircon dates from the young and intermediate tonalites overlap with the younger zircons of the older intrusion(s) (Fig. 5). This overlap, together with the largely overlapping trace element systematics recorded by zircons, is used here to infer zircon crystallisation within the same mid- to upper crustal magma reservoir that sourced magmas forming the three tonalitic porphyry stocks but most likely also volatiles and metals to form the porphyry Cu–Au deposit. High-precision geochronology records a total duration of zircon crystallisation of 336±27 kyr, which is also a minimum estimate for the lifetime of the deeper reservoir underlying Batu Hijau. The first exposed and highly mineralised tonalite intrusion (old tonalite) was injected into the upper crust 246±28 kyr after the onset of zircon crystallisation. Emplacement of the three tonalites occurred within 90±32 kyr. Emplacement of the old tonalite was followed by the emplacement of the intermediate tonalite after 39±29 kyr, and the young tonalite was emplaced after a further 51±28 kyr.

The maximum duration of ore formation is defined by the time span between the emplacement of the pre- to syn-mineralisation old tonalite and the post-mineralisation young tonalite (Fig. 3d, e) and can be therefore constrained to less than 122 kyr. This maximum duration is in good agreement with previous geochronological studies indicating timescales of ore formation from <100 to <29 kyr (Figs. 8, 9; von Quadt et al., 2011; Buret et al., 2016; Tapster et al., 2016). It is also coherent with results from thermal modelling studies (Cathles, 1977; Weis et al., 2012) and modelling of diffusive fluid–rock equilibration (Cathles and Shannon, 2007; Mercer et al., 2015; Cernuschi et al., 2018) suggesting timescales of ore formation between a few thousand years and 100 kyr. Strongly elevated Cu and Au grades in the old tonalite and somewhat lower, but still economic, grades within the intermediate tonalite (Clode, 1999; Garwin, 2000; Arif and Baker, 2004) together with cross-cutting relationships (Fig. 3d, e) indicate that mineralisation occurred within at least two but possibly more pulses: (i) one strong mineralisation pulse associated with or slightly post-dating the emplacement of the old tonalite but predating the injection of the intermediate tonalite (Fig. 3d) and (ii) a second pulse is bracketed by the intrusion of the intermediate and the young tonalite (Fig. 3e). More than one episode of mineralisation is also inferred based on detailed mineralogy and vein petrography (see Geological background section; Arif and Baker, 2004; Zwyer, 2011). This further strengthens the hypothesis that individual ore-forming hydrothermal pulses are relatively short events, possibly on the millennial or sub-millennial scale (Cathles, 1977; Weis et al., 2012; Mercer et al., 2015), but that the formation of large economic Cu–Au deposits occurs in several pulses over tens of thousands of years but ≤100 kyr (von Quadt et al., 2011; Weis et al., 2012; Cernuschi et al., 2018).

5.2 Reconstructing the chemical and physical evolution of a porphyry-forming magma reservoir

Th∕U ratios and Hf concentrations are commonly used as proxies for the degree of crystal fractionation within a magma reservoir (e.g. Claiborne et al., 2006, 2010b; Samperton et al., 2015). The decreasing Th∕U ratios and increasing Hf concentrations between samples and from cores to rims (Fig. 4) are indicative of progressive melt differentiation during zircon crystallisation. The good correlation of these melt evolution proxies with decreasing Ti contents (Fig. 4) further suggests progressive cooling during differentiation. Ratios of HREE over MREE or LREE (e.g. Yb/Dy, Yb/Nd) can be utilised to make inferences about the co-crystallising mineral assemblage. Titanite for example preferentially depletes the melt in MREE, resulting in distinct trace element patterns recorded by co-crystallising zircon (e.g. Reid et al., 2011; Wotzlaw et al., 2013; Samperton et al., 2015; Loader et al., 2017). The systematically higher HREE (e.g. Yb) over MREE (e.g. Dy) and LREE (e.g. Nd) contents in the rims of most zircons relative to their cores (Fig. 4a) thus indicate zircon crystallisation from a fractionally crystallising magma with co-crystallisation of minerals that preferentially incorporate MREE and LREE (e.g. apatite, titanite). At Batu Hijau apatites were petrographically identified, whereas magmatic titanite occurs very subordinately. The apparent lack of magmatic titanite is unusual as it is reported as a common accessory phase in many other porphyry Cu deposits (e.g. Bajo de la Alumbrera, El Salvador, Ok Tedi, Oyu Tolgoi). The absence of euhedral titanite within the mineral separates could be the result of dissolution during intense hydrothermal alteration (van Dongen et al., 2010).

Trace element compositions of zircons from the equigranular quartz diorite suggest crystallisation within a hotter and less evolved magma than the zircons from the tonalites (Fig. 4). In principle, this might indicate that all zircons analysed in this study have crystallised from the same reservoir, where the zircons from the equigranular quartz diorite reflect the earliest crystallised zircons from least evolved melt. However, the >2 Myr time gap is longer than the thermal lifetime of any recognised upper-crustal magmatic body (e.g. Schoene et al., 2012; Wotzlaw et al., 2013; Caricchi et al., 2014; Samperton et al., 2015; Eddy et al., 2016; Karakas et al., 2017) and longer than considered possible based on thermal modelling (Jaeger, 1957; Annen, 2009; Barboni et al., 2015). We therefore consider the zircons within the equigranular diorite to be part of a separate crustal magmatic system not directly related to the ore-forming system that sourced the three tonalitic intrusions.

Trace element populations of zircons from the three tonalites demonstrate that the crystallising magma at the time of emplacement of the old tonalite was hotter and less fractionated (Fig. 4) than at the time of emplacement of the younger intermediate and young tonalite (i.e. 39±29 and 90±32 kyr after emplacement of the old tonalite, respectively). The good correlation of proxies indicating progressive differentiation (Th∕U and Hf) with decreasing Ti concentrations (Fig. 4d) indicates that the magma reservoir cooled during concurrent crystallisation and melt evolution. In situ analyses of cores and rims are evidence of an evolving magma reservoir over the course of individual zircon crystallisation (decreasing Hf; Fig. 4h). Core–rim systematics of zircons from the old tonalite further demonstrate cooling during protracted zircon growth (Fig. 4f). Rarely recorded coherent zircon trace element systematics recording melt differentiation over time are commonly inferred to result from zircon crystallisation within a homogeneous magma that evolved continuously (e.g. Wotzlaw et al., 2013; Large et al., 2018). The lack of such systematic temporal changes in the chemistry of the zircons (Fig. 6) indicates that the magma reservoir at Batu Hijau was not evolving homogenously. This could be explained by incremental recharge or assembly of the magma reservoir. However, this would imply at least partial resetting of the intra-grain systematics recorded in zircons from the old tonalite (Buret et al., 2016; Large et al., 2018). To explain the intra-grain and inter-sample systematics but absence of temporal trends (Figs. 4, 6), we favour different degrees of crystallinity in the magma reservoir. Overall, the reservoir is generally hotter and less evolved at the time of emplacement of the old tonalite than thereafter (Fig. 4). We therefore suggest that the magma reservoir underlying Batu Hijau progressively but heterogeneously cooled and crystallised over at least 246±28 kyr, with potential incremental recharges until emplacement of the old tonalite.

A change from a differentiating, crystallising and cooling magma reservoir to a state of chemical and thermal stability is recorded between emplacement of the old and young tonalite (separated by 90±32 kyr) as demonstrated by the trace element systematics of the intermediate and young tonalite porphyries. The indistinguishable, highly fractionated and low-temperature zircon characteristics (Fig. 4) indicate that the magma reservoir remained in near-steady-state conditions between emplacement of the old and young tonalite as coherent intra-grain systematics are not pronounced (Hf) or are absent (Ti) in zircons from the younger tonalites (Fig. 4f, h).

Irregular zircon trace element systematics in other intrusive magmatic settings have been associated with crystallisation in non-homogenised and small melt batches sometimes, with contemporaneous incremental magma addition to the mushy magma reservoir (e.g. Schoene et al., 2012; Buret et al., 2016; Tapster et al., 2016). Geochemically similar zircon chemistries of the intermediate and young tonalite could also result from chemical stability as the magma reservoir reached the “petrological trap” at a crystallinity of ∼55 %–65 % (Caricchi and Blundy, 2015), where the crystal fraction does not change over a broad temperature interval. Rim analyses that plot more outside the mineral co-crystallisation trends than the respective core analyses (Fig. 4) could suggest late-stage crystallisation within a nearly solidified magma that can be characterised by unsystematically variable trace element systematics (Buret et al., 2016; Lee et al., 2017). Alternatively, they could indicate thermal and possibly chemical rejuvenation of the magma (Buret et al., 2016). The latter would help explain the recorded thermal stability over tens of thousands of years. It is not possible to unambiguously identify one of the two mechanisms as dominant, and a concurrence of both is feasible. We therefore propose that in between emplacement of the old and young tonalite the underlying magma reservoir was in a thermally and chemically stable and crystal-rich state and was most likely affected by incremental magma recharge or underplating.

Our data of a porphyry-Cu-fertile magmatic system constrain a heterogeneous magma reservoir that was initially dominated by cooling and melt differentiation and evolved into a thermally and chemically stable, crystal-rich magma that possibly experienced incremental recharge. The likely transitional change of reservoir behaviour can be temporally constrained to have occurred between emplacement of the old and young tonalites and coincides with the formation of a world-class Cu–Au reserve. This suggests that porphyry Cu–Au deposits form after a few hundred thousand years of cooling and crystallisation, potentially within an originally melt-rich magma reservoir.

5.3 Different timescales of processes related to porphyry Cu–Au ore formation

To date no clear relationship between the duration of magmatic–hydrothermal activity and the size of porphyry deposits can be identified from studies applying high-precision CA–ID–TIMS geochronology. Comparison of published datasets (Buret et al., 2016; Tapster et al., 2016; Large et al., 2018) reveals maximum durations of metal forming events of between tens of thousands of years and 10⁵ yr (Fig. 9). Although these studies are so far constrained to deposits of <10 Mt of contained Cu, they range over at least 1 order of magnitude in size (Koloula vs. Batu Hijau). A correlation between the duration of the mineralising or magmatic event and the total mass of deposited copper had been previously suggested based on compilations of different geochronological datasets (Chelle-Michou et al., 2017; Chiaradia and Caricchi, 2017; Chelle-Michou and Schaltegger, 2018; Chiaradia, 2020). High durations of ore formation (>1 Myr) were suggested based on Re–Os geochronology on molybdenite at the giant porphyry deposits and deposit clusters in Chile (>50 Gt Cu; El Teniente, Cannell et al., 2005 and Maksaev et al., 2004; Rio Blanco, Deckart et al., 2012; and Chuquicamata, Barra et al., 2013). Mineralising timescales of copper (and gold) were calculated by subtracting the youngest from the oldest Re–Os date. However, recent Re–Os dates from El Teniente (Spencer et al., 2015) indicate that the spread in dates is more consistent with several short (≤200 kyr) hydrothermal events separated by hiatuses of ∼500 kyr. Thus, the large tonnage of these deposits could be the result of the superimposition of several ore-forming mid- to upper crustal magmatic systems. As the correlation of deposit size and timescales of shallow magmatic–hydrothermal systems is currently ambiguous, we would argue that other variables could be the dominant factors controlling the deposit size, such as magma reservoir size, magma or fluid chemistry, fluid release, and focusing mechanisms or the metal precipitation efficiency.

Zircon crystallisation over ∼200 kyr before the onset of porphyry ore formation recorded at Batu Hijau is consistent with other high-precision geochronological studies on porphyry deposits (Figs. 8, 9: Buret et al., 2016; Tapster et al., 2016; Large et al., 2018). The lack of variation observed in these deposits suggests the necessity of a long-lived and continuously crystallising magma reservoir preceding economic ore formation. The recorded ∼200 kyr of protracted zircon crystallisation could indicate a period of volatile enrichment as a result of fractional crystallisation and cooling of the magma reservoir before porphyry emplacement.

Figure 8 Compilation of high-precision datasets on several pre-, syn- and post-ore intrusions at magmatic–hydrothermal Cu–Au deposits. Data for Ok Tedi, Bajo de la Alumbrera and Koloula are from Large et al. (2018) Buret et al. (2016) and Tapster et al. (2016), respectively. Coloured vertical bars are individual analyses including analytical uncertainty (2σ). Intrusions are categorised as pre-ore, pre-/syn-ore and post-ore intrusion. Decreasing deposit size from left to right (tonnages from Cooke et al., 2005).

The geochronological data from the Batu Hijau district are further evidence that rapid porphyry emplacement and ore formation (<100 kyr) are the product of a longer-term evolution (a few 100 kyr) of a large magma reservoir underlying the porphyry deposit that is the main driver of ore formation (von Quadt et al., 2011; Chelle-Michou et al., 2014; Buret et al., 2016, 2017; Tapster et al., 2016; Large et al., 2018). Magma reservoirs capable of forming porphyry deposits are in turn part of a longer-term (several Myr) evolution of lithosphere-scale magma systems (Sasso, 1998; Rohrlach et al., 2005; Longo et al., 2010; Rezeau et al., 2016), which is consistent with the ≫2 Myr record of intrusive rocks preceding porphyry emplacement and ore formation recorded in the Batu Hijau district (Garwin, 2000; Wawryk and Foden, 2017).

Figure 9 Overview of high-precision geochronology studies on porphyry deposits. Data for Ok Tedi, Bajo de la Alumbrera and Koloula are from Large et al. (2018), Buret et al. (2016) and Tapster et al. (2016), respectively. Coloured vertical bars are emplacement ages of different intrusive rocks described from these deposits. Intrusions are categorised as pre-ore, pre- and syn-ore, and post-ore intrusion. Deposit size decreases from left to right (tonnages from Cooke et al., 2005). Diagrams are aligned so that the onset of zircon crystallisation overlaps in all deposits. Grey bars indicate the recorded duration of zircon crystallisation. Yellow bars illustrate the maximum durations of total ore formation or individual ore formation pulses. The yellow bar fading out downwards indicates the absence of a post-ore intrusion and the inability to constrain the total duration of ore formation. Note that we excluded the sample X176 from Koloula as it is not related to magmatic history leading to ore formation (Tapster et al., 2016).

5.4 Resolving lower crustal magmatic processes from zircon petrochronology

The lack of inheritance within the zircon record at Batu Hijau suggests that the crustal magmas experienced very minor crustal assimilation. Typically, magmas that are associated with porphyry ore formation contain diverse suites of inherited zircons (e.g. Tapster et al., 2016; Lee et al., 2017; Large et al., 2018), which have been interpreted to represent extended interaction with arc lithologies (Miller et al., 2007). This apparent lack of crustal contamination is consistent with the juvenile isotopic signatures (Pb–Pb, Sm–Nd, Rb–Sr) of intrusions in the Batu Hijau district (Garwin, 2000; Fiorentini and Garwin, 2010). The juvenile and “porphyry-fertile” magmas at Batu Hijau have been explained by asthenospheric mantle upwelling through a tear in the subducting slab that resulted from the collision with the Roo rise (Garwin, 2000; Fiorentini and Garwin, 2010). This would also explain why the only mined porphyry deposit in the Sunda–Banda arc (Batu Hijau) and the most promising prospects (Elang and Tumpangpitu) are located above the inferred margin of the subducting Roo rise (Fig. 1).

The formation of porphyry Cu(–Au) deposits has been commonly associated with the fractionation of amphibole ± garnet in thickened crust (e.g. Rohrlach et al., 2005; Lee and Tang, 2020) within lower crustal magma reservoirs that are active over several million years (Rohrlach et al., 2005). Zircons have been suggested to directly track this extended lower crustal history (Rohrlach et al., 2005). At Batu Hijau no zircon was identified that crystallised in a resolvable way before the main crystallisation period, which we consider to have occurred in the mid- to upper crust (Fig. 5, and discussion above). Non-zoned cores surrounded by oscillatory zoned rims (Fig. 3f) could be interpreted to reflect a two-stage crystallisation process; however, the depth of these two processes cannot be resolved, and they would have occurred within the few hundred thousand years of recorded zircon crystallisation (Fig. 5). As most crystals within a mount are not polished exactly to their centre, the non-zoned cores could equally likely represent a polishing effect where the surface of one zone appears as an non-zoned core. Therefore, it is highly speculative to directly relate zircon textures to a locus or style of zircon crystallisation.

In the case of Batu Hijau, petrochronological data were used to reconstruct the mid- to upper crustal magma evolution, but the data can only provide indirect information about the lower crustal processes involved in the formation of the deposit. For example, the overall elevated Eu ∕ Eu^* of the investigated zircons (0.4–0.7; cf. Loader et al., 2017) could be the result of amphibole fractionation in the lower crust, which would have relatively enriched the residual melt Eu compared to the other REE. This would be analogous to elevated whole-rock Sr ∕ Y ratios in exposed rocks being indicative of the lower crustal fractionating assemblage (Rohrlach et al., 2005; Chiaradia, 2015). The intra-crystal and intra-sample trends of decreasing Eu ∕ Eu^* discussed above describe the evolution within the mid- to upper crustal magma reservoir that was dominated by plagioclase crystallisation and do not reflect any lower crustal process. Zircon can thus directly record the mid- to upper crustal magma evolution, but the information about lower crustal processes is limited to potentially identifying the chemistry of melt and magma that was injected from below into the mid- to upper crust, where zircon started crystallising.

5.5 An assessment of the accuracy and precision of CA–ID–TIMS and in situ U–Pb zircon geochronology

The U–Pb dataset from Batu Hijau allows a critical comparison of the two zircon U–Pb geochronology techniques (LA-ICP-MS, CA–ID–TIMS) that have different analytical precision and can analyse samples on varying spatial scales. Previous investigation of the same lithologies by SHRIMP (Garwin, 2000) allows further comparison. The spatially resolved and fast in situ U–Pb geochronology techniques – LA-ICP-MS or secondary-ion mass spectrometry (SIMS)–SHRIMP – allow the investigation of different crystal domains, whereas the much more time-consuming CA–ID–TIMS analysis of zircons or zircon fragments provides the highest analytical precision. The in situ techniques can discriminate between different zircon populations within single crystals (e.g. inheritance), whereas CA–ID–TIMS geochronology allows for an analytical precision for individual grains that is more than 10 times better, which is required to resolve rapid geochronological events. To increase precision of the in situ techniques, large numbers of individual dates that are considered to represent the same geological event are commonly used to calculate a weighted mean date and standard error of the mean (Wendt and Carl, 1991). On the other hand, the CA–ID–TIMS community has started to measure only small zircon fragments to increase spatial resolution (e.g. Samperton et al., 2015; Smith et al., 2019). Here, the comparison of the different U–Pb zircon techniques applied to the same rock suite allows an assessment of the accuracy of the techniques and of the effect of statistical treatment on the accuracy and precision of the different techniques.

At Batu Hijau, the youngest individual CA–ID–TIMS U–Pb date of each sample is used as our best estimate for the emplacement age of the respective porphyry based on the low number of analyses and high ratio of crystallisation duration to individual uncertainty (Keller et al., 2018). The resulting porphyry emplacement ages are 3.736±0.023, 3.697±0.018 and 3.646±0.022 Myr for the old, intermediate and young tonalite, respectively (Fig. 5). Other statistical means to determine the emplacement ages would be the calculation of a weighted mean for the youngest zircon population (e.g. Buret et al., 2016) or a stochastical sampling approach (Keller et al., 2018). We calculated each alternative and found that the emplacement ages overlap within uncertainty so that the timescales and conclusions derived from this study remain identical (see Appendix). The extended range of concordant zircon dates obtained by CA–ID–TIMS does not allow us to distinguish between different stages of zircon crystallisation within each sample (e.g. inherited vs. autocrystic), but the common geochemical trends support crystallisation within one common magma reservoir (see above). Thus, the range in zircon dates preceding emplacement is interpreted to represent zircon crystallisation within the underlying source magma reservoir over parts or, depending on the onset of zircon saturation, the entirety of its lifetime. The total recorded duration of zircon crystallisation is 336±27 kyr.

Even more so than the CA–ID–TIMS dates, in situ analyses by LA-ICP-MS illustrate an extended range of zircon dates that cannot be separated into different stages of zircon crystallisation. However, the span in zircon dates is about a magnitude higher for the in situ analyses (1.41±0.5–2.1±1.1 Myr) than that obtained by CA–ID–TIMS (0.171±0.026–0.246±0.028 Myr). It could be argued that such LA-ICP-MS data in isolation might indicate an extended period of zircon crystallisation not covered by CA–ID–TIMS data, potentially due to sampling bias or by one sample-set being inaccurate, but we suggest (see Sect. 5.6 below) that the span within the in situ data is the result of analytical scatter and potentially minor amounts of Pb loss or common Pb.

Sampling bias in the selection of the zircons for CA–ID–TIMS geochronology can be excluded as the analyses were conducted on chemically abraded zircons (Mattinson, 2005) that cover the oldest and youngest dates obtained by LA-ICP-MS (Fig. 10). High accuracy of both CA–ID–TIMS and LA-ICP-MS datasets is suggested by routine measurements of secondary standards during the LA-ICP-MS analytical run (see Supplement) and regular measurements of zircon standards by CA–ID–TIMS over the period of data acquisition including the Pleistocene Aus_Z7_5 (von Quadt et al., 2016; Wotzlaw et al., 2017). The distributions of the zircon dates of each sample, as illustrated by probability density plots (Fig. 11), illustrate that the peak of the LA-ICP-MS and SHRIMP dates falls within the mean of zircon crystallisation as defined by the CA–ID–TIMS dataset. This suggests that all datasets are accurate but that the in situ data display more scatter and lower precision. LA-ICP-MS analyses record younger zircon dates for core analyses than rim analyses in 17 of 47 cases; however the LA-ICP-MS dates for core and rim always overlap within uncertainty. Direct comparison of U–Pb dates from the same zircon crystals by the two techniques (Fig. 10) reveals that the less-precise LA-ICP-MS data are not correlated with the more precise TIMS ages, and the suggested dates from the two techniques do not overlap within uncertainty in some cases (6 in 47 cases for rim analyses; Fig. 10b). This could indicate that uncertainties calculated for the LA-ICP-MS data have been underestimated in relation to the achieved precision of the technique. However, due to the high number of analyses it is more likely that it is purely an effect of analytical scatter, where 5 % of the data do not fall within the 95 % confidence interval. This is corroborated by ∼7 % (39/554) of LA-ICP-MS dates not overlapping with the minimum overall duration of zircon crystallisation identified by CA–ID–TIMS dates from all porphyries (336±27 kyr). Additionally, minor amounts of Pb loss or common Pb not identified during data-screening could account for some older and younger outliers, but the overlapping peaks of the in situ populations and the mean of CA–ID–TIMS dates support no systematic bias of the bulk population to younger or older dates. It is therefore concluded that all three techniques are accurate and represent the ∼300–350 kyr of zircon crystallisation. The high number of analyses obtained by LA-ICP-MS together with the lower-precision results in extreme outliers that extend the apparent duration of zircon crystallisation but can be regarded purely as a statistical sampling result that does not indicate a more extended duration of zircon crystallisation. Therefore, we consider alternative methods of estimating age uncertainties in the following section.

Figure 10 Comparison of in situ LA-ICP-MS dates and CA–ID–TIMS dates on the same zircons. In the upper panel each CA–ID–TIMS date is aligned with the rim (filled) and core (empty) LA-ICP-MS date of the same zircon. Coloured bars indicate individual CA–ID–TIMS analysis, including analytical uncertainty (2σ). Downward-pointing black arrows indicate that core analyses are older than rim analyses of the respective zircon, whereas upward-pointing red arrows indicate the opposite. Note that CA–ID–TIMS dates can be plotted several times, with core and rim analyses of the same zircon. The lower panel compares the CA–ID–TIMS date with the respective rim analysis.

5.6 Determining geological ages, uncertainties and rates from in situ U–Pb data

Understanding the timing of magma emplacement, crystallisation or eruption is essential for determining dates and rates of magmatic processes and those directly related to or bracketed by them. Where high-precision CA–ID–TIMS data are not available, porphyry emplacement ages are commonly inferred by calculating a weighted mean and standard error from the youngest overlapping population of in situ U–Pb dates (e.g. Correa et al., 2016; Rezeau et al., 2016; Lee et al., 2017). In the case of Batu Hijau such a calculation would include all LA-ICP-MS zircon dates for each sample as there is no apparent reason to exclude parts of the dataset due to inheritance, common Pb or Pb loss (Fig. 7). The resulting weighted mean dates for the old, intermediate and young tonalite are $\mathrm{3.879}\pm \mathrm{0.027}/\mathrm{0.064}$ (MSWD =2.1, n=207), $\mathrm{3.783}\pm \mathrm{0.023}/\mathrm{0.061}$ (MSWD =2.5, n=189) and $\mathrm{3.751}\pm \mathrm{0.023}/\mathrm{0.060}$ (MSWD =2.6, n=158), where the first stated uncertainty is the standard error including internal uncertainties and those associated with tracer calibration (Schoene, 2014), and the second includes the added 1.5 % external uncertainty as suggested by Horstwood et al. (2016) to account for excess variance. The MSWD for each dataset (2.1–2.6) is elevated in respect to the sample size (n=150–200) based on the formulation by Wendt and Carl (1991). This suggests either an underestimation of the individual uncertainties or that the data do not represent a normal distribution, e.g. by prolonged zircon crystallisation. However, there is no obvious treatment of the data to obtain more appropriate MSWDs. Under these conditions weighted means and standard errors should not be considered geologically meaningful according to Wendt and Carl (1991), but we nevertheless calculate these numbers to illustrate a few points below. Analogously, the weighted mean and standard error of all zircons analysed by SHRIMP from each sample results in weighted means of 3.74±0.14 Myr (MSWD =1.2, n=8), 3.843±0.094 Myr (MSWD =1.2, n=18) and 3.81±0.2 Myr (MSWD =2.35, n=7) for the old, intermediate and young tonalite, respectively (Fig. 11). The weighted means of the different tonalites obtained by LA-ICP-MS would be in accordance with cross-cutting relationships, whereas the SHRIMP dates overlap within uncertainty. The calculated standard errors for the LA-ICP-MS dates are significantly smaller than for the SHRIMP data. The decrease in the standard errors is directly correlated with the increasing sample size (Wendt and Carl, 1991; McLean et al., 2011b), implying that a comparably high number of SHRIMP analyses would result in similarly low standard errors. Irrespective of the different standard errors, the calculated weighted means by SHRIMP and LA-ICP-MS overlap within uncertainty, thus seemingly suggesting that both are accurate or similarly inaccurate (see below).

At Batu Hijau, emplacement ages determined by CA–ID–TIMS geochronology are systematically younger than the weighted mean dates calculated from in situ data (100–150 kyr, except CA–ID–TIMS and SHRIMP for the old tonalite), and the emplacement ages determined by CA–ID–TIMS do not overlap with the LA-ICP-MS mean ages within their attributed uncertainties (Fig. 11). Indeed, disparities between different U–Pb datasets on the same porphyry samples have been noted in several studies comparing high-precision CA–ID–TIMS data with in situ data (von Quadt et al., 2011; Chiaradia et al., 2013, 2014; Chelle-Michou et al., 2014; Correa et al., 2016). As there is no evidence for any systematic errors jeopardising the accuracy of any of the three methods in this study, all methods are considered accurate, and the reason for apparent discrepancy must lie in the geological interpretation of the statistical uncertainties.

Figure 11 Probability density plots for the geochronology data for each analytical technique. All dates of each technique are combined in panel (a). Plots in panels (b), (c), and (d) are constructed from the data of the old, intermediate and young tonalite. Dashed lines indicate the youngest and oldest zircon crystallisation age as determined by CA–ID–TIMS for the respective investigated dataset. Weighted means, standard error, standard error + 1.5 % (Horstwood et al., 2016) and standard deviation are shown by bars with three colours to the right of each diagram. These are identical to those in Fig. 5.

The protracted zircon crystallisation identified at Batu Hijau has broader implications for the determination of magma emplacement, crystallisation or eruption ages. Extended magma reservoir lifetimes are not unique to Batu Hijau but a commonly described feature (e.g. Miller et al., 2007; Claiborne et al., 2010a; Reid et al., 2011; Buret et al., 2016). A weighted mean is a measure to quantify the mean of a population while emphasising the importance of values with low uncertainties over those with high uncertainties (Reiners et al., 2017) and is only valid in cases where the data are normally distributed around one expected value (Wendt and Carl, 1991). This is rarely the case when investigating geological processes, and indeed the presented in situ datasets record protracted zircon crystallisation (>300 kyr) in the magma reservoir that results in zircon population distributions that cannot be easily defined statistically (e.g. by a normal distribution; Figs. 5, 7; cf. Keller et al., 2018). The calculated weighted mean thus does not represent the porphyry emplacement age but the average zircon crystallisation. This is corroborated by the weighted means of the LA-ICP-MS and SHRIMP dates approximately describing the mean of the zircon populations defined by CA–ID–TIMS (Fig. 11). Therefore, the calculated weighted mean dates of datasets with a large number of dates but low precision do not necessarily describe any specific geological event, especially as the uncertainties indicated by the standard error for the LA-ICP-MS data are too small to even cover the entire recorded duration of zircon crystallisation.

Traditionally, problems related to the oversimplification associated with calculating weighted means and their standard errors in geochronology were hidden by the larger uncertainties resulting from larger analytical uncertainties and smaller sample sizes. Due to rapid data acquisition nowadays by in situ techniques, calculated standard errors can result in uncertainty envelopes of <0.1 % for a sample. In the case of the LA-ICP-MS dates from the Pliocene Batu Hijau porphyry intrusions, the standard error of the weighted mean (∼1 %, ∼40 kyr) is on the same order of magnitude as an individual CA–ID–TIMS date and therefore smaller than the geological spread of zircon crystallisation ages. The combination of using a weighted mean to describe a non-Gaussian sample distribution with the very small uncertainties attributed to a weighted mean results in highly precise dates that may have no relation to a specific geological event.

The MSWD (a reduced chi-square statistic) of a dataset provides a first measure to indicate whether dates are normally distributed around an expected value and thus whether the calculated weighted mean and standard error are of significance ($\surd (\frac{\mathrm{2}}{n-\mathrm{1}}$) rule by Wendt and Carl, 1991). As discussed before, the MSWDs for the LA-ICP-MS data are elevated, implying that weighted means and standard errors should not be calculated to characterise a geological event. However, the MSWDs for the SHRIMP zircon analyses of the intermediate and old tonalite are within the 95 % confidence interval of the MSWD, mainly due to the larger individual uncertainties. Still, they are similarly affected by protracted zircon crystallisation, which biases the weighted mean to higher values (Fig. 11). Furthermore, overestimated individual uncertainties can result in acceptable MSWDs but similarly inaccurate dates and low standard errors. For example, increasing individual uncertainties for the ²⁰⁶Pb–²³⁸U dates obtained for the old tonalite by LA-ICP-MS by a factor of 1.5 would result in an acceptable MSWD (0.95), but the weighted mean and standard error would be nearly identically precise to those calculated with the actual uncertainties but inaccurate (3.880±0.041, 2SE). Based on the presented data it is advised not to characterise a geological event by a weighted mean with an associated standard error if the MSWD is elevated (Wendt and Carl, 1991).

Based on the presented data we recommend using an uncertainty attributed to the weighted mean that is more representative of the uncertainty of the individual analyses so that it will most likely cover the event actually dated. Here, we tested the standard deviation of zircon dates from each sample as a measure for the uncertainty of the weighted mean. This approach provides a more realistic estimation of the uncertainty associated with calculating a weighted mean of a geologically young dataset as it describes the variability in the measurements, (0.29–0.62 Ma, Figs. 7, 11), and, importantly, it would be independent of the number of analyses. The resulting values at least for the Pliocene Batu Hijau deposit result in appropriate uncertainties for the weighted mean as it would cover an appreciable part of the range of in situ dates and thus the >300 kyr of zircon crystallisation and the emplacement age.

The presented data highlights the importance of CA–ID–TIMS zircon U–Pb geochronology to resolve complex zircon crystallisation patterns. In the case of porphyry research, high-precision ID-TIMS dates are required to resolve the durations of porphyry emplacement and hydrothermal processes, but in situ data can reliably reconstruct a timeline of magma emplacement events within porphyry districts over million-year timescales (e.g. Rezeau et al., 2016). Furthermore, combination with in situ petrochronology techniques (i.e. U–Pb isotope and geochemical data from the same analyte) allows us to screen zircons for inheritance and more importantly provides spatially resolved geochemical information that can be integrated with high-precision dates.

6 Conclusions

High-precision zircon geochronology by CA–ID–TIMS combined with in situ zircon geochemistry provides valuable datasets that allow the reconstruction of geological processes with the highest temporal resolution. At Batu Hijau zircons record the magmatic-to-hydrothermal evolution of the world-class Batu Hijau porphyry Cu–Au deposit from the onset of zircon crystallisation to emplacement of the post-ore young tonalite. The magma reservoir that sourced the tonalites and the Cu–Au mineralising fluids records zircon crystallisation over 336±27 kyr. Emplacement of the first exposed tonalite at the Batu Hijau deposit (old tonalite) occurred after 246±28 kyr of uninterrupted zircon crystallisation in this subjacent reservoir. Zircon trace element signatures support a dominantly crystallising and cooling magma reservoir over 285±24 kyr until emplacement of the intermediate tonalite. After emplacement of the intermediate tonalite the chemistry of the reservoir remained in rather steady conditions for 51±28 kyr, during which it could have been disturbed by magmatic recharge or underplating until final emplacement of the young tonalite. Ore formation is most probably associated with the last stages of the chemically and thermally evolving magma reservoir. The maximum duration of ore formation can be constrained to <122 kyr by the emplacement ages of pre-to syn-ore old tonalite and the post-ore young tonalite. This maximum duration of ore formation covers different pulses of mineralisation that could have lasted only a few thousand years. We record a magmatic system that was active over ∼250 kyr before emplacement of the first porphyry intrusion and onset of several pulses of hydrothermal activity forming the world-class ore reserve in less than 100 kyr.

Comparison between in situ LA-ICP-MS and SHRIMP as well as CA–ID–TIMS U–Pb geochronology reveals that all techniques provide accurate individual dates (within the stated confidence interval). However, statistical treatment of in situ data by calculating a weighted mean and standard error can result in highly precise but inaccurate older ages of questionable geological significance with apparent uncertainties that do not provide an accurate measure for the uncertainty of emplacement age. The tempo of magma evolution and hydrothermal processes associated with magmatic–hydrothermal systems, such as porphyry deposits, is too fast to be reliably resolved by currently available in situ U–Pb geochronology and requires ID-TIMS geochronology. Combination of high-precision geochronology with in situ or TIMS–TEA (thermal ionization mass spectrometry with trace element analysis) geochemistry is currently the most powerful tool in deciphering these geologically rapid processes.