Multi-scale isotopic heterogeneity reveals a complex magmatic evolution: An example from the wallundry suite granitoids of the lachlan fold belt, Australia

Open-system magmatic processes are expected to impart various sorts of isotopic heterogeneity upon the igneous rocks they produce. The range of processes under the “open-system” umbrella (e.g., simple two-component mixing, magma mingling, assimilation with fractional crystallization) cannot usually be uniquely identified using data from a single isotope system. The use of bulk-rock, mineral separate and in situ techniques and multiple isotope systems allows the characterization of isotopic variability at different sampling scales, illuminating details of the petrogenesis of a magmatic system. This approach has been applied to granitoids of the Wallundry Suite in the Lachlan Fold Belt, Australia. The Wallundry Suite exhibits variations in mineral assemblage, mineral composition and trends in bulk-rock major and trace element compositions consistent with the involvement of liquid-crystal sorting processes such as fractional crystallization. In situ paired O-Hf isotope data from zircon in six samples show an array indicating the isotopic evolution of the melt phase. Similarly, bulk-rock Sr-Nd-Hf isotope arrays support open-system magma evolution. These data combined with the petrographic observations and major and trace element geochemical variations suggest some form of assimilation-fractional crystallization process in the petrogenesis of the Wallundry Suite. Added complexity is revealed by two observations: 1) the isotopic variations are only weakly coupled to the lithology and major element compositions of the samples; and 2) there are distinguishable differences between the Hf isotope compositions of bulk-rock samples and those of the magmatic zircons they host. To varying degrees the rocks consistently show negative ΔεHfbulk-zrc values (i.e., the bulk-rock compositions have less radiogenic Hf isotope values than their coexisting zircons). The preservation of distinctly low Nd and Hf isotope ratios in an Fe-Ti oxide mineral separate suggests that the bulk-rock vs. zircon discrepancy is caused by the presence of unmelted components derived from a contaminant of continental origin (i.e., a rock with low Sm/Nd and Lu/Hf and thus unradiogenic Nd and Hf). Evidently, a complex interplay of assimilation, crystallization and melt segregation is required to account for the data. This investigation demonstrates that such complexity can, nevertheless, be disentangled through comparison of complementary isotope data at multiple sampling scales.


