Zircon from Altered Monzonite Rocks Provides Insights into Magmatic and Mineralizing Processes at the Douay Au Project, Abitibi Greenstone Belt

Abstract

Zircon provides essential information on the age and oxidation state of magmatic systems and can be used to characterize magmatic-hydrothermal Au mineralizing systems. Using the Douay intrusion-related gold system (IRGS) as a type example of Neoarchean syenite-associated mineralization (Abitibi greenstone belt), we demonstrate that zircon from altered quartz-monzonite rocks can also be used to infer the age of a magmatic-hydrothermal event. Here, zircon chemistry is used to identify the following sequence of events at the Douay exploration project: (1) the crystallization of zircon at ~2690 Ma in evolved residual melts with distinct U-contents (quartz-monzonite magma); (2) the extensive radiation damage for the U-rich grains over a period of ~10–15 My; and (3) the alteration of zircon grains at ~2676 Ma by interaction with magmatic-hydrothermal mineralizing fluids derived from syenite and carbonatite intrusive phases. This study also distinguishes extensively altered zircon grains from pristine to least-altered zircon formed in distinct magmatic environments using a Th/U vs. U discrimination diagram.

1. Introduction

Magmatic-hydrothermal mineralizing systems are now frequently recognized in the Abitibi and other greenstone belts [1,2,3,4]. This includes early porphyry Au-(Cu) and Cu-Au deposits, e.g., the Central Camp and Côté Gold deposits [4,5,6,7], as well as Au-dominated systems formed prior and during the main deformation phase (syntectonic period) that led to the stabilization of the craton. These latter systems correspond to porphyries, e.g., Bachelor deposit [8], and intrusions-related gold systems (IRGS), also referred to as syenite-associated Au deposits [1]. In the gold-endowed Abitibi greenstone belt, Au was also concentrated by orogenic gold systems (OGS), which correspond to hydrothermal systems dominated by metamorphic fluids [9,10]. Additionally, OGS tend to postdate magmatic-hydrothermal systems [1]. Most of the gold endowment of the Abitibi greenstone belt is associated to OGS [3], even if the involvement and importance of magmatic fluids remain debated for many deposits [2].

Distinguishing IRGS from OGS and hybrid (e.g., IRGS overprinted by an OGS) systems is essential for exploration strategies, as well as for our understanding of fertilizing processes and of the distribution of gold in the upper Neoarchean crust. Moreover, the mineralizing style is an important parameter. For example, mineralization at the Douay exploration project displays characteristics generally attributed to IRGS [1], such as a polymetallic signature, disseminated pyrite, and spatial association to a syenite intrusion and to high temperature alteration such as K-metasomatism [11]. However, the mineralization is located next to a large-scale fault and part of the mineralization is structurally controlled [11], with these latter characteristics being more generally attributed to OGS [9]. The Douay exploration project is a representative example of a gold mineralization of the syntectonic period that may correspond to an IRGS, an OGS, or a hybrid system [2,3,11].

This contribution aims to determine the timing of magmatic and hydrothermal events at the Douay project using the chemistry of zircon from an intrusive phase of the Douay intrusion. Zircon chemistry has been used to constrain the oxidation state of magmatic systems and their gold potential [12,13,14], while evaluating the maximum amount of Au that can be dissolved by a magma requires insights into the oxygen fugacity parameter (fO₂) and volatile (S, Cl) content of the melt [15]. Zircon, in a mineralized alkaline system, can be U-rich, metamict, and extensively altered, as is the case here. This study demonstrates that altered zircon grains can provide insights into the crystallization age and oxidation state of a magma, and that additional minerals (apatite, amphibole) are required to fully characterize the volatile content. Zircon chemistry, however, can be used to discuss the characteristics and timing of magmatic and hydrothermal events. Zircon is thus used here to gain insights into these processes at the Douay exploration project.

2. Geological Setting

Alkaline-carbonatite intrusive complexes, including the Douay intrusion, are recognized in the northern part of the Abitibi greenstone belt, which is located in the southern part of the Superior craton (Figure 1). As is common in many other greenstone belts [16], the tectono-magmatic evolution of the Abitibi belt can be summarized in two main stages, i.e., the synvolcanic and syntectonic stages [17,18]. In the northern part of the Abitibi greenstone belt, the synvolcanic stage initiates at or before ca. 2.80 Ga, as ~2.79 Ga mafic and felsic volcanic units [19] and ~2.80 Ga tonalite intrusion [14] emplace in the upper crust. The synvolcanic stage was followed by the syntectonic stage, which extends from ca. 2.70 to 2.69 Ga in the northern part of the Abitibi greenstone belt [1,20,21]. The syntectonic stage follows, or initiates, during an early compressional event (D₁) that is followed by the development of Porcupine basins (flysch-like sedimentary sequences) and Timiskaming basins (conglomerate-bearing) [1]. Coeval with basin development during an extensional stage, alkaline magmatism and the related syenite-associated (or IRGS) style of mineralization then developed in the upper crust [1,2]. Following basin development, the main N–S compressional event (D₂) and its waning stage (D₃) produced the main E-W faults, foliation planes, and folds observed in the field, and shortening is coeval with the peak of greenschist facies metamorphism [1,22].

The Douay intrusion is located in the northwestern part of the Abitibi greenstone belt, 55 km SW of the town of Matagami. The Douay area is an exploration project owned by Maple Gold Mines Ltd. and Agnico Eagle Mines Ltd. The gold mineralization is mostly contained in the Douay intrusion and its host rocks, i.e., the Cartwright Hills and Taibi groups (Figure 1). A syenite phase of the Douay intrusive complex is dated at 2676 +5/−6 Ma (upper intercept ²⁰⁷Pb/²⁰⁶Pb age, obtained by regression of three titanite grains and one zircon fragment) [23]. The Douay complex mostly consists of aphyric to feldspar-phyric and quartz-feldspar-phyric syenite, as well as aplite, pegmatite, and carbonatite dykes [11]. The Cartwright Hills Group (2721.4 ± 3.1 Ma; [24]) is dominated by tholeiitic pillow basalts and gabbro sills, and also comprises ultramafic lava flows, as well as interflow felsic pyroclastic and sedimentary rocks such as chert and siltstone [11,24]. The Taibi Group (<2685 ± 2.6 Ma; [25]) is a basin-restricted sedimentary unit formed in deep water toward the end of the synvolcanic period [3]. The sedimentary rocks consist of siltstone, mudstone, graywacke, polymict conglomerates, and iron formations, and are interbedded with mafic to felsic lava flows and pyroclastic rocks.

The Douay gold project is spatially associated to the Casa Berardi Tectonic Zone (CBTZ), which overlaps the Cartwright Group to the south and the Taibi Group to the north [11]. The CBTZ is an E–W trending fault that accommodated brittle-ductile movement evolving from reverse (during the D₂ deformation stage) toward dextral-transpressive motion (D₃ stage) [11,26]. Mineralization at the Douay project comprises disseminated and structurally controlled styles, hydrothermal breccia, and hydrothermal alteration induced K-metasomatism, albitization, hematization, carbonatization, silicification (quartz-carbonate veins included), and fluoritization [11].

3. Materials and Methods

3.1. Petrography and Whole Rock Chemistry

Only a single outcrop of the alkaline intrusive complex can be found on the Douay property. Five samples were taken from this outcrop in August 2020 (

Table 1. Sample location (WGS84 datum) and general description.

Sample	Latitude	Longitude	Description
D1A	49.531720	−78.142188	coarse-grained intrusive phase (carbonatized)
D2	49.531588	−78.142307	coarse-grained intrusive phase (carbonatized)
D3	49.531649	−78.142220	aplite dyke
D4	49.531677	−78.142260	strongly altered intrusive phase (‘red rock’)
D5	49.531620	−78.142139	coarse-grained feldspar-phyric intrusive phase (carbonatized)

). Thick (50 μm) polished sections were obtained from these five samples. Petrographic observations were performed using an Olympus BX53M petrographic microscope equipped with a motorized stage and a SC50 camera. Further characterization was performed using a petrographic microscope (Olympus Instruments, Olympus Corporation, Tokyo, Japan) and cathodoluminescence (CL8200 Mk5-1 Optical, CITL, Hertfordshire, United Kingdom) at UQAC.

Whole rock analyses were performed by ALS Canada Ltd. on the samples. The major elements were analyzed by inductively coupled plasma (ICP) atomic emission spectroscopy (AES) following acid digestion and the detection limit is 0.01%. Trace elements were analyzed by ICP mass spectrometry (MS) using lithium borate fusion and the detection limits are between 0.01 and 5 ppm (Supplementary Material S1). Quality assurance and quality control procedures, including the use of blanks, duplicates, and standards, were maintained rigorously to obtain precise and accurate results.

The chemistry of the Douay intrusion is compared to the chemistry of another alkaline complex of the Abitibi greenstone belt, i.e., the Beattie syenite. The Beattie syenite is selected because its chemistry is well documented and because this syenite intrusion is genetically related to an IRGS [27]. The chemical composition of the Beattie syenite (n = 32) is compiled from the literature [27,28,29].

3.2. Zircon Chemistry and Geochronology

Zircon concentrates were produced for two samples (D1A and D3) by the Overburden Drilling Management Ltd. (ODM) laboratory (Ottawa, ON, Canada). At ODM Ltd., the rocks were crushed using the electric pulse disaggregation method. Heavy minerals were separated using a shaking table and heavy liquids with densities of 3.0 and 3.3 g·cm⁻³. Pan concentrates containing the zircon fraction were poured on a tap, mounted onto epoxy resin and fine-polished at UQAC. In both samples, zircon grains have distinct sizes and the largest (heaviest) fraction separated gravitationally from the smaller zircon grains during the preparation of epoxy mounts. The large grains (~70–200 µm) and small grains (<100 µm) were mounted separately. The mounts were imaged in transmitted light using an Olympus BX53M microscope at UQAC and by backscattered electron (BSE) and cathodoluminescence (CL) imaging using a Tescan Vega 3 SEM at the MERC-Isotope Geochemistry Laboratory (MERC-IGL) of Laurentian University. The images were used to identify altered and fractured grains, as well as apatite and other types of inclusions, which are avoided during subsequent laser ablation analysis.

Zircon U-Th-Pb isotopes and trace element analyses were conducted on ~300 zircons per sample at MERC-IGL, Laurentian University. Laser ablation sampling was performed using a Photon Machines Analyte G2 ArF excimer laser, with 193 nm wavelength, <4 ns pulse width, and HelEx II cell. Ablation durations were 30 s, with a laser fluence of 2 J/cm² and 7 Hz repetition rate, leaving estimated ablation pit depths of <15 μm. Sixty seconds of background were measured at the beginning and end of the analytical session, with 30 s of background measured between each ablation. Helium carrier gas flows through the ablation cell were 0.525–0.55 L/min (cup) and 0.1 L/min (cell), with 1.68 L/min Ar and 0.015 mL/min N₂ makeup gas added downstream of the cell. The ablated aerosol was split downstream of the sample cell, so that U-Th-Pb and trace elements measurements could be conducted simultaneously (LASS; e.g., [30]).

The U-Th-Pb isotopes measurements were conducted using a Thermo Scientific Neptune Plus multicollector (MC) ICP-MS, equipped with a Jet interface and nine Faraday cups. All MC-ICP-MS analyses were conducted in a low resolution, static mode to ensure maximum sensitivity and stability. All cups were coupled with 10¹¹ Ω amplifiers, with the exception of L2, which was coupled with a 10¹² Ω amplifier for better precision on the low intensity ²⁰⁷Pb signal. Ion beam drift is typically <0.05 amu, with a flat top peak width typically ~0.2 amu, ensuring stable on-peak measurements throughout long-duration analytical sessions. Cool gas and Ar auxiliary gas flows were set to 0.07 L/min, respectively, and RF power was set to 1200 W. Trace element measurements were conducted using a Thermo Scientific iCap-TQ ICP-MS in single quad mode to ensure maximum sensitivity on the low to intermediate mass range. The masses measured and associated dwell times are shown Table S4.

The raw U-Th-Pb and trace element data were synchronized in Iolite, along with the laser log file, prior to processing, so that the same timeslices (selections) for each analysis were used in the calculations. Baseline subtraction, instrumental drift, and downhole fractionation corrections were performed with the U-Pb Geochronology reduction scheme (DRS) implemented within Iolite v4 [31,32], with U-Th-Pb isotope ratios normalized to the zircon reference material OG1 (OG1; 3465 ± 1 Ma; [33]). The OG1 standard was analyzed three times at the beginning and end of the session, and once every ten unknowns throughout the session. Typically, 3 s at the beginning and 1 s at end of the ablation period were excluded from the selections in order to minimize potential fractionation effects, leaving ~26 s of signal for integration. The within run variance in the measured ratios for OG1 (i.e., the additional percent error required to achieve MSWD = 1) was propagated into the 2 standard errors (2SE) uncertainty for all unknowns.

Verification reference materials DD91-1/DM18 (2682 ± 1 Ma; [34]) and Maniitsoq (3008.7 ± 0.66 Ma; [35]) were analyzed repeatedly throughout the session to ensure accuracy of the U-Pb ratios, and measured dates identical to the published TIMS ages and low uncertainty and MSWD were obtained. Trace element data were processed using the Internal Standard DRS method within Iolite v4 and normalized to the synthetic glass NIST610. An assumed stoichiometric concentration of 15.284 wt% Si was used for all zircon, to account for the differing ablation characteristics between the glass and zircon. The synthetic glasses NIST610 and NIST612, and zircon 91500 (1065 Ma; [36,37]), were used to verify the accuracy of the trace elements analyses. Analytical results and operating parameters are provided in Supplementary Material Tables S1–S4.

4. Results

4.1. Fieldwork

The Douay alkaline-carbonatite intrusive complex is covered in glacial material and is mostly known from drill core. A significant volume of rock was taken from the single outcrop that exposes the alkaline intrusive phases in the study area. This outcrop is located in the Porphyry zone delineated by exploration (Figure 1) [11]. This natural outcrop, enlarged by stripping, mostly exposes a coarse-grained equigranular to feldspar-phyric alkaline intrusive phase that is pervasively hematized and carbonatized. This intrusive phase comprises red patches of strongly altered rock (‘red rock’) that are elongated in the E–W to ESE–WNW direction over 1–5 m. Thin and elongated silicified zones are also observed (Figure 2a). The main intrusive phase is cut by several 10 cm to 1 m thick aplite dykes, oriented NW–SE (Figure 2b). Samples were taken from the main intrusive phase, including altered (samples D1A, D2, and D5; Figure 2c,d,g) and strongly altered (D4; Figure 2f) samples. A sample was also taken from the thickest aplite dyke observed in the field (D3; Figure 2e).

4.2. Petrography

Samples D1A, D2, and D5 display similar mineral assemblages. They are dominated by euhedral, 0.5–2 mm long, feldspars (aspect ratio = 0.3–0.6) with irregular, possibly albitized, contacts (Figure 3a,b). Feldspar displays coarsened or irregular (chessboard texture) perthite exsolutions (Figure 3c,d,e,h). Feldspar is dark to orange-pink, or green in CL, possibly due to irregularly distributed Na and/or K-metasomatism, and hematization. The D1A, D2, and D5 samples also contain a minor amount (2–5 vol%) of ~50 µm rounded quartz included in feldspar and, less frequently, occurring as small (˂10 µm) interstitial grains. Samples D1A, D2, and D5 also contain clusters of strongly altered mafic minerals that occupy 5–10 vol% of the rock (Figure 3a,b). These clusters consist of biotite, chlorite, carbonate, epidote, and/or feldspar. According to CL imaging, apatite grains are most abundant in these mafic clusters. Small (˂5–10 µm) disseminated magnetite (˂1 vol%) is also observed, which is mostly restricted to the clusters of mafic minerals. The magnetite grains are rounded or anhedral with irregular outlines.