Introduction
A critical aspect of granitoid research has been identifying sources that contribute to the production of melt in the source region, and processes involved throughout magma evolution until its final crystallisation in the upper crust. Constraining these parameters relies on elucidating the origin of observed mineralogical and geochemical variations, often at different scales (see, for example, Chappell and White, 1974;Reid et al., 1983;Watson and Harrison, 1983;Williams et al., 1983;Chappell, 1984;Gray, 1984;Wyborn and Chappell, 1986;Chappell et al., 1987;Clemens and Vielzeuf, 1987;Chappell and Stephens, 1988;White and Chappell, 1988;Vernon, 1990;Patiño Douce and Johnston, 1991;Vielzeuf and Montel, 1994;Patiño Douce, 1995;Collins, 1996;Sawyer, 1996;Singh and Johannes, 1996;Beard et al., 2005;Brown, 2007;Kemp et al., 2007;Sawyer et al., 2011;Clemens and Stevens, 2012;McLeod et al., 2012;Paterson et al., 2016). The many models for granitoid genesis range from closedsystem, single-source (be they crust or mantle) models to open-system, multi-source models. These include crystallisation-differentiation (fractional crystallisation), mixing of pure melts or of (crystalbearing) magmas, solid entrainment models (i.e., restite unmixing or peritectic assemblage entrainment, PAE), and assimilation through bulk-digestion or coupled to fractional crystallisation (AFC). In different contexts and by different researchers, assorted petrogenetic models are viewed as competing hypotheses, envisaged as complementary processes, or invoked in combination with each other.
Fractional crystallisation has long been regarded as an important process in the evolution of both intrusive and extrusive rocks, and is in many ways the classic form of magmatic differentiation (for a historical overview, see Marsh, 2006;Wilson 1993). Bowen (1947) emphasised the ability of fractional crystallisation to generate a range of compositions from basaltic to granitic along a liquid line of descent. It has also long been recognised that separating a residual granitic liquid from the considerable quantity of crystalline material required to produce it in this fashion is problematic, and many potential mechanisms to achieve this have been proposed (Bowen, 1919;Bowen, 1947;McKenzie, 1985;Marsh, 2006). Marsh (2006) suggested punctuated differentiation, in which extensive settling of phenocrysts (formerly carried by ascending magma) occurs upon emplacement, to segregate relatively large volumes of fractionated liquid. The modelling of Dufec and Buchmann (2010) suggests melt extraction is optimised in particular crystallinity windows (~50-70%) such that a fractionated liquid is likely to escape from a crystalline mush of a distinctly different composition (explaining compositional gaps observed in many volcanic settings). Thus, those authors proposed that fractional crystallisation could produce a step-wise compositional evolution of both volcanic and plutonic rocks, which is consistent with the broadly layered structure of the continental crust.
For virtually as long as fractional crystallisation has been considered an important form of magmatic differentiation, geologists have been cognisant of the possibility that this process may be accompanied by assimilation of the surrounding rocks (McBirney, 1979;Wilson 1993; and references therein). Numerous quantitative models for AFC have been derived (e.g., Allègre and Minister, 1978;Taylor, 1980;DePaolo, 1981;Nielsen, 1989;Hagen and Neumann, 1990;Aitcheson and Forrest, 1994;Bohrson and Spera, 2001;Spera and Bohrson, 2001;Bohrson et al., 2014;Bohrson et al., 2020;Heinonen et al., 2020), showing that isotope systems and some trace elements are likely to be strongly affected by progressive assimilation. A review of the theory and implementation of such models was undertaken by Heinonen et al. (2021), including a comparison of possible outputs for different modelling approaches. Typically (and critically) such models have emphasised departures for trends in isotope-isotope and trace element-isotope space expected to arise through AFC from those expected for simple liquid-crystal sorting or binary mixing. This will especially be the case when magmatic system modelling is thermodynamically constrained and the trace element and isotopic evolutions are tied to phase equilibria modelling, as is the case with the Magma Chamber Simulator (MCS) modelling program Heinonen et al., 2020).
Regardless of whether any particular model was initially developed to examine intrusive or extrusive samples, fundamentally they describe the evolution of a crystallising magma chamber and may, therefore, be applied to the formation of granitoids. Taylor (1980), for example, modelled Sr and O isotope variations during combined assimilation and fractional crystallisation to the role of such a process in the formation of both the Adamello Massif and the Roccamonfina volcanics. Furthermore, Iles and Heinonen (2022) modelled bulk rock major and trace element trends using the MCS to investigate the formation of the Jindabyne Suite granitoids from Australia.
A large body of literature focussing on Australian granites debates the roles of restite unmixing (a closed-system process) versus a whole family of open-system models (i.e., various forms of magma mixing). Fractional crystallisation is viewed as an accompanying process in many magma mixing models concerning the granites of the Lachlan Fold Belt (LFB), Australia (Collins, 1996;Keay et al., 1997;Kemp et al., 2007;Kemp et al., 2009). In contrast, proponents of restite unmixing (White and Chappell, 1977) have seen fractional crystallisation as a complementary process. Chappell et al. (1987) suggested that fractional crystallisation might replace restite unmixing as the differentiation process in a felsic liquid from which essentially all restite has been removed. Furthermore, it was recognised that some magmas (which might include basaltic melts) would leave their source region free of residual solids and ascend through the crust as essentially pure melts (Chappell et al., 1987). The mechanism for producing petrographic and chemical variations in granite (sensu lato) suites that result from such cases would be fractional crystallisation.
In this context, few LFB granitoids (those classified as hightemperature I-types according to the scheme of Chappell et al., 1998) have escaped the contention associated with magma mixing versus restite unmixing, and are considered to have evolved via fractional crystallisation. One such example from the LFB is the Wallundry Suite. Using these granitoids as an example, this contribution demonstrates what insights into open-system magmatic process can be gained by examining isotopic heterogeneity at the inter-mineral and inter-sample scales. Particularly, this approach can track the interplay between crystallisation and assimilation and probe the physical nature of the assimilant.
The granitoid Wallundry Suite is dominated by the Middledale Gabbroic Diorite (MGD), which hosts the TEMORA 1 and 2 zircon reference materials (Wormald, 1993;Black et al., 2003;Black et al., 2004). The age of the MGD is well established by various studies. Black et al. (2004) showed that although the accepted ages of the TEMORA 1 and 2 zircons ( 206 Pb/ 238 U ages of 416.8 ±1.1 Ma and 416.8 ±1.3 Ma, respectively, as determined by isotope dilution thermal ionisation mass spectrometry, ID-TIMS) are identical. Iles et al. (2015) confirmed the homogeneity of U-Pb ages (at the precision of ionmicroprobe techniques) across the pluton. The best estimate TEMORA 2 age is the chemical abrasion (CA)-ID-TIMS 206 Pb/ 238 U age of 417.68 ±0.15 Ma reported by Ickert et al. (2015). The MGD exhibits isotopic heterogeneity for both O (Black et al., 2004;Allen and Campbell, 2012;Hervé et al., 2013;Yakymchuk et al., 2013) and Hf (Woodhead and Hergt, 2005;Iles et al., 2015). Iles et al. (2015) documented zircon O and Hf isotope compositions and covariation in the MGD consistent with variable crustal contamination of a basaltic parent magma, with major and trace element trends attesting to fractional crystallisation. Iles et al. (2015) demonstrated that Zr saturation is reached late in the crystallisation history and the gabbroic samples crystallise zircon in interstitial melt pools aided by late-stage alteration of hornblende. Whereas most previous research focussed on the MGD, its zircon isotope characteristics and zircon formation, here the petrogenesis of the whole Wallundry Suite is investigated by combining petrographic and geochemical observations with Sr-Nd-Hf isotope data. Bulk-rock and mineral isotope analyses serve to explore the ways that assimilation and fractional crystallisation operated to produce petrographic and geochemical variations.
The MGD is dominated by subequal amounts of euhedral to anhedral plagioclase and hornblende. Minor amounts of clinopyroxene, magnetite, ilmenite, apatite, titanite and quartz are observed in varying proportions across the pluton. Olivine and orthopyroxene are rarely seen in petrographic thin section. Biotite is also rare, as it has been replaced by chlorite in most samples. Actinolite and, to a lesser extent, epidote are alteration products of Map of the Wallundry Suite, the regional geology and the broad geological setting with southeastern Australia (modified after Figure 1 of Iles et al., 2015). (A) Map of the orogenic zones of southeastern Australia, including the Lachlan Fold Belt, which is divided into the Western, Central and Eastern Lachlan by the Stawell-Ararat Fault Zone (1), the Governor Fault Zone (2) and the Gilmore Fault Zone (3). The Parkes Thrust (4) and the Coolac-Narromine Fault Zone (5) bound the region outlined in red, which is expanded in (B) to show the regional geological setting of the Wallundry Suite. (C) Map of the Wallundry Suite, showing the sample sites from this study (A-P, S, U, Y, WG, WQM and MG) and those from Wormald (1993;numbered sites). The TEMORA 1, two and three sample sites are located to the north of Site A, the east of Site B and the north-east of Site A, respectively, although these boulders are not in situ.
The MGD is dominantly medium-grained (2-4 mm), but variations occur such as the contrast between the fine-grained TEMORA 1 host rocks and the coarser-grained host of TEMORA 2 (Black et al., 2004). Samples along a central NNW trending ridge are fine-medium-grained (1-2 mm), whereas samples located topographically lower to the west and east are medium-coarsegrained (2-4 mm). The north of the pluton is generally more mafic than the south, where plagioclase is typically euhedral and has a higher modal abundance. Figure 2 illustrates examples of the textural and mineralogical variations observed in thin sections. The texture of hornblende is also spatially variable. Patches of hornblende that enclose small plagioclase crystals (<1.5 mm) and even smaller clinopyroxene and magnetite crystals (<0.5 mm) become increasingly prominent southwards along the central ridge and away from it into the coarser-grained western and eastern samples. These usually consist of one or two large hornblende crystals (1-6 mm varying in proportion to the average grainsize of the sample) and smaller (~1 mm) crystals of hornblende, actinolite, biotite and chlorite aggregated to form the patch. Textural relationships and mineral chemistry suggest that pyroxenes were the dominant earlycrystallising phases, followed by the onset of plagioclase accumulation and later the crystallisation of hornblende (which also partially replaced pyroxene) then biotite (Iles et al., 2015).
The WQM, located to the southwest of the MGD (Figure 1), is fine-medium-grained (~1 mm average grainsize). It exhibits an This fine-medium-grained sample (from the ridge, approximately 600 m south of A contains hornblende grains that enclose partially altered clinopyroxene (cpx) and plagioclase (plag). (D) Chlorite (chl) has replaced part of a large hornblende that has grown around multiple plagioclase crystals.

Frontiers in Earth Science
frontiersin.org 04 equigranular fabric of interlocking, quartz, plagioclase and K-feldspar, with minor hornblende, biotite, magnetite and titanite. To the north of the WQM, the WG has a medium-grained, equigranular appearance (most grains are~2 mm wide). In addition to quartz, plagioclase and K-feldspar, the mineral assemblage consists of magnetite, biotite, chlorite and minor apatite, titanite, epidote and calcite. Examination of petrographic thin section reveals pseudomorphs of biotite and chlorite after hornblende and iddingsite after olivine. A felsic dyke (sample U) that intrudes the southern portion of the MGD is petrographically identical to the WG. In both WQM and WG, there is partial alteration of plagioclase and K-feldspar. The MG is a felsic coarse-grained (3-6 mm, with some K-feldspar up to 15 mm) alkali feldspar granite. Only a highly weathered outcrop could be accessed for this study; however, Wormald (1993) reports that biotite and both primary and secondary muscovite are present, along with zircon, apatite, iddingsite and rutile as accessory phases. Furthermore, 2-4 cm cavities infilled with tourmaline and quartz are also observed in the MG.

Analytical methods
Splits of the samples collected for this study ( Figure 1) were shared for sample preparation at Geoscience Australia (GA) and the School of Geography, Earth and Atmospheric Sciences at the University of Melbourne (UM). Zircons were extracted from a subset of the samples at GA. For all samples except Sites S, U, and MG, X-ray fluorescence spectroscopy (XRF) at GA was used to acquire bulk-rock major and trace element data. Trace element analyses of these samples were also conducted using an Agilent 7500ce quadrupole inductivelycoupled plasma mass spectrometer (ICPMS) at GA. Duplicate trace element data were acquired from the UM splits using an Agilent 7700x quadrupole ICPMS. Based on an assessment of inter-laboratory bias and analytical constraints such as detection limits for some elements (Supplementary Material) the most reliable data have been compiled for use in figures and the discussion. The UM splits were used for bulkrock isotope analysis and mineral separate geochemical analyses.
The zircon U-Pb geochronology and O-Hf isotope ratio determinations presented in this study were conducted in the same analytical sessions as the study of Iles et al. (2015), wherein the analytical methods are reported in detail. Briefly: zircon U-Pb geochronology was performed using the SHRIMP II at GA; the O isotope compositions of zircons were determined using the SHRIMP II at the Research School of Earth Sciences, Australian National University (ANU); zircon Hf isotope analyses was performed using laser ablation multicollector inductively coupled plasma mass spectrometry (LA-MC-ICPMS) at UM; cathodoluminescence images of zircons were used to guide analytical spot selection and ensure that, following U-Pb analyses, the O isotope then Hf isotope analyses were performed in the same domains for each technique.
Bulk-rock sample dissolution and subsequent element separation chemistry was conducted at UM according to the methods described in Iles et al. (2020); however, the dissolution procedure was simplified in the first run as the samples were prepared without isotopic tracers (i.e., Sr, Nd and Hf isotopes were analysed 'unspiked'). In addition to the bulk-rock powders, mineral separates of two Wallundry Suite samples (Site J from the MGD and the WG sample) were selected for isotopic analysis. Mineral separation sought to obtain hornblende from Site J, mafic aggregates (clots dominated by chlorite and biotite, likely pseudomorphs after amphibole) from the WG sample, heavy weak-to non-magnetic accessories (zircon removed from the fraction, leaving titanite and epidote) from both, and ilmenite from Site J. Having been crushed, separated by conventional magnetic and heavy liquid techniques and hand-picked, the grains were cleaned by heating in dilute nitric acid and rinsed multiple times in purified water before dissolution. The grains were then prepared for analysis in same way for the powdered samples. Trace element analyses of splits of the dissolved separates were also performed, using the quadrupole ICPMS at UM.

Major and trace element geochemistry
When SiO 2 is plotted against MgO, it is clear that, although SiO 2 overall increases as MgO decreases, the majority of samples show little variation in SiO 2 (45.5-49.8 wt%) compared with the changes in MgO content (~3 and 6 wt%). Thus, MgO content provides the better proxy for progressive crystallisation and magma evolution. A general observation from plots of major and trace elements versus MgO is that the data form coherent trends and the oxides of mobile elements (K 2 O and Na 2 O) retain this magmatic coherency and do not show a significant degree of scatter ( Figure 3).
Many of the major element variation diagrams show kinked trends ( Figure 3). The variation in the major element oxide contents with MgO allows the geochemical data to refine the paragenetic sequence suggested by petrographic evidence. The plot of TiO 2 vs. MgO ( Figure 3B) shows that the data define an approximately linear array between 0 and~4 wt% MgO, and that this divides into two clear groups at higher MgO. Of the two groups, those samples with higher TiO 2 are broadly co-linear with the main array, whereas the samples with lower TiO 2 appear to define a trend of increasing TiO 2 with decreasing MgO. In a number of other variation diagrams (e.g., those illustrating variations in Al 2 O 3 , CaO, Na 2 O, total Fe) there is a clear inflection point at~3 wt% MgO. Apart from the samples with the highest TiO 2 contents, the variations in TiO 2 , Fe 2 O 3 and CaO among the most mafic samples (MgO 4-6 wt%) are consistent with the Frontiers in Earth Science frontiersin.org removal of Mg-rich, Ti-poor minerals. The removal of olivine (reported as a rare occurrence) would change the SiO 2 and total Fe contents only subtly, but increase the CaO and TiO 2 content of fractionated magmas. At 4 wt% MgO, the replacement of olivine by orthopyroxene (also a rarely observed phase) and appearance of clinopyroxene would generate the inflections illustrated. Assemblages dominated by pyroxene (both clinopyroxene and orthopyroxene) would buffer the CaO content but drive the Al 2 O 3 contents to higher values. At MgO contents below~3 wt% it is apparent that plagioclase becomes the dominant fractionating phase (both CaO and Al 2 O 3 display a marked decline) but with decreasing MgO, indicating clinopyroxene removal is also important. The total Fe content also decreases at this point, reflecting the more Fe-rich nature of the clinopyroxene being removed, as well as the possible involvement of Fe-bearing oxides.
Magnetite is among the opaque oxides observed to occur in these lithologies of intermediate bulk-rock composition.
The above interpretation is consistent with the initial sharp decrease in Ni (not shown) in the most magnesian samples with lower TiO 2 contents. The removal of olivine dominated assemblages would be expected to result in the initial sharp decrease in Ni in the fractionated magmas. In addition, although the variations are subtle, it would also appear that Al 2 O 3 decreases in these same samples ( Figure 3C), and the behaviour of Cr (~60 ppm in one sample at MgO~6 wt%, Cr <15 ppm in all others) is consistent with the removal of Cr-bearing spinel. The samples located on an extension of the main TiO 2 -MgO trend to high MgO contents (beyond~4 wt%) and elevated TiO 2 are interpreted here to represent the accumulation of pyroxenerich assemblages at the point where olivine is replaced by orthopyroxene (see Figure 3B). This is consistent with the observation that the highest clinopyroxene modes and the rare preservation of orthopyroxene are found in this group of samples.
From~3 wt% to 0 wt% MgO, the pattern of increasing Na 2 O changes with values decreasing in the most felsic samples. This is

FIGURE 3
Variation diagrams for all Wallundry Suite samples, with elemental oxides plotted against MgO. (A) SiO 2 varies little for most of the MGD samples, but increases markedly with decreasing MgO for the felsic samples. (B) TiO 2 shows two trends that intersect at~4 wt% MgO. (C) With decreasing MgO content, Al 2 O 3 begins to increase sharply at~4 wt% MgO, then decreases from~3 wt% MgO. As for Al 2 O 3 , (D), (E) and (F) show that~3 wt% MgO marks an inflection point for many major elements (Ca, Na, Fe). The trend in K 2 O is similar to that of SiO 2 . The five samples with highest TiO 2 are shown using a different symbol (diamonds with a heavier red outline) and are believed to represent pyroxene-rich cumulate rocks. The arrows on panel (D) illustrate the general shape of a possible liquid line of descent and are labelled with the main phases believed to control the fractional crystallisation assemblage (see text for details).

Frontiers in Earth Science
frontiersin.org 06 consistent with the increasing Na content of the plagioclase as fractional crystallisation proceeds and is supported by the variation in plagioclase composition reported by Wormald (1993). In the Napoor sample 35 (MGD, 2.92 wt% Na 2 O, 4.45 wt% MgO) cumulate plagioclase compositions range from An 6o (core) to An 53 (rim), whereas sample 312 (WQM, 5.05 wt% Na 2 O, 0.50 wt% MgO) contains cores from An 49 to An 16 and rims from An 36 to An 8 .
Given the fluid mobility of large ion lithophile elements (LILEs) and the petrographic evidence of alteration, scatter of these elements on variation diagrams might be expected; however, variation diagrams for CaO, Na 2 O (mobile elements) as well as K 2 O (and other LILEs) versus MgO display coherency, and scatter that might have indicated open-system alteration is lacking ( Figure 4). Barium and Rb behave in much the same way as the K 2 O and are the least scattered LILEs; Cs is generally similar though it preserves more scatter; and, Sr displays a trend very similar to that of CaO with respect to MgO, with a maximum concentration (610 ppm, sample J) and an inflection point at~3 wt % MgO. Thus, concentrations are preserved on a hand specimen scale despite the local, grain scale redistribution by alteration discussed in Iles et al. (2015). Figure 4E shows that the majority of the MGD samples have broadly low Zr (88-132 ppm), but the Zr concentration steps up in the most felsic samples (WQM, WG and sample J from the MGD). The variation in Zr concentration in the more felsic samples can be further examined by plotting the data against SiO 2 ( Figure 4F). This shows that the Wallundry Suite reaches a maximum Zr concentration in the WQM (sample 169 of Wormald, 1993) then decreases with increasing SiO 2 .

Zircon isotope geochemistry
In addition to the U-Pb ages and O-Hf isotope data for zircons from the four samples reported by Iles et al. (2015), equivalent analyses (B and C) Rb and Ba show the least scatter, and have trends similar to K 2 O. (D) Cs is the most scattered of the LILEs and its trend is similar to that of Rb, but enrichment occurs only in one sample (WQM). (E and F) Zr concentration plotted against MgO and SiO 2 respectively. Plotting against MgO shows the uniformity of Zr for a range of mafic compositions (all bar one of the MGD samples). Plotting against SiO 2 and including the data of Wormald (1993) and Black et al. (2004) shows the detail of the variation of Zr in the more felsic samples of the Wallundry Suite.

Frontiers in Earth Science
frontiersin.org of zircons from one sample each of the WQM and WG are reported here. Mean 206 Pb/ 238 U ages of 413.9 ± 2.0 Ma (n = 29) and 412.9 ± 2.0 Ma (n = 29) were determined for zircons from WG and WQM, respectively, which are within uncertainty (quoted at the 95% confidence limits) of the four MGD age determinations reported by Iles et al. (2015). The age for WQM is, however, slightly outside of uncertainty with respect to the accepted ages of the TEMORA 1 and 2 (quoted above; Black et al., 2003;Black et al., 2004). Additionally, a minor inherited population of morphologically distinct zircons ("large grains") with a mean 206 Pb/ 238 U age of 427.4 ± 4.3 Ma (n = 7) is found in the WQM sample. In contrast to the "small grains" (mostly 50-100 µm long) that define the crystallisation age, the "large grains" are mostly 100-200 µm long and have a uniform dark appearance in CL and higher U concentrations and Th/U ratios (1239-2760 ppm and Th/U = 0.83-1.68 compared to 143-1040 ppm and Th/U =0.32-0.74).
In situ O isotope analysis of zircons from WG yield mean δ 18 O values of 8.33 ± 0.23 (1s.d.) for a coherent group (n = 28) and 8.30 ± 0.30 if an anomalously low outlier (δ 18 O = 7.27 ± 0.12, 1σ) is included. Analysis of zircons from WQM yields a mean δ 18 O value of 8.36 ± 0.28 for a coherent group (n = 26) and 8.25 ± 0.54 if two outliers (δ 18 O = 5.99 ± 0.08 and 7.66 ± 0.06, 1σ) are included. Both samples are thus within uncertainty of the coarse-grained MGD samples from Iles et al. (2015). Interestingly, analyses from the 'large grain' group yield δ 18 O values showing some of the least deviations from the mean for the sample and, as a consequence, removal of the 'large grain' group does not significantly alter the mean for the WQM ( Figure 5). The δ 18 O values for the two age groups present in this sample thus are indistinguishable.
For Hf isotope compositions, because sample sizes are smaller and the is greater uncertainty on individual analyses relative to scatter, a more robust measure of central tendency has been used. A weighted mean 176 Hf/ 177 Hf ratio of 0.282761 ± 0.000021 (95% confidence about the Tukey's biweight mean, n = 17) is preserved in the zircons from WG. This includes one anomalously high outlier, but such values are down-weighted in Tukey's biweight mean. As noted for the O isotope values, the mean 176 Hf/ 177 Hf ratios of the 'large grain' and 'small grain' groups of WQM are indistinguishable (each 0.28273 ± 0.00010 from five analyses). If the two groups are combined, they yield weighted mean 176 Hf/ 177 Hf ratio of 0.282726 ± 0.000039 (n = 10). The zircons of the 'large grain' group have higher 176 Lu/ 177 Hf ratios than those of the 'small grain' group; therefore, the 176 Hf/ 177 Hf ratios of the two groups diverge slightly after age correction but are still within uncertainty (0.28271 ± 0.00010 for the 'small grains' and 0.282689 ± 0.000099 for the 'large grains'). The O and Hf isotope data discussed above are illustrated in Figure 5. These values correspond to εHf i of 7.7 ± 0.7 and 6.2 ± 1.4 (95% confidence) for WG and WQM, respectively, using the 416.8 Ma age of Black et al. (2004). Furthermore, the paired O and Hf isotope analyses (outliers excluded for this purpose) for all six zircon populations are plotted in Figure 6, demonstrating the overlap between the coarse-grained MGD samples and the WG and WQM samples in O-Hf isotope space and the contrast between these and the fine-grained MGD sample.

Bulk-rock isotope geochemistry
The results of the bulk-rock Rb-Sr, Sm-Nd and Lu-Hf isotope analyses are reported in Table 1 and Table 2. Initial isotope ratios have Frontiers in Earth Science frontiersin.org been calculated using the U-Pb zircon ages of 416.8 Ma for the Wallundry Suite (Black et al., 2004). Plots of the Sr, Nd and Hf isotope compositions of bulk-rock samples (Figure 7) show that the suite forms an array in Sr-Nd-Hf isotope space, which might be expected given the previously described trends in zircon O-Hf isotope space for Wallundry. The samples appear to form clusters within a broader Nd-Hf isotope array. Samples A, H, L and S from the MGD define a more isotopically primitive group with εHf i ranging from 7.2 to 9.2 (2 s.d. range) and εNd i ranging from 3.1 to 5.1. Samples B, J, N and Y from the MGD and the WG, WQM and felsic dyke samples constitute a more evolved group (εHf i 3.4-6.4 and εNd i 0.9-2.9) with an internal trend congruent with the overall Nd-Hf trend of the suite. The most evolved sample is the MG, with an εHf i value of 2.40 ±0.53 and an εNd i value of −0.89 ±0.38. Although such clustering is not as obvious for the Sr-Nd isotope data, the same division between the two groups defined from the Nd-Hf isotope plot is nevertheless evident in Figure 7A. The evolution of bulkrock Sr-Nd-Hf isotope compositions does not strongly follow the lithological and chemical evolution of the suite. This is particularly well illustrated by the essentially identical Sr, Nd and Hf isotope compositions of MGD Site N and the WQM sample ( 87 Sr/ 86 Sr of 0.70472, εNd i of 1.7, εHf i of 4.2), despite different major element chemistry, petrography and trace element compositions (e.g., 23.5 ppm Rb and 341 ppm Sr for Site N and 121 ppm Rb and 127 ppm Sr for WQM).

Hf isotope composition of the melt phase
An important and unexpected early observation of this study is the discrepancy between the magmatic zircon Hf isotope compositions and the bulk-rock Hf isotope composition for each of Wallundry Suite samples. This is illustrated in Figure 8. In order to rule out any possible cause resulting from analytical problems, multiple repeat analyses were performed on bulk-rock samples. During these experiments it was confirmed that the observation is genuine (i.e., the phenomenon is not analytically induced) and further noted that the bulk-rock Hf compositions of different aliquots of the same sample varied significantly.
Following the approach described in Iles et al. (2020), the difference between the initial εHf values of the bulk-rock samples and the means of their magmatic zircon populations is represented by their ΔεHf bulk-zrc values. Intriguingly, although the low-temperature I-types of Iles et al. (2020) consistently exhibited positive ΔεHf bulk-zrc values, the high-temperature I-type Wallundry Suite behaves differently. Although some samples of the Wallundry Suite exhibit greater variability from aliquot to aliquot (poorer reproducibility) than others (2.0ε units for Site J compared to 0.5ε units for Site L), the bulk-rock Hf isotope measurements produce age-corrected εHf i values that are all within uncertainty of or lower than the εHf i values of the zircon population. There is no clear relationship between the ΔεHf bulk-zrc values and the chemical or petrographic character in the samples. For example, although the largest ΔεHf bulk-zrc value occurs in the fine-grained, mafic sample A, the ΔεHf bulk-zrc values are smallest among the more isotopically evolved samples, which range from mafic to felsic in composition.

Geochemistry of mineral separates
That the bulk-rocks and zircons do not preserve the same Hf isotope compositions indicates the presence of inter-mineral isotopic heterogeneity in the samples, a component of which must balance the composition of the zircon. This component would need to be less radiogenic in Hf than the bulk-rock samples. Isotopic analyses of mineral separates from two samples (Site J of the MGD and WG) explore this further, to uncover the identity of the less radiogenic component in the bulk-rock samples. The Sm-Nd and Lu-Hf isotope data are presented in Table 3. Trace element compositions have been included in the Supplementary Material.
The titanite-epidote separate from WG has MREE and HREE abundances similar to the range of published data from titanites,   Note the weak coupling between the evolution of the isotope compositions and the evolution from mafic to felsic samples.   (Black et al., 2004;Iles, 2012  The Ti concentration measured in the Fe-Ti oxide separate from sample J is~19 wt%, which is approximately two-thirds of the concentration expected for pure ilmenite, but could reflect hematite-ilmenite solid solution in the grains. Alternatively, minor inclusions have been observed in the ilmenite grains when viewed under reflected light and these may also be responsible for 'diluting' the measured Ti content. Certainly, the Ca concentration of~8000 ppm (higher than any other trace element measured in this separate) is unusual and is possibly attributed to titanite inclusions. Chalcopyrite is observed in petrographic thin section, sometimes attached to or included in ilmenite; however, the Cu concentration (39 ppm) does not suggest the presence of significant chalcopyrite in the separate. The trace element pattern for the Fe-Ti oxide separate is broadly similar to ilmenite compositions when plotted on a "spider diagram"; however, the separate is enriched in Th, La, Ce and Nd and to a lesser extent Ba, U, Sr and Y (Supplementary Material).
For both samples, the εHf i values of hornblende (biotite-chlorite pseudomorphs after primary amphibole in the case of WG) and titanite-epidote separates (where available) are bracketed between the zircon LA-ICPMS values and the lowest repeat bulk-rock * Numbers refer to the aliquot numbers used in Tables 1 and Table 2 for bulk-rock powders, and to distinguish different mineral separates.
Frontiers in Earth Science frontiersin.org values (Table 3). Similarly, the same mineral separates record εNd i within one to two ε units of the bulk-rock compositions (one biotitechlorite separate from WG, which is 1.9 εNd units lower than the lowest bulk-rock analysis, deviates the most). Crucially, one mineral separate, the Fe-Ti oxides from sample J, has Nd and Hf isotope compositions unlike any other analysis. With εHf i = −9.1 and εNd i = −13.6, the Fe-Ti oxide phases are distinctly out of equilibrium with the zircon and other melt precipitated phases, and their compositions would account for the discrepancy observed between bulk-rock and zircon Hf isotope compositions (Figure 9).

The complex geochemistry of the wallundry suite
Many features of the Wallundry Suite petrography and geochemistry indicate that fractional crystallisation was an important process in the formation of this suite: cumulusintercumulus relationships between plagioclase and amphiboles ( Figure 2D); relicts of early-formed pyroxenes ( Figure 2C); inflections and turning points in major and trace elements trends (Figures 3, 4), suggesting sequential appearance of crystallising phases. The isotopic variability exhibited by the magmatic zircon attests to an isotopically evolving melt, suggestive of crustal contamination in the petrogenetic history.
The isotopic similarity of the zircon populations from coarsegrained MGD, WG and WQM samples and the contrast with the more scattered data from the fine-grained, marginal lithology (Site A) could indicate that the dominant compositional group represent a parental magma that underwent fractional crystallisation to produce the range of samples observed, and that Site A represents a portion of the pluton affected by an additional process (such as interaction with meteoric water) that altered the isotope compositions and generated the marked scatter. This would, however, require a process that caused scatter in O-isotopic compositions and Hf-isotopic compositions in zircons. Although the presence of strongly sericitised plagioclase in sample A suggests the scatter in O isotopes could be due to alteration processes, such processes would not be expected to generate the observed scatter in the Hf isotope compositions, due to the slow diffusion and effectively immobile nature of Hf (Cherniak et al., 1997). Alternatively, the samples could be part of an assimilationfractional crystallisation trend in which the low δ 18 O and high 176 Hf/ 177 Hf of the zircons from the fine-grained sample represent a less contaminated (but heterogeneous) melt composition than those represented by the zircon populations, with the higher δ 18 O and lower 176 Hf/ 177 Hf, found in the coarser-grained samples. Either of the explanations for the diverse zircon O and Hf isotope compositions of the suite must also account for the apparent entrainment of older zircons into the WQM alone and their remarkable isotopic similarity to the dominant, younger WQM zircons. The older 'large grain' zircon group may have originated from an earlier magmatic pulse from the same source as that which produced the dominant zircon O-Hf isotope compositions of the Wallundry Suite; or their entrainment represents the contamination of the magma that formed the WQM, causing the melt to reach an isotopic composition essentially equivalent (coincidentally) to that of the zircons in the contaminant.
The newly acquired Sr-Nd-Hf bulk-rock isotope data further demonstrate the isotopic variability of the suite, and the correlation between these isotope systems adds to the evidence suggesting progressive evolution of the magma system by assimilation processes (Figure 7). It becomes apparent, however, by comparing the various sources of isotope data with the major and trace element data, that a straightforward combination of assimilation and fractional crystallisation (with isotope and trace element chemistry changing continuously with major element variations) is inadequate to explain the data. Two granitoids (WG and WQM) are isotopically similar to some of the most felsic samples of the MGD, despite being very different chemically. Sample L (52% SiO 2 ) is isotopically similar to the most mafic samples (A and H,~47% SiO 2 ). The MG seems to have evolved both isotopically and chemically, being the most felsic sample and having the lowest εNd i and εHf i values. The grouping of samples into 'mafic', 'intermediate' and 'felsic' groups in Figure 7 highlights this apparent decoupling between the evolution of Sr and Nd isotopes and the changes in major element compositions in the suite.
This is also born out in the Hf isotope data. Plotting εHf i against MgO content ( Figure 10A) shows that a wide range of samples have very similar Hf isotope compositions (εHf i around 4-5), and sample L again plots off-trend as an isotopically primitive but chemically  Figure 10B) hints at the presence of at least two trends. From the most mafic and isotopically primitive samples (A and H) a low Hf trend seems to extend through mafic and intermediate samples to the felsic MG and a high Hf trend seems to extend through sample L and into the other granitoids (WG, WQM and the felsic dyke). Alternatively, the low Hf samples could reflect an AFC trend, with the high Hf samples resulting from later fractional crystallisation of different batches along that trend. Thus, isotope data combined with major and trace element data point to different degrees of contamination in different samples, different degrees of fractional crystallisation and some decoupling between the two. A possible scenario based on this would be a two-level petrogenetic model. Firstly, assimilation accompanied by some fractional crystallisation occurs below the emplacement level, producing diverse magma batches. Secondary processing via fractional crystallisation alone might occur multiple times as these batches ascend and are emplaced. The somewhat random distribution of granitoids around the gabbro-diorite body, the discovery of a felsic dyke in the south of the MGD that is similar to the WG and hand specimens (albeit not in situ) that show mingling textures support intrusion in batches. Although such a model mostly reconciles the field observations, petrography and bulk-rock major and trace element data with the bulk-rock isotope data, it does not account for the disequilibrium between bulk-rock and magmatic zircon Hf isotope composition and the intra-sample Nd-Hf isotope heterogeneity. The highly evolved Nd and Hf isotope of the Fe-Ti oxide separate explains why bulk-rock samples would have less radiogenic Hf isotope compositions than the zircon, but also adds complexity to the petrogenesis. With εHf i = −9.1 and εNd i = −13.6, the Fe-Ti oxide must be derived from an unmelted (restitic) component of crustal material (possibly sedimentary in origin), which provides strong support for the role of assimilation, but complicating it by requiring that the contaminant is only partially melted and that the solids do not equilibrate with the mixed melt. A petrogenetic model for the Wallundry Suite must, The graph of εHf i versus Hf concentration demonstrates that there is no simple relationship between the Hf isotope composition and trace and major element geochemistry. This graph hints at the presence of more than one trend in the evolution of the suite. There is a main "low Hf" trend from mafic, isotopically primitive samples through to the felsic, isotopically evolved Mayfair Granite, but there are also samples with higher Hf concentrations that suggest additional magma evolution pathways.

FIGURE 11
Magma-melt relationships for the simple mixing of basaltic melt with bulk incorporation of a partially melted sediment. (A) Relationships between mixing endmembers in SiO 2 -Hf isotope space. Overall the mixture is between basalt and the sediment ("Source"). The sediment has a melt component that can mix with the basaltic melt. Two different possible melts (connected to the bulk-sediment by the "source-melt lines") are shown along with the two different "melt mixing lines" that show the composition of the mixed melt at increments that match those of the overall magma mixing. The compositions of the sedimentary components are modelled by applying the disequilibrium melting approach of Iles et al. (2018) to biotite dehydration (Vielzeuf and Holloway, 1988) of a metapelite (Bea and Montero, 1999). These relationships generate negative ΔεHf i magma-melt values (mostly) because the sediment is much less radiogenic than the basalt and the mixed magma contains both solid (contaminant restite) and melt components of the sediment whereas the mixed melt only contains the melt component. Note that for very high fractions of sediment in the magma mixture of the second melt case (ΔεHf source-melt = −0.2, 41% melt) the ΔεHf i magma-melt values are very small negative values. (B) For the two different melt cases, the relationship between the Hf isotope composition of the magma and its melt is shown, relative to a 1:1 line. Note that the maximum differences between magma and melt are at the intermediate degrees of mixing.

Frontiers in Earth Science
frontiersin.org therefore, account for the complex combination of assimilation and fractional crystallisation represented by the trends in major element, trace element and isotope data, the partially molten nature of the contaminant and the mismatch between the bulk and zircon Hf isotope data. This is well beyond the assumptions and framework of classical DePaolo (1981) AFC or EC-AFC Spera and Bohrson, 2001). The diverse options that can be explored with the Magma Chamber Simulator (Bohrson et al., 2014;Bohrson et al., 2020;Heinonen et al., 2020) could potentially address some aspects of the Wallundry Suite petrogenesis; however, MCS does not account for the behaviour of zircon in the magma and wallrock subsystems.
Crucially, the petrogenetic model must address why the Wallundry Suite samples have negative ΔεHf bulk-zrc values-that is, why would the melt from which zircon crystallises have more radiogenic Hf than the magma dominated by earlier crystallised minerals given that assimilation of common continental crust lithologies (with generally much less radiogenic Hf than the most primitive samples) should cause melt evolution from more radiogenic to less radiogenic Hf isotope compositions? Keeping in mind the evidence (Nd and Hf isotope compositions of the Fe-Ti oxide) that unmelted mature material from continental crustal rocks is present in the samples, this question can be examined through the simplified model illustrated in Figure 11. Setting aside the role of fractional crystallisation, assimilation of partially molten sediment by a basaltic magma is simplified to the three-component mixing of basaltic melt, sedimentary melt and sedimentary restite. The proportion of sedimentary restite relative to sedimentary melt is fixed by the degree of melting, but the ratio of basalt to sedimentary contaminant (melt + restite) is a free parameter. Previous work (Tang et al., 2014;Iles et al., 2017;Iles et al., 2018;Wang et al., 2018;Wolf et al., 2019) has demonstrated that partial melting in the crust may generate melt-restite Hf isotopic disequilibrium; therefore, the sedimentary melt and unmelted components need not have the same isotope composition. With the exception of extreme positive ΔεHf source-melt values, a sedimentary restite would be far less radiogenic than basaltic melt, and thus also less radiogenic than a mix between sediment melt and basalt; therefore, the εHf value of a mixed magma (basalt + sedimentary restite + sedimentary melt) would be lower than the εHf value of the mixed melt (basaltic + sedimentary melt). This holds even at ΔεHf source-melt~0 , and negative ΔεHf magma-melt values will arise also for positive ΔεHf source-melt values, as long as the proportion of basaltic melt in the mixture is high enough (the basalt fraction required is dependent on the Hf concentrations of the components and the positive ΔεHf source-melt value).
This simplified case provides a conceptual understanding of the way the bulk-rock-zircon Hf isotope discrepancy arises and is tied to the assimilation of a partially molten crustal rock. It does not, however, account for the role of fractional crystallisation. Furthermore, the bulk-rock-zircon discrepancy in the granitic members of the suite is outside of what can be accounted for by this simple mixing because they have higher SiO 2 contents than common sediments (Wyborn and Chappell, 1983;Vielzeuf and Holloway, 1988;Patiño Douce and Johnston, 1991;Turner et al., 1993;Vielzeuf and Montel, 1994). The felsic nature and high Hf concentrations of the WQM and WG, in particular, require that fractional crystallisation be integrated into the model. The following section discusses some additional considerations and constraints on these complexities.

Constraints of assimilation and fractional crystallisation
Considering the possibility that the mixing illustrated in Figure 11 may be accompanied by crystallisation and, by extension, liquidcrystal sorting, there are various factors that would play a role in the production of the observed bulk-rock vs. zircon Hf isotopic disequilibrium. These are considered here.
One factor already considered in Figure 11 is the melt fraction (f = 0-1) of the contaminant (sedimentary rock in the mixing example). As illustrated, with increasing degree on melting, the degree of magmamelt isotopic disequilibrium decreases (although this will be complicated by disequilibrium in contaminant melting). In the example, for f = 0.41 the maximum ΔεHf i magma-melt value is~5, whereas the largest ΔεHf i bulk-zrc values recorded by the samples arẽ 3, suggesting a similar or higher degree of melting of the contaminant is required for the Wallundry Suite. A ratio, R, defines the mass of melt crystallised incrementally during the process relative to the mass of sediment incrementally assimilated, and this would have an important impact on these relationships.
Crystallisation will affect the melt compositional evolution, but the magma trend in Figure 11 will not change unless crystals are removed from the magma. Without crystal segregation, the magma will increase in mass and crystallinity and decrease in melt mass as crystallisation progresses (unless R and f are such that more melt is added from the contaminant than is crystallised at each increment). The melt composition will be affected by the compositions (Hf contents and Hf isotope ratios) of the assimilated melt and the primary magma, as well as by the partitioning behaviour of Hf between the melt and the crystals. Additionally, the partitioning is affected by phase variations in the crystal assemblage and whether or not equilibrium is achieved between melt and unsegregated early crystals (fractional crystallisation of the melt phase could occur via chemical isolation of earlier-formed grains rather than crystal settling).
Although a series of magmas could be produced through assimilation-crystallisation without crystal segregation (and this would be analogous to the mixing relations shown in Figure 11), batches of this magma evolution could subsequently undergo partial segregation. Compositions between the magma trend and the melt trend would represent magma batches that have undergone some degree of melt segregation (simple liquid-crystal fractionation). This would lower the magnitude of melt-magma disequilibrium, since magmas would become more melt rich. Samples would then plot between the 1:1 line and the mixing trends in Figure 11B. Alternatively, any composition along the magma trend could be selected to become the parent magma for fractional crystallisation at the emplacement level to yield chemically evolved compositions (with similar implications for the positions of samples in Figure 11B). Similarly, crystal segregation concurrent with crystallisation adds another layer of complexity to the system evolution and the possible relationships between magma and melt major and trace element and isotope compositions.
Some more specific considerations can also be made. Given the compositions of the Wallundry samples, the primary melt is assumed to be a hot, basaltic melt, in which Zr and Hf bulk partition coefficients (D Zr and D Hf ) would be initially below unity (Zr and Hf are incompatible). The AFC process causes Zr concentration in the melt to increase while D Zr <1, but this would eventually lead to Zr saturation, causing the onset of zircon crystallisation and a change in Frontiers in Earth Science frontiersin.org Zr and Hf partitioning; therefore, it is assumed that both D Zr and D Hf should increase to values greater than one at some point during the AFC process. Zircon crystallisation would act to deplete the melt in Hf. An important effect of the decreasing Hf content is that the melt in the contaminant has more weight in the Hf isotope budget, causing progressively larger decreases in the εHf value of the mixed melt (until the mixed melt approaches the contaminant melt in isotope composition, when the decreases in εHf value become smaller). This highlights both the importance of the role of zircon in future attempts to model this system and the possibility for nuances that cannot be captured in the simple mixing model. As noted previously the more granitic members of the suite (WG, WQM, the felsic dyke & MG) are too felsic for the mixing model in Figure 11 to account for their compositions. Furthermore, Figure 10 suggests some may belong to a different trend than the MGD, or involve additional fractionation. Two points should be highlighted. Firstly, the magmatic zircon in WQM is much more heterogeneous in Hf isotope composition than the zircon in WG (and in most MGD samples) and is actually within uncertainty of the WG zircon (at the 95% confidence level). Such heterogeneity implies that the sample records a range of melt compositions, and thus should not be modelled as a single contaminated magma and coexisting melt. Secondly, The MG is likely the product of more extensive fractional crystallisation of a contaminated magma or melt, which could be consistent with the low Nd content (~3.4 ppm), high SiO 2 (~75 wt%) and extreme Rb/ Sr (26.7;Wormald, 1993).
The oxygen isotope data may provide additional constraints upon the system's evolution. As a first approximation, the oxygen isotope data could be modelled as resulting from binary mixing of mantle-derived and sedimentary endmembers. Assuming a δ 18 O value of 5.5‰ for basaltic magma (Eiler, 2001;Valley et al., 2005) and δ 18 O = 13.6 ±1.6‰ (2s.d.) for the sedimentary contaminant (the average of metasedimentary rocks associated with the Cooma Granodiorite; Munksgaard, 1988), the proportion of contamination required to produce the oxygen isotope compositions of the bulk-rock samples can be estimated. The oxygen isotope compositions for magmatic zircon from the Wallundry suite are shown in Table 4, as are the corresponding bulk-rock δ 18 O values, which range from 7.87 ±0.12‰ to 10.20 ±0.20‰ (2σ external precision; see Table  caption for details). As shown in Table 4, the calculated proportion of contaminant involved each sample ranges from 29 ±17% for sample A to 58 ±25% for the WQM sample. It is expected that higher degrees of assimilation will be calculated assuming binary mixing compared with AFC (Taylor, 1980;Aitcheson and Forrest, 1994); therefore, the calculated amounts of assimilation are likely overestimates, particularly for the most evolved samples.

Conclusion
Zircon O-Hf and bulk-rock Sr-Nd-Hf isotope data indicate opensystem magmatic processes in the formation of the Wallundry Suite. This magmatic system evolved via assimilation and liquid-crystal fractionation processes (i.e., fractional crystallisation and variations thereof) with decoupling between the two at some stages.
Field observations of the Wallundry Suite indicate that it was constructed via multiple magma injections. Petrographic observations in concert with major and trace element data indicate the involvement of fractional crystallisation in generating the different members of the suite (MGD, WQM, WG and MG) and the inter-sample heterogeneity (particularly in the MGD). The zircon O-Hf data demonstrate that the suite underwent some degree of contamination during its formation; however, these data are only weakly related to the lithological and chemical character of the samples. Likewise, the newly acquired bulkrock Sr-Nd-Hf isotope data indicate different degrees of assimilation of crustal material by a basaltic parent melt and different degrees of fractional crystallisation together produce lithologically similar samples with diverse isotopic compositions, as well as isotopically similar samples with distinct lithologies and chemical compositions.
Comparison of bulk-rock and magmatic zircon Hf isotope compositions permit deeper insights into the magmatic system. The data reveal isotopic disequilibrium that requires that the samples were mixtures of isotopically diverse solids and melts. Particularly, the negative ΔεHf magma-melt values in the samples led our investigation to uncover details about assimilation that could not be inferred from O-Hf-Sr-Nd-Hf trends. Isotopic analyses of mineral separates demonstrate that the more isotopically evolved component required to account for this disequilibrium between melt (recorded by zircon) and magma (captured by the bulk-rock) was provided by Fe-Ti oxide crystals with Nd-Hf isotope compositions that can only be explained as unmelted mature crustal material. Thus, the isotope data have indicated that the material assimilated by the basaltic parent was only partially melted. This study illustrates that combined assimilation and fractional crystallisation (with variable degrees of coupling/decoupling) can give rise to highly complex geochemical relationships, especially considering the possibility that the material being assimilated need neither be completely melted nor in isotopic equilibrium. While this potential complexity makes open system processes challenging probe geochemically, we demonstrate that combining bulk-rock, mineral separate and in situ zircon isotope data can disentangle this complexity. The comparison of bulk rock and magmatic zircon Hf isotope composition is particularly useful in this regard.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.