Quartz and carbonate veinlets (up to 0.5 cm thick in sample D1A) occupy ˂2 vol% of the rock. Carbonate stringers are most abundant in the broken-up rock located next to the 0.5 cm thick carbonate vein observed in sample D1A (Figure 3c). Shown by CL imaging, carbonate is found in veinlets and is most abundant in mafic minerals clusters. Carbonate is also found disseminated in the rock, within feldspar and at the contact between feldspars (Figure 3e). Trace amounts of minute chalcopyrite and small pyrite grains are also disseminated in samples D1A, D2, and D5.

Sample D4 is dominated by feldspars and contains 10 vol% of clusters of mafic minerals and trace amount of disseminated pyrite. Most feldspars are broken-up to thinly brecciated in sample D4, and fractures are filled with fine-grained phyllosilicates, i.e., white mica and possibly chlorite. This sample also contains 1–2 vol% of disseminated carbonates (Figure 3g). The aplite sample (D3) is fine-grained, equigranular, and dominated by feldspar and ~10 vol% rounded quartz grains (Figure 3b,f). This sample also contains 2 vol% of mafic clusters (magnetite included), 1 vol% apatite, and 1 vol% of disseminated carbonates.

4.3. Whole Rock Chemistry

The four samples from the main intrusive phase (samples D1A, D2, D4, and D5) came from a single outcrop and display similar chemistry (Figure 4). Compared to these samples, the D3 sample (aplite) is silica-enriched and trace elements-depleted (Figure 4). All samples contain a moderate amount of water, i.e., 0.1 to 0.36 wt% H₂O (Table S1). According to the location of samples on the Total Alkali Silica (TAS) diagram, the main intrusive phase is a syenite and the aplite dyke is a granite (Figure 4a)—this preliminary result should, however, be interpreted with care as the TAS diagram is inadequate for altered rocks. On the aluminium saturation index (ASI) vs. Al/NK diagram, the bulk of samples plot between the fields for peralkaline and metaluminous rocks (Figure 4b).

The REE and multielement diagrams confirm the chemical similarities between the studied samples (Figure 4c,d). The REE profiles are fractionated and lack an Eu anomaly, as is generally the case for alkaline intrusions of the Abitibi greenstone belt (Figure 4c). Samples from the Douay project also display the negative Nb, Ta, Zr-Hf, and Ti anomalies observed in other alkaline intrusions such as, for example, the Beattie syenite (Figure 4d) and in tonalite-dominated synvolcanic intrusions [38]. The main difference between the intrusive phase of the Douay intrusion described here and syenite intrusions (e.g., the Beattie syenite) is the general enrichment in trace elements, with the Douay intrusion being depleted in most elements (Nb, Ta, and Ti excepted) compared to syenite intrusions (Figure 4c,d). The aplite dyke is the most depleted in trace elements and displays more pronounced negative Nb, Ta, Zr-Hf, and Ti anomalies compared to the main intrusive phase (Figure 4d).

The method of Trépanier et al. [39] was applied to the studied rocks to calculate pre-alteration major elements composition (i.e., precursor composition). This method is based on an artificial neural network that relies on Al, Ti, Zr, Cr, Y, Nb, and Th to predict precursor composition (i.e., major elements concentrations) for igneous rocks. The method of Trépanier et al. [39] also performs mass balance calculations (

Table 2. Results of mass balance calculation provided in g per 100 g of precursor (absolute mass changes).

Sample	CaO	FeOT	K2O	MgO	Na2O	SiO2
D1A	0.84 1	−0.53	1.31	−0.89	2.12	−4.71
D2	0.14	−0.74	0.58	−0.74	1.11	−10.85
D3	calculation failed (invalid Al/Ti ratio)
D4	0.56	−0.77	−0.07	−0.63	2.80	−15.17
D5	1.13	−0.21	0.60	−0.60	1.28	−11.51

¹ Mass changes are considered negligible when they are between −1 and +1 g for most elements, and between −10 and +10 for silica [39,40].

). This method failed for sample D3 because the Al/Ti ratio of this sample falls outside of the range of values used to train the neural network [39]. For the main intrusive phase, results indicate negligeable mass change for most elements. Exceptions are samples D1A and D5 that gained a small amount of K and Ca, respectively. More significantly, the four samples gained a moderate amount of Na and lost a small amount of Si. The main alterations are thus Na-metasomatism, moderate silicification, as well as K-metasomatism (samples D1A) and, possibly, calcite veining (D5).

The mass changes reported here are moderate, but these mass changes are, nonetheless, sufficient to modify the chemistry of the rocks and to change the location of the samples on the TAS diagram, i.e., the main intrusive phase is likely a quartz-monzonite (Figure 4a). To complete the evaluation of hydrothermal alteration processes, carbonatization is discussed. Most samples are enriched in CO₂ (3.08 to 6.96 wt% CO₂), the D3 sample excepted (0.04 wt% CO₂). The carbonate saturation index (CO₂/[FeO + MgO + MnO + CaO]·100 molar) returns values of 7 for the D3 sample, and of 72–97 for the other samples, pointing to weak carbonatization in the aplite dyke and to intense carbonatization in the main intrusive phase.

Figure 4. Chemical composition of five rock samples from the Douay intrusive complex. Rocks are displayed on (a) the Total Alkali Silicate (TAS) binary diagram and on the (b) Al saturation index (ASI) vs. Al/NK binary diagram [41]. Rocks are also displayed on (c) a primitive mantle-normalized [42] REE plot and (d) a primitive mantle-normalized multi-element arachnid diagram modified from Pearce [43]. The chemical composition of the Beattie syenite (n = 32) is compiled from the literature [27,28,29].

Figure 4. Chemical composition of five rock samples from the Douay intrusive complex. Rocks are displayed on (a) the Total Alkali Silicate (TAS) binary diagram and on the (b) Al saturation index (ASI) vs. Al/NK binary diagram [41]. Rocks are also displayed on (c) a primitive mantle-normalized [42] REE plot and (d) a primitive mantle-normalized multi-element arachnid diagram modified from Pearce [43]. The chemical composition of the Beattie syenite (n = 32) is compiled from the literature [27,28,29].

Whole rock chemical analyses are also used to estimate the zircon saturation temperature and the activity of TiO₂ (a_TiO2), which will be required to apply the Ti-in-zircon thermometer (see below). The solubility of Zr is evaluated using the method of Watson and Harrison [44], which return saturation temperatures of 570 °C (sample D3) and 670–720 °C (other samples). A more recent method was developed by Borisov and Aranovich [45]. This method relies on a large number of compositional parameters, as well as on experiments performed for melts containing >155 ppm Zr. The studied rocks, however, contain <155 ppm Zr, i.e., 8 ppm Zr (sample D3) and 104–152 ppm Zr (other samples). For this reason, unrealistically low temperatures of 340 °C (sample D3) and 413–512 °C (other samples) are obtained, indicating that the method of Borisov and Aranovich [45] cannot be extrapolated to Zr-poor melts.

The a_TiO2 parameter was estimated using the method of Hayden and Watson [46]. Using a generic temperature of 600 °C (sample D3) and 800 °C (other samples), values of 0.5 (D3) and 0.66–0.79 (other) are obtained. Repeating the calculation with the temperatures calculated using the method of Watson and Harrison [44] (see above), values of 0.82 (D3) and 1.86–3.96 (other samples) are obtained for the a_TiO2 parameter. The studied rocks are Zr-poor and zircon likely crystallized after Ti-minerals and from Ti-poor residual melts. Hence, at the time of zircon crystallization, the a_TiO2 parameter was likely closer to the lower values obtained here, such as 0.5–0.66. Typical a_TiO2 values are 0.6 or more in silicic melts at typical magmatic temperatures [46], and the value 0.6 will be used to calculate temperature using the Ti-in-zircon thermometer (see below).

4.4. Texture of Zircon Grains and Ages

Zircon separates were obtained for the D1A and D3 samples. In the D1A sample, broken zircon grains and whole grains are 50–150 μm long (Figure 5). Incomplete zircons (fragments) are most frequent (aspect ratio = 0.5 to 0.7). Complete zircon grains are also found and are elongated (aspect ratio = 0.3 to 0.4). Most grains are subhedral and have {110} (prismatic form) and {101} (pyramidal crystal form) as their dominant form, which is characteristic of zircon crystallized in water-rich granites and pegmatites [47,48]. In sample D1A, most grains are dark in CL and can be porous (Figure 5). A minority of grains (˂2%) display crude oscillatory zoning and irregular zoning. An additional 10% of grains display relics of poorly preserved oscillatory zoning, with a bright core region that can be preserved in some grains. Other grains are homogeneously dark in CL, with more or less expressed porosity (Figure 5). Bright CL irregular and discontinuous rims formed by coupled dissolution-reprecipitation processes are poorly developed. Zircon also contains bright CL inclusions of apatite and of minute secondary minerals (possibly thorite or coffinite) located along cracks.

In the D3 sample, zircon grains are incomplete fragments that are 60–180 μm long, with an aspect ratio of 0.5 to 1 (Figure 6). These grains display thin oscillatory zoning, with irregular zoning (up to thick transgressive zones) being more developed in some grains. Core regions can only be distinguished in 2% of the grains (Figure 6). Cracks can be observed that contain minute bright CL inclusions (possibly thorite). Bright CL outer regions formed by coupled dissolution-reprecipitation processes are poorly developed. Apatite inclusions are small (1–2 μm) and observed in a minority (10–20%) of grains.

Geochronological data for the D1A sample show significant common lead enrichment (Figure 7a), indicating that these zircon grains behave as completely open systems. The addition and/or loss of U and/or Pb during alteration shifted data location in the Tera–Wassserburg plot (Figure 7a) and the trend defined by the dataset has no age significance. As a consequence, ²⁰⁷Pb/²⁰⁶Pb data span a range of values (Figure 7b) and no crystallization age data can be extracted from this sample. For the D3 sample, grains are mostly concordant and analyzed grains can be enriched in Sr and La (Figure 7c,d). Applying a filter criterion for pristine grains (i.e., Sr ˂ 1 ppm), an intercept age of 2690.1 ± 0.96 Ma (MSDW = 0.86) is obtained (Figure 7e,f).

4.5. Chemistry of Zircon Grains

Zircon grains from the D1A and D3 samples display distinct chemical characteristics. In sample D1A, zircon is enriched in trace elements, including Sr and REE (Eu excepted), compared to sample D3 (Figure 8,

Table 3. General chemical composition of zircon grains (median values).

	Sample D3		Sample D1A
	˂3 ppm Sr 1	˃3 ppm Sr	˂32 ppm Sr	˃32 ppm Sr
P (ppm)	451.88	670.62	1210.88	1548.60
Ca	not detected	2274.88	4958.04	6298.40
Sc	220.17	227.61	529.38	1019.42
Ti	2.57	5.93	12.49	19.89
Fe	114.02	235.30	580.58	1110.12
Sr	0.83	9.84	19.57	56.73
Y	228.16	591.78	2091.35	2627.67
Zr	52,0568	508,821	434,006	433,971
Nb	0.88	1.46	153.38	221.65
La	0.20	1.88	7.28	18.22
Ce	22.11	54.18	78.03	167.98
Pr	0.47	4.92	6.14	18.34
Nd	4.26	31.16	25.88	82.25
Sm	3.93	20.58	25.45	68.28
Eu	1.83	8.19	0.05	0.07
Gd	9.59	37.62	59.02	121.49
Tb	2.23	8.11	28.57	52.14
Dy	21.52	68.82	354.73	594.33
Ho	7.00	19.43	123.94	185.67
Er	32.27	80.66	682.79	949.87
Tm	6.70	14.92	183.55	248.54
Yb	55.53	115.15	1973.78	2760.02
Lu	8.68	17.46	403.29	561.33
Hf	5811.33	5457.94	27,517.40	29,958.15
Th	42.29	133.44	492.80	1606.59
Ta	0.09	0.13	134.90	190.82
U	110.07	209.28	5388.84	8383.87

¹ 3 ppm Sr (D3 sample) and 32 ppm Sr (D1A sample) correspond to median Sr values for these samples.

). The Sr and LREE enrichment, for zircon of sample D1A, point to extensive alteration, which obscures the positive Ce anomaly (Figure 8a) and questions the capacity of such zircon grains to provide reliable estimates of the fO₂ parameter (see below). Zircon of sample D1A also display a pronounced negative Eu anomaly, while such an anomaly is lacking in sample D3. In both samples, the range of HREE values is narrower than the range of LREE values. Additionally, zircon of sample D3 display a positive Ce anomaly, which is attenuated in the LREE-enriched grains (Figure 8b).

The chemistry of zircon is further discussed using binary diagrams (Figure 9). The stochiometric Zr value for zircon is 431,400 ppm Zr when using an internal standard value of 147,600 ppm Si (webmineral.com, accessed on January 2022. For sample D1A, the Zr-content of zircon is close to expected stochiometric values, while zircon of sample D3 displays a higher than expected Zr/Si ratio (Figure 9a). The Zr values should, however, be interpreted with care as the accuracy and reproducibility of Zr concentration is not quantitatively evaluated (Zr is analyzed to monitor for inclusions). Zircon, in the D1A sample and in part of the D3 sample is enriched in non-stochiometric elements, and only 3% and 45% of grains in, respectively, the D1A and D3 samples, contain ˂10 ppm Fe + Ca + Sr + La and correspond to least altered zircon (Figure 9a).

In sample D3, Th and Lu are positively correlated (Figure 9b), suggesting co-incorporation in zircon during crystallization. This is possibly the result of a simple substitution mechanism involving the substitution of Zr for trace elements such as Th, Hf [49], and possibly Lu, with Lu being representative of all REE and some other trace elements [50]. In sample D1A, Lu and Th are not correlated for most grains (Figure 9b).

The «xenotime» substitution (Zr + Si substituted for REE + P) [50,51] is also investigated (Figure 9c). A correlation between P and LREE is observed for the D1A sample. For the D3 sample, this correlation is investigated for the grains with a P content that is above the detection limit (i.e., n = 138 out of 292), and the correlation between P and LREE is observed for 53% of the grains. This suggests that the coupled «xenotime» substitution controlled the incorporation of LREE in most grains. The P-enrichment observed for part of the grains of the D3 sample could be related to minute apatite inclusions incorporated into spot analyses (Figure 9c). By contrast, in sample D1A, apatite inclusions are sufficiently large to be identified on CL images and could be more easily avoided during analysis (Figure 9c). The diagram of Hoskin [52], established for pristine and hydrothermally altered zircon from TTG intrusive suites, indicates that most of the studied zircons are altered, with grains from the D1A sample being the most altered (Figure 9d).

Figure 9. Chemical composition of the studied zircon grains displayed on the binary diagrams (a) Zr vs. Fe + Ca + Sr + La, (b) Th vs. Lu, (c) sum LREE vs. P, and (d) La_N vs. (Sm/La)_N. On the latter diagram, from Hoskin [52], values are normalized to chondrite.

Figure 9. Chemical composition of the studied zircon grains displayed on the binary diagrams (a) Zr vs. Fe + Ca + Sr + La, (b) Th vs. Lu, (c) sum LREE vs. P, and (d) La_N vs. (Sm/La)_N. On the latter diagram, from Hoskin [52], values are normalized to chondrite.

4.6. Calculating the fO₂ and Temperature Parameters

The fO₂ and temperature parameters are calculated for the D1A and D3 samples (Figure 10a) using the method of Smythe and Brenan [53]. The method mostly relies on the Ce content in melt and zircon, as well as on the activity of silica (a_SiO2) and titanium (a_TiO2), and on the composition of the melt. The a_SiO2 parameter is set to 1 because the studied rocks contain interstitial quartz. The a_TiO2 parameter, presented above, is set to 0.6 [46]. The composition of the melt is unknown, as melt inclusions are lacking in the studied zircon grains and is approximated using whole rock chemical composition. The water content of the melt is unknown. The water content of the melt is 4–5 wt% H₂O for other intrusive systems of the Chibougamau area (i.e., alkaline intrusion, Opémisca granodiorite-monzodiorite intrusion, Chevrillon granodiorite-monzonite pluton, and Muscocho granodiorite pluton) for which the water content of the melt could be estimated using amphibole chemistry [13,14]. The water content parameter is thus set to 4 wt% H₂O for samples D1A and D3. Tests are also performed using higher H₂O values (not shown) and indicate that increasing the water content parameter shifts fO₂ estimates toward higher values.

For the D1A sample, the La content increases as the estimated fO₂ value increases (Figure 10a). As pristine grains tend to be La-poor, zircon grains with the lowest La values likely return the most reliable estimates of the fO₂ parameter. In sample D1A, the calculated oxidation state is FMQ −3.5 ± 1.3 (median ± standard deviation − std) for the n = 140 grains with La ˂ 10 ppm, and FMQ −2.1 ± 2.7 for the n = 4 grains with La ˂ 1 ppm. Zircon, in sample D3, points to reducing conditions estimated at FMQ −3.94 ± 0.61 using the 149 grains that contain ˂1 ppm La.

The fO₂ and temperature parameters are also calculated using the method of Loucks et al. [54] (Figure 10c). This empirical method requires insights into the U, Ce, and Ti content of zircon, and knowledge of the crystallization age. Dating is used to estimate the eU parameter, which is the age-corrected initial U content of zircon. The age parameter is set to 2690 Ma for samples D1A and D3. In sample D1A, the calculated oxidation state is FMQ −0.58 ± 1.3 (median ± std) for the n = 140 grains with La ˂ 10 ppm, and FMQ +0.80 ± 2.0 for the n = 4 grains with La ˂ 1 ppm. In sample D3, an oxidation state of FMQ +2.42 ± 0.64 is obtained for the 149 grains that contain ˂1 ppm La (Figure 10c).

Using the Ti-in-zircon thermometer [55] implemented by the method of Smythe and Brenan [53], and using a value of 0.6 for the a_TiO2 parameter, a temperature estimate of 672 ± 53 °C (median ± std) is obtained for the zircon grains of sample D3 that contain ˂1 ppm La (n = 149 grains) (Figure 10b). The zircon grains of sample D1A return temperatures comprised between 651 °C and 1488 °C, with median and std values of 849 ± 131 °C (Figure 10b). The La-poor grains in sample D1A, return temperatures of 800 ± 106 °C (for the n = 140 grains with La ˂ 10 ppm) and 768 ± 46 °C (for the n = 4 grains with La ˂ 1 ppm).

4.7. Comparison with Another Carbonatite-Bearing Intrusive Complex

The chemistry of the studied zircon is further discussed using a comparison to zircon from a carbonate-nepheline-bearing syenite from the Crevier carbonatite-alkaline complex studied by Groulier et al. [56]. This intrusive complex is located in the Proterozoic Grenville orogeny, in Québec, and counts as one of the few alkaline-carbonatite complexes of the North American continent for which zircon data are available. Zircon grains, in the Crevier complex, are classified as: (1) euhedral, which are grains enriched in Ca and Fe and crystallized from a carbonatite-saturated (pre-immiscibility) alkaline melt; and as (2) anhedral, which correspond to grains least enriched in Ca and Fe and that may have crystallized in a post-immiscibility syenite melt [56].

Displaying chemical data as a function of the non-stochiometric elements Ca and Fe (Figure 11) shows that most zircon grains in sample D3, are similar to the anhedral (Ca- and Fe-poorest) zircon from the Crevier complex, i.e., their Th/U ratio, Sr content, and Ti-in-zircon temperature are comparable (Figure 11a,c,d). The Douay zircon grains (D3 sample) are, however, enriched in LREE compared to the Crevier zircon grains (Figure 11b,e). For the Douay intrusion, the fO₂ parameter is uncorrelated to the Ca + Fe content (Figure 11f) of the grains. The fO₂ parameter calculated using the method of Loucks et al. [54] is higher in sample D3 compared to the sample from the Crevier complex (Figure 11f). For the D1A sample, zircon grains span the range of Ca + Fe values displayed by the zircons of the Crevier complex and have little in common with the anhedral and euhedral Crevier grains (Figure 11). The positive correlation between Ca + Fe and Th/U, LREE, Sr, temperature, and (Sm/Yb)_N observed at Crevier is not observed for sample D1A (Figure 11).

5. Discussion

5.1. Rock Type and Intrusive Phases

The studied rocks (except sample D3) have the chemistry and petrography of syenite (Figure 3 and Figure 4); however, petrographic observations indicate that feldspar has been modified by alteration and the pre-alteration rock may have contained plagioclase, now replaced by alkali feldspars. In addition, mass balance calculations indicate that the alkali content has been increased by hydrothermal alteration, and this modified the location of samples on the TAS diagram (Figure 4a). Another consequence of alteration-induced Na-K-gains is that it would have decreased the ASI value, i.e., precursors may have been metaluminous. Using the method of Trépanier et al. [39], it can be concluded that the precursors have the chemistry of quartz-monzonite. The studied rocks are also less enriched in incompatible elements compared to the syenite intrusions of the Abitibi greenstone belt (Figure 4c,d). We thus conclude that the studied rocks (samples D1A, D2, D4, and D5) are quartz-monzonite and not syenite. As other parts of the Douay intrusion consist of syenite and minor carbonatite phases [11], the identification of a quartz-monzonite phase in the mineralized Porphyry zone (Figure 1) indicates that Douay is a polyphase intrusion or intrusive suite. Magmas emplaced in the area of the Douay project are thus moderately alkaline or high-K calc-alkaline (altered quartz-monzonite phase) and alkaline (syenite and carbonatite phases). Douay is similar to other intrusive suites of the Abitibi greenstone belt (e.g., Kirkland Lake area), which formed during the syntectonic period and consist of quartz-monzonite, syenite, and other intrusive phases [3].

Hydrothermal alteration, for the quartz-monzonite phase, is characterized by Na-metasomatism and minor K-metasomatism, as well as moderate to negligeable silica loss. During the alteration process, magmatic plagioclase and ±quartz were likely destroyed as hydrothermal alkali feldspar, which is enriched in silica compared to plagioclase, crystallized. Another important hydrothermal alteration is carbonatization, which produced disseminated carbonate and carbonate veinlets, and that induced negligeable to moderate Ca-gain. Sulfidation formed a minor amount of pyrite and rare chalcopyrite on the studied outcrop, because this outcrop is located at the periphery of the mineralization of the Porphyry zone. Within the nearby deposit, pyrite ranges from 0.5% to >10%, and chalcopyrite is a minor phase (El Rassi [11] and F. Speidel pers. comm.). Alkali metasomatism is consistent with the circulation of magmatic fluids in the studied rocks. Given the context, it is unclear whether these fluids were exsolved from the studied quartz-monzonite intrusive phase (i.e., autometasomatism), or whether these fluids were produced by syenite and carbonatite phases. Other alterations, such as silicification (quartz veinlets observed in the field) and carbonatization, are consistent with both magmatic-hydrothermal and orogenic gold mineralizing systems, and, thus, have an unclear origin.

In quartz-monzonite rocks, hydrothermal alteration also had a significant impact on mafic magmatic minerals, such as amphibole and/or clinopyroxene. These minerals have been altered to clusters of biotite (indicating local K-gain), carbonate (C-gain), chlorite, epidote, feldspar, and magnetite (hydration and local reorganization of Fe, Mg, Ca, Al, and Si). Given the extent of alteration, it is unclear whether magnetite is a secondary mineral or whether it corresponds to a recrystallized magmatic mineral. The pre-alteration mineralogy of the quartz-monzonite phase thus comprises two feldspars, quartz, and ±amphibole (or clinopyroxene?), apatite, zircon, and possibly magnetite.

The aplite dyke (sample D3) is less carbonatized than the quartz-monzonite intrusive phase. Mass balance calculation could not be performed for this sample and petrographic observation suggests that Na-metasomatism modified the outer margins of feldspars (Figure 3f). Sample D3 is a granite, and whether its precursor was an alkaline or a sub-alkaline rock that gained a significant amount of Na during alteration is unclear. The trace element profile of sample D3 is parallel to those of the quartz-monzonite phases. The aplite dyke may thus correspond to a volatile-rich residual melt derived from the quartz-monzonite magma. Sample D3 contains a lesser amount of trace elements compared to the other samples (Figure 4e,d), possibly because the residual melt separated after prolonged crystallization of trace elements-enriched phases such as zircon, apatite, and amphibole. The pre-alteration mineral assemblage of sample D3 likely consisted of two feldspars, quartz, and ±amphibole (or biotite?), apatite, zircon, and possibly magnetite.

5.2. Timing of Zircon Crystallization

In the studied rocks, most minerals have been modified by hydrothermal fluids and zircon is the main magmatic mineral that can be used to further document the nature and physical-chemical properties of the studied magma. Prior to examining zircon chemistry, the timing of zircon crystallization is discussed. Whole rock chemistry points to Zr-poor magmas, which is consistent with the possible metaluminous nature of the melt (see previous section), as experimental studies indicate that metaluminous (and peraluminous) melts generally contain ~100 ppm Zr or less [57,58], whereas peralkaline melts can dissolve ˃1 wt% Zr [59,60]. Zircon saturation likely occurred late as a consequence of the limited availability of Zr in the magmatic system.

In sample D1A, zircon grains are euhedral and anhedral (incomplete fragments), with this latter crystal form indicating that zircon crystallized in melt pockets with a relatively limited volume [47]. This is consistent with the crystallization of zircon during the last stage of the evolution of the magmatic system. In addition, complete grains are strongly elongated (aspect ratio = 0.3 to 0.4), suggesting a fast growth rate [47]—this is expected for the feldspar-phyric quartz-monzonite phase (Figure 2a), which likely emplaced at a shallow depth as a feldspar-bearing magma and cooled rapidly.

In zircon of the D1A sample, core regions can be distinguished on CL images, that may correspond to inherited grains, and relics of oscillatory zoning are observed, confirming that these zircons have a magmatic origin (Figure 5). Zircon in sample D1A, has also been impacted by post-magmatic processes that coarsened and blurred the oscillatory zoning in some grains [61]. In addition, most grains are dark on CL images, porous, and have lost their oscillatory zoning, indicating that zircon is partly amorphous. This is a likely consequence of radiation damage to the crystal structure resulting from U-decay [62]. Zircon grains in sample D1A, behave as completely open systems and are unsuitable for dating (Figure 7).

The aplite dyke (sample D3) was likely emplaced at temperature conditions of ˃570 °C. Zircon is recovered as incomplete fragments with partially developed crystalline faces, indicating that it crystallized from late pools of highly fractionated residual melt [47]. Zircon grains from sample D3 lack core regions that can be easily distinguished from the rest of the grain, and they lack multiple age populations (Figure 7e), i.e., they lack inherited cores. This is a common feature of aplite dykes, which formed from late residual melts that are unlikely to contain inherited zircon [63]. Oscillatory zoning is irregular in some grains and can be cut by irregular domains with homogeneous composition indicating a partial recrystallization of the grains. Recrystallization expelled incompatible (large ionic radius) trace elements in parts of the grains and concentrated impurities in trace elements-rich convolute zones, producing bright and dark CL domains (Figure 6). Such recrystallization features can be promoted by the circulation of magmatic fluids [47], as is likely the case here.

In addition, the bright CL corroded rims observed for both samples are typical of zircon that has interacted with an aqueous solution [62]. These rims correspond to transgressive recrystallization patches from which trace elements have been expelled [61]. In both samples, zircon likely interacted with magmatic fluids either exsolved from the studied intrusion (autometasomatism) or from other intrusive phases, and subsequent interaction with metamorphic fluids cannot be excluded. The main difference between samples D1A and D3 is that zircon in sample D1A was likely metamict prior to alteration, while those of sample D3 were pristine. Metamictization, in sample D1A, would have promoted the destruction of the crystal and exchanges of elements with the fluids, inducing the observed enrichment in trace elements.

As indicated previously, the aplite dyke (sample D3) may correspond to a residual melt extracted from the quartz-monzonite intrusive phase (other samples), implying that zircon grains in samples D1A and D3 have similar crystallization ages. That a same alteration event impacted metamict (sample D1A) and pristine (sample D3) zircon grains with similar crystallization ages implies that metamictization occurred most rapidly in sample D1A.

It is possible to estimate the time required for complete metamictization of zircon. Here, the calculation is performed using Equation (5) defined by Meldrum et al. [64]. This equation is derived from textural observations performed on synthetic and natural zircon with known age and U content [64]. The method of Meldrum et al. [64] calculates the minimum amount of eU (age-corrected initial U content of zircon) required for a zircon grain to become entirely metamict within a given time interval and at a defined temperature. Several time intervals—15 My, 50 My, and 100 My—and temperatures—50 °C and 100 °C—are tested (Figure 12a). The temperature at which the ~2690 Ma intrusions of the Douay intrusive suite were maintained after complete cooling of the magmatic system are unconstrained because the emplacement depth of the intrusion, as well as the paleo-geothermal gradient, are unknown. Temperatures of 50 °C and 100 °C are considered realistic for an intrusive complex likely emplaced a few kilometers below the surface. Further constrains on the temperature parameter are unnecessary, because temperature variation has a limited impact on the calculated eU parameter (Figure 12). The method of Meldrum et al. [64] shows that zircon from the D3 sample would have required ˃˃ 100 My to become metamict. Most zircon grains from sample D1A, however, would be entirely metamict under 100 My, and would likely be partially metamict within a shorter time interval of 15 My (Figure 12a).

The amount of radiation damage in zircon can also be estimated using the radiation dose in displacements-per-atom (D_dpa). The D_dpa parameter is calculated using the equations reported by Gao et al. [65]. This parameter combines U and Th concentrations, as well as the time span of decay, i.e., 15 My in the example shown here (Figure 12b). The results show that, for a 15 My time interval, the D_dpa parameter is >0.01 for most zircon of sample D1A, and <0.005 for sample D3 (Figure 12b). To further interpret the D_dpa parameter, it is compared to values calculated for the Suzhou granite, whose zircon grains show a ‘high D_dpa effect’ at D_dpa > 0.03 [65]. For sample D1A, most grains return D_dpa values that are also close to, or above 0.03 (Figure 12b). Metamictization, in sample D1A, may thus have been significant a few Myrs after complete cooling of the quartz monzonite phase and dissipation of the magmatic fluids exsolved from this intrusion, and prior to subsequent interaction with hydrothermal fluids. Sample D3 indicates that the aplite dyke, as well as the likely co-magmatic quartz-monzonite phase, emplaced at ~2690 Ma, i.e., 10–20 My prior to the emplacement of the syenite and carbonatite phases dated at 2676 +5/−6 Ma [23]. We propose that the hydrothermal process that has altered the studied zircon results from interaction with magmatic fluids, exsolved from the alkaline phases (syenite and carbonatite). Additional alteration induced by the subsequent circulation of metamorphic fluids during the D₂–D₃ deformation events (or after) cannot be excluded.

The quartz-monzonite and granite phases (main intrusion and aplite dyke) and the syenite-carbonatite intrusive phases at Douay, correspond to magmatic events with distinct ages and chemistry, and these two intrusive stages are unlikely to be co-magmatic. However, the occurrence of successive intrusive events over a restricted area points to a structural configuration that favored magma injection in the upper crust during the syntectonic stage. This is important from a metallogenic point of view as structures that can facilitate the transport of magma through the crust can also channel mineralizing fluids. This study also has implications for our understanding of magmatism in the Abitibi greenstone belt. Indeed, the quartz-monzonite of the Douay intrusive suite is part of a broader magmatic event in the northern part of the belt. This quartz-monzonite phase is coeval with several alkaline and subalkaline intrusions, such as the 2691 +5/–3 Ma Lac Shortt intrusion (see above), the 2692.1 ± 2 Ma O’Brien syenite-monzonite-granodiorite stock [3,8], the 2688 ± 8 Ma Saussure syenite [14], and the Chevrillon quartz-monzonite pluton [66]. The ~2676 Ma syenite-carbonatite intrusive phases of the Douay intrusion emplaced late and possibly toward the end of (or after) the syntectonic period, which is generally considered to extend from 2700 to 2690 Ma in the northern part of the Abitibi greenstone belt [1,20,21]. These syenite-carbonatite intrusive phases may be coeval with the Dolodau intrusion, which is approximately dated at 2630–2680 Ma [67,68,69]. The syenite-carbonatite intrusive phases of the Douay suite is also coeval with alkaline and subalkaline magmatism in the southern part of the belt, where the syntectonic period extends from 2690 to 2665 Ma [1,17].

5.3. Physical-Chemical Properties

Zircon saturation temperatures of 667 °C (sample D1A) and 570 °C (D3) can be obtained using a method [44] that relies on whole rock chemistry. These temperatures are unrealistic because they are close to below the wet granite solidus [70]. The method of Watson and Harrison [44] is thus imprecise in the study area, possibly because bulk rock composition has been modified by hydrothermal alteration.

Zircon crystallization temperature can also be calculated using a Ti-in-zircon thermometer [53], given that the a_TiO2 parameter can be estimated. This parameter is likely overestimated when calculated from whole rock analysis [44], as whole rock is unrepresentative of the chemical composition of the residual melt pockets out of which the studied zircon crystallized, and we elected to use a generic value of 0.6 for the a_TiO2 parameter [46]. The Ti-in-zircon thermometer also requires data on the Ti content of zircon, analyzed here using the LA-ICP-MS instrument. Finally, the thermometer method requires that Ti be incorporated in zircon following Henry’s Law behavior and that its distribution is not modified by alteration. This latter assumption may be questioned given the extent of alteration observed in the studied rocks. To limit the impact of alteration on the results, the temperature is estimated for the least-altered (˂1 ppm La) zircon grains of the D3 and D1A samples, which return temperatures of 672 ± 53 °C and 768 ± 46 °C, respectively. The Ti-in-zircon thermometer thus indicates that zircon crystallized at ~700–800 °C in the quartz-monzonite phase, prior to crystalizing at ~600–700 °C from the most evolved residual melt (aplite dyke) extracted from the quartz-monzonite phase.

To estimate the fO₂ parameter, we used the method of Smythe and Brenan [53] that relies on several parameters, with the most important being: (1) the Ce content of zircon; and (2) the estimate of the water content of the melt (i.e., 4 wt% H₂O). Other parameters, such as a_SiO2, a_TiO2, and melt composition (water excepted) have a lesser impact on final results [13,14]. The fO₂ parameter is also estimated using the method of Loucks et al. [54] that is based on the U, Ce, and Ti content of zircon. Both methods return fO₂ values that are highest for La-rich than for La-poor grains (Figure 10). This indicates that alteration may affect the concentration of Ti, U, and/or Ce. Temperature estimates obtained using the Ti-in-zircon thermometer are also correlated to La content (Figure 10b), suggesting that Ti is indeed affected by alteration. Uranium is strongly compatible in zircon, as it is incorporated by a simple substitution mechanism (U(IV) = Zr(IV); [49]), and may thus be difficult to remove from zircon during an alteration process; however, U can be mobilized by F-Cl-bearing fluids [71]. Such fluids may have circulated in the Douay mineralized system, where fluoritization is reported and where base metals-bearing sulfides (chalcopyrite, molybdenite, sphalerite, galena) are observed [11] that are indicative of Cl-bearing fluids. Modification of the initial U content of the studied zircon cannot be excluded. Cerium, on which the fO₂ calculation relies extensively, has a ionic radius close to this of HREE and, like HREE, is abundant in zircon [49]. Results presented here show that the LREE content varies more than the Ce and HREE contents (Figure 8), suggesting that the alteration has a limited impact on Ce concentration. In addition, Ce(IV) is strongly compatible in zircon because it likely substitutes for Zr(IV), and Ce(IV) is unlikely to be easily removed from zircon and to be mobile in the presence of a fluid.

For sample D1A, the least-altered (˂1 ppm La) zircon grains return values of FMQ −2.1 ± 2.7 (method of Smythe and Brenan [53]) and of FMQ +0.8 ± 2 (method of Loucks et al. [54]). The fO₂ value obtained using the method of Loucks et al. [54] compares to the FMQ 0 to +1 values obtained for several syntectonic intrusions of the Chibougamau area (i.e., the Saussure, Opémisca, Chevrillon, and Muscocho plutons) [13,14] and also compares to the oxidizing conditions generally reported for alkaline systems of the Abitibi greenstone belt, e.g., the Murdock Creek intrusion [72]. Moderately oxidizing conditions of FMQ +0.8 ± 2 are also consistent with the occurrence of magnetite in the Douay samples (see Section 4.3). We thus conclude that the method of Loucks et al. [54] returns results that are more realistic than the low fO₂ values returned by the method of Smythe and Brenan [53]. A possible explanation for the poor performance of this latter method is an imprecise estimate of the water content of the melt. Additional calculations (not shown) indicate that increasing the water content of the melt to 7–10 wt% H₂O is sufficient for the method of Smythe and Brenan [53] to return realistic estimates of the oxidation state of ~FMQ 0 to +1. We thus conclude that the magma was enriched in water, and this is consistent with the shape (i.e., {110}) and {101} as the dominant form) of zircon in the D1A sample.

For the D3 sample, the method of Loucks et al. [54] points to oxidizing conditions of FMQ +2.42 ± 0.64, suggesting an increase in the oxidation state in the highly evolved residual melt injected to form aplite dykes. To obtain similar results using the method of Smythe and Brenan [53], water content of up to 30–40 wt% H₂O is required. Even if a high water content is expected for a volatile-rich residual melt injected as an aplite dyke, values of 30–40 wt% H₂O are considered to be unrealistically high. Modification of the melt composition parameter (e.g., decrease of the alkali content) also shifts the results provided by the method of Smythe and Brenan [53] toward higher fO₂ values. This latter hypothesis is favored as zircon in sample D3 crystallized within isolated pools of residual melt, whose chemistry is likely poorly represented by whole-rock chemistry.

Results show that, even in the strongly altered zircon population of sample D1A, a few grains have retained their initial Ce content and return a fO₂ value that is close to the value of FMQ 0 to +1 obtained for several syntectonic intrusions of the Chibougamau area [13,14]. We thus conclude that the Ce content of pristine to moderately altered zircon can be used to document the fO₂ parameter. In altered rocks, constraining the parameters required by the method of Smythe and Brenan [53] is difficult, while the method of Loucks et al. [54] performs best for such rocks.

5.4. Chemical Composition of Zircon

Zircon from sample D1A is significantly enriched in all trace elements (Eu excepted) compared to zircon from the D3 sample (Table 3). This trace elements enrichment in sample D1A, is likely a consequence of hydrothermal alteration that impacted U-rich and radiation-damaged (metamict) zircon grains, which is consistent with alteration through a diffusion-reaction process [62]. Zircon of sample D1A is also characterized by a complete re-equilibration of the isotopic system (Figure 7a), inclusions of thorite-coffinite, and micrometer-sized pores, which are characteristics generally attributed to the dissolution-reprecipitation process [62]. The diffusion-reaction process is considered more likely here, because of the extreme enrichment in non-formula elements observed in sample D1A (Table 3) and because those zircon grains were likely metamict prior to alteration. Zircon in sample D3 was likely altered through a dissolution-reprecipitation process, which impacts pristine grains [62]. In this section, we will examine additional impacts of alteration on the chemistry of zircon grains.

A first consequence of the incorporation of trace elements is to modify the Zr/Si ratio of zircon, which is higher in the D3 than in the D1A sample (Figure 9a). In addition, zircon in sample D1A is enriched in common Pb, which is a known consequence of alteration [73]. Zircon is also Eu-depleted in sample D1A (Figure 8a). This is unlikely to be a consequence of crystallization from an Eu-depleted melt and/or to crystallization from a reduced (Eu(II)-dominated) magma, as the studied rocks lack negative Eu anomalies (Figure 4) and as zircon from the most evolved melt (sample D3) also lack negative Eu anomalies (Figure 8b). A more likely explanation for sample D1A is that Eu and possibly other trace elements were preferentially removed from the radiation-damaged grains before being re-introduced by hydrothermal fluids unable to transport Eu(III) efficiently. These fluids may have transported Eu(II) but, as only Eu(III) can substitute into the zircon lattice [74], zircon grains in sample D1A maintained the Eu-depletion that they acquired during metamictization.

The Eu(III)-depleted fluids discussed here likely correspond to high temperature fluids. Indeed, high temperature fluids mobilize Eu(II) and not Eu(III) according to experiments performed using aqueous and chloride-rich solutions [75]. Additional experiments are, however, required to evaluate the mechanism of Eu transportation in the oxidized, and likely CO₂-, F-, and LREE-bearing magmatic fluids exsolved from the ~2676 Ma syenite and carbonatite phases of the Douay complex, and which correspond to the fluids that likely re-introduced trace elements in zircon of the D1A sample.

The positive correlation between Th and Lu, which is a consequence of a simple substitution mechanism between Zr and trace elements in pristine zircon, is lost in sample D1A. Zircon in this sample displays enrichment in both Th and Lu, with Lu ˃ Th (Figure 9b). This is likely a consequence of the low mobility of Th in a hydrothermal fluid [76], which is more likely to migrate toward thorite micro-inclusions within zircon, while the syenite-carbonatite complex emplaced at ~2676 Ma may have produced Th-poor and REE-bearing fluids, rendering Lu available to the radiation-damaged zircon of sample D1A.

Compared to zircon from the Crevier intrusive complex, zircon from sample D3 has similar chemistry to that of zircon interpreted to have crystalized from a carbonate-poor syenite melt of the Crevier complex [56]. Only the LREE enrichment is more pronounced at Douay, inducing higher values for the Sm/Yb ratio (Figure 11b,e). The REE content of zircon depends on the REE content of the melt and on co-crystalizing minerals [77], and the melt composition and crystallization sequence may have been different at Douay and Crevier. Moreover, the LREE-enrichment displayed by zircon of sample D3 could be a consequence of alteration induced by REE-enriched magmatic fluids at Douay. Regarding the D1A sample, the chemistry of altered zircon is different from that of the zircon grains of the Crevier complex that display similar Fe + Ca contents. Alteration in the D1A sample induced a significant enrichment in Ca, Fe, Sr, and LREE that can be compared to that observed at Crevier; however, the Th/U ratio, the Sr vs. Fe + Ca positive correlation, and the Ti-in-zircon temperature are not comparable (Figure 11) and can be used to distinguish strongly altered zircon (sample D1A) from zircon that crystallized from a carbonate-rich alkaline magma (Crevier complex).

5.5. Discrimination Diagrams

The impact of the diffusion-reaction process on the chemistry of zircon in the D1A sample is further investigated using a comparison with the D3 sample, and with samples from other syntectonic and synvolcanic plutons of the Chibougamau area and compiled from recent studies [13,14]. This comparison is presented in detail in the electronic Supplementary Material Figure S1, and results indicate that U, Th, Pb, Y, Sr, and Eu are diagnostic of magma type and alteration processes. The Th/U ratio is also useful to distinguish crystallization environments (Figure 11a). The Th/U ratio is well characterized and generally takes values of 1, or 0.2–0.8, in magmatic zircons [78], while it takes values of ˂0.1 in metamorphic zircon due to the expulsion of Th from the structure of zircon [77,79]. In sample D1A, the Th/U ratio is within the range expected for magmatic zircon but the U, Ca, and Fe contents are much higher than in most magmatic zircon, enabling a clear distinction between extensively altered and least-altered zircon (Figure 11a).

Following these observations, several discrimination diagrams are tested to facilitate the interpretation of zircon chemistry. The Th/U vs. La binary diagram (Figure 13a), which is equivalent to a Th/U vs. Sr diagram (not shown), can discriminate between zircon that crystallized in carbonate-nepheline-bearing syenitic magma (Crevier carbonatite-alkaline complex), zircon grains that were modified by a diffusion-reaction process following interactions with magmatic fluids exsolved from a carbonatite-alkaline complex (Douay intrusion), zircon modified by metamorphism (with Th/U ˂ 0.1, see above), and other zircons, i.e., pristine and altered zircon crystalized in alkaline and sub-alkaline magmas. This latter category can be further discriminated using a Th/U vs. U binary diagram, as sub-alkaline magma produced zircons that contain ~100–1000 ppm U and with a Th/U ratio comprised between 0.2 and 1 (field A in Figure 13b), while alkaline magma (Crevier complex and Saussure intrusion) produce zircon with a chemistry that falls in field B in Figure 13b. Zircon in sample D3, falls in fields A and B (Figure 13b), possibly because it crystallized from a sub-alkaline magma and was subsequently partially altered by magmatic fluids exsolved from an alkaline magma.

6. Conclusions

The Douay intrusive complex, Abitibi greenstone belt, hosts part of the gold of the Douay exploration project. The petrogenesis of the Douay intrusion, as well as the metallogenic processes that led to the concentration of Au in the upper crust in and near this intrusion, remains debated. The aim of this study is to demonstrate that zircon from altered intrusive rocks can be used to infer the age of magmatic and hydrothermal events. The Douay intrusion was known as a ~2676 Ma syenite-carbonatite intrusive complex and this study demonstrates that Douay is an intrusive suite that also comprises a ~2690 Ma quartz-monzonite intrusive phase. Both intrusive phases are unlikely to be co-magmatic given the age and chemical differences between these magmatic events. Investigation of zircon chemistry shows that the following sequence of events occurred: (1) crystallization of zircon at ~2690 Ma in an evolved residual melt enriched in trace elements, U-included (sample D1A), and in a highly evolved residual melt depleted in trace elements (sample D3); (2) extensive radiation damage for the U-rich grains (sample D1A) over a period of ~10–15 My; and (3) alteration of zircon grains by interaction with magmatic fluids likely exsolved from the ~2676 Ma syenite-carbonatite Douay complex. Alteration is controlled by the diffusion-reaction process (sample D1A) and the dissolution-reprecipitation process (sample D3). Characterization of zircon chemistry led to the identification of the Th/U vs. U discrimination diagram that can be used to distinguish extensively altered zircon grains from pristine to least-altered zircon formed in distinct magmatic environments (Figure 13b).