Incremental Growth of Layered Mafic-Ultramafic Intrusions Through Melt Replenishment Into a Crystal Mush Zone Traced by Fe-Hf Isotope Systematics

Layered mafic-ultramafic intrusions (LMI) are among the largest igneous bodies on Earth, and represent aggregations of large volumes of mantle- and some crustal-derived melts. Melts are emplaced over time-intervals of less than 1 million years, predominantly through multiple pulses of injections into pre-existing melt-crystal slurries. The dynamic interaction of physical processes, including density-driven separation and mixing of different components, within a solidifying magma chamber leads to such extreme chemical diversity between cumulate rock units that no unified model currently explains all aspects of the genesis of these intrusions. Here we present whole-rock stable Fe isotope data (expressed in ‰ variations as δ57Fe relative to IRMM-014) for samples of drill core taken from the stratified paleo-magma chamber of the Upper Zone of the late-Archean Windimurra Igneous Complex, Western Australia. Variations from near chondritic (δ57Fe ∼ + 0‰) to heavy (δ57Fe∼ + 0.2‰) values show a co-variation with initial radiogenic Hf isotope data that is unique to the Windimurra Upper Zone. The systematic isotopic variations from the roof to the base of the Upper Zone are best explained by an intricate sequence of events that included fractional crystallization and physical mixing. We propose that melt freshly sourced from the mantle was injected into and inflated a pre-existing crystal-melt mush comprising the upper Middle Zone. Re-establishment of crystal layering after replenishment introduced a chemical stratification with the formation of what became the Upper Zone and a crystal-interstitial melt ratio decreasing from roof to base. Basal, vanadiferous magnetitite horizons crystallized from Fe out of the liquid-in-residence. Variable degrees of perturbation and chaotic stirring of crystals with imperfect mixing of new and old components was followed by rapid crystal settling and subsequent cumulate stratification. Such melt rejuvenation, proposed here to be the cause for a newly established Upper Zone, leaves no unique petrologic fingerprint and forges a range in hybrid magma compositions throughout the Upper Zone rather than one single parental melt. This incremental growth of the intrusion can explain not only the observed coupled Fe-Hf isotope systematics, but also mineral disequilibria and cryptic layering in layered intrusions. In the context of incremental pluton assembly, two-component mixing paired with source heterogeneity can, at least in parts, potentially explain zircon-hosted Hf isotope heterogeneity so often observed in larger magmatic bodies.


INTRODUCTION
Layered mafic-ultramafic intrusions (LMI) are among the largest igneous bodies on Earth (Cawthorn, 1996;Ernst and Buchan, 2001) and have been subject to detailed investigations for decades, yet their petrogenesis remains a matter of ongoing debate. Most LMI are emplaced into the crust via pulses of mantle-derived magma over a relatively short period of time (Zeh et al., 2015), the magma considered by many to be the result of partial melting within an uprising mantle plume at the base of cratonic, proto-lithospheric roots (Hatton, 1995;Buchan, 1997, 2001;Nebel et al., 2013a). Subdivisions of the magmatic body into stratiform rock packages reflect a broad succession from an ultramafic base towards an evolved roof section, as represented for example, by the Bushveld-, Windimurra-, Stillwater-, and Duluth complexes, and the Sept Ile Intrusion.
The stratiform packages, or "zones, " are characterized by rhythmic or intermittent modal layering of mineral abundances, which relate to complex magma chamber processes. These processes can include, but are not limited to: (i) physical crystalliquid separation through gravity settling (Naslund et al., 1991); (ii) compaction and expulsion, or filter pressing (Berger et al., 2017), often in conjunction with (iii) convection, including double-diffusive convection (Cawthorn and McCarthy, 1980;McBirney, 1995;Bons et al., 2015); (iv) melt recharge and mixing (Kruger, 2005;Yuan et al., 2017), or not (Tegner et al., 2006); (v) liquid immiscibility (VanTongeren and Mathez, 2012;Charlier et al., 2013;Vukmanovic et al., 2018) and/or accompanied by (vi) crustal contamination. Notable is that most minerals have gradual compositional zoning, which has typically been related to growth from late-stage liquids trapped within inter-grain space (McBirney, 1995). The bulk rock composition of each layer thus mirrors a combination of modal mineral abundances and interstitial, late-stage liquids, which are not necessarily in chemical equilibrium. Additional petrologic complexities are revealed by cryptic layering that is manifest in trace element variations within mineral species both within and between layers. This trace element heterogeneity may be consequent to melt rejuvenation by magmatic recharge involving two liquids that are similar in their major element systematics (Tanner et al., 2013). Due to the cumulate nature of the layered series, and often unknown mineral-melt ratios (cumulate phases vs. interstitial liquid), bulk rock compositions are also often deemed to be of limited use for modeling crystal fractionation as they do not relate to a simple parental liquid composition and cannot straightforwardly be related to a liquid line of descent (Morse, 1976). From the sum of these complex observations, it is also clear that most minerals are not in equilibrium with each other, both compositionally (Dunham and Wadsworth, 1978;Czamanske and Loferski, 1996) and isotopically (Chen et al., 2018), rendering parental liquid calculations from mineral compositions specious.
Compositional variations in major and trace elements of both cumulus and intercumulus minerals may also be accompanied by large variations in radiogenic isotope ratios, such as 87 Sr/ 86 Sr (Kruger, 1994;Maier et al., 2000;Yang et al., 2013). This isotopic heterogeneity indicates some degree of magma mixing either through magma recharge by isotopically diverse magmas or by late-stage fluid-wall rock interaction and associated contamination. These external factors must operate in addition to, and in concert with, internal magma chamber processes, such that LMI are usually treated as open systems. To this end, the sum of petrologic and geochemical observations has led to a variety of models that can explain some individual features of intrusions, but it seems that no single model can capture the complexity of all these magmatic bodies.
Here, we endeavor to constrain the processes that produced the chemical diversity found within the Upper Zone of the ca. 2.8 Ga old Windimurra Igneous Complex of the Yilgarn Craton, Western Australia (Ivanic et al., 2018). We use stable Fe isotopes measured in Upper Zone rocks in order to constrain those magmatic process that led to the formation of the layering. A unique feature of the Windimurra Igneous Complex is that the Upper Zone shows apparent initial Hf isotope heterogeneity among whole-rock samples that varies with intrusion depth. In addition, mineral separates show extreme isotope signatures with ε Hf (2.8Ga) >+20 (Nebel et al., 2013a). These variations have been interpreted as resulting from the influx of a new magma batch into a pre-existing crystal mush with an apparent two-component mixing within the Upper Zone.
We present whole-rock stable Fe isotope data for samples taken from within the Upper Zone, sampled via drill core material. Four cores collectively sample several hundreds of meters of the Upper Zone, which has characteristic layering on the cm scale, including repetitive, intermittent crystal cumulate horizons. These horizons are mainly troctolitic to gabbroic in composition, with locally concentrated Fe-Ti oxides, predominantly found at the base of the Upper Zone. A detailed description of the rocks is given in Nebel et al. (2013b).
Previous Fe isotope studies on other LMI revealed complex Fe isotope systematics in neighbouring minerals within cumulate layered intrusive rocks Liu et al., 2014), yet without any systematic variations, and even in the case of the Bushveld Complex, without notable isotope fractionation (Bilenker et al., 2017). Here, we measured Fe isotope compositions in cumulate rocks with known variable whole-rock Hf isotope systematics that relate to magma mixing, which vary with intrusion depth (Nebel et al., 2013a). We propose a hypothesis of layered intrusion magmatic recharge based on the relationship between stable Fe and radiogenic Hf isotopes that, in tandem, can explain some cryptic features of layered intrusions.

Geologic Background and Sample Description
The Windimurra Igneous Complex is a layered mafic-ultramafic intrusion, located in the central Youanmi Terrane of the Western Australian Yilgarn Craton (Figure 1; Ivanic et al., 2010;Wingate et al., 2012). It is composed of several zones (Mathison and Booth, 1990;Mathison et al., 1991;Ivanic et al., 2018) that broadly resemble those of other layered intrusions and most likely represent multiple injections of juvenile, mantle-derived melts into larger magma chambers (Nebel et al., 2013b). A detailed description of the stratigraphy of the Windimurra Igneous Complex is given in Ivanic et al. (2018). In brief, the Complex is similar in size to the Duluth Complex (United States) and has a thickness, determined from seismic imaging, of about 11 km (Ivanic et al., 2018). The base of the Complex is formed by an Ultramafic Zone (∼3 km thickness) with, in parts, serpentinized olivine-chromite horizons, overlain by a Lower Zone (∼5 km) composed of leuco-gabbro-norites and olivinegabbros. The succeeding Middle Zone (∼1.5 km) contains abundant gabbros with the first appearance of Ti-Fe oxides. The Upper Zone (∼ 1 km) of the Complex, the focus of this study, is composed of a series of layers that contain dominantly either feldspar (noritic), olivine-pyroxene (gabbroic) or Fe-Ti-oxides (magnetitic). However, each layer is a continuum between these end-member compositions with different proportions of olivine, pyroxene, plagioclase, and magnetite-ilmenite. Fresh drill core samples in this study represent various macroscopic horizons, which formed as mineral cumulate layers (Nebel et al., 2013b). In Figure 2, Harker diagrams illustrate the variability in major elements in the Upper Zone cumulates, and also show a grouping of samples into three distinct subsets of samples. The grouping is based on the respective distributions of samples in the Harker plots, which reflect the modal mineralogy of the rocks, dominated by either olivine and pyroxene, feldspar, or Fe-Ti-oxides. Drill cores MNDD001-MNDD004 represent four sequences that penetrate a large portion of the Upper Zone, and are correlated, based on marker horizons, resulting in approximately 500 meters stratigraphic thickness of the Upper Zone.

Analytical Procedures
All samples were crushed by hydraulic press and subsequently milled by an agate planet mill. Sample descriptions, and analytical details of both major and trace element data analyses, and whole-rock Lu-Hf isotope analyses can be found in previous publications (Nebel et al., 2013a,b).
For Fe isotope analyses, ca. 25 mg of whole rock powder was dissolved in a mixture of distilled, concentrated HF-HNO 3 in 15 ml Savillex vials at 130 degrees Celsius for 48 hours on a hotplate. After evaporation to dryness, each sample was redissolved in concentrated nitric acid and H 2 O 2 , and dried again until a clear solution was obtained. Finally, each sample was taken up in 9 M HCl, and Fe was purified using BIORAD AG-MP-1M chromatographic extraction resin. The details of the procedure and additional standard reference materials are given in Cheng et al. (2014).
Mass spectrometric analyses follow the protocols detailed in Sossi et al. (2015). In brief, ca. 6 ppm of sample material was spiked with 20 ppm Ni ICP standard solution for external mass bias control. All samples were analyzed three times, bracketed with the IRMM-014 standard reference material. Deviations in the per mil notation (δ 57 Fe) are reported as an average of repeated analyses. The external precision of the method is ±0.03 (2 s.d. of the mean) deduced from repeated analyzed of international reference materials as reported in Sossi et al. (2015), and is applied to all samples unless internal errors are larger. The internal error of each analyses is the standard deviation of the mean of the average of the three repeated analyses. The external error is the two standard deviation of the reference material. As a quality control measure the BCR-2 standard reference material was analyzed five times with a mean of δ 57 Fe = 0.14 ± 0.02 (identical to that reported in Sossi et al., 2015) with the error being the two standard deviation of the mean. Blanks were below detection limit and negligible. Table 1) vary from δ 57 Fe = −0.09 to +0.21 , and thus are spread over igneous rock compositions reported for other mafic systems (Williams et al., 2005;Schuessler et al., 2007;Weyer and Ionov, 2007;Poitrasson et al., 2013;Teng et al., 2013;Foden et al., 2015). A similar spread has also been reported for other layered intrusions, although some intrusions show a larger, whereas others smaller variations. There is no observed correlation between Fe isotope compositions and major element abundances (Ca-Al-Fe-Mg-Si-Ti), rock type or associated modal mineral abundance. Iron isotopes, however, show systematic variations with intrusion depth (Figure 3), with progressively lighter isotopes measured toward the base of the Upper Zone. These systematics are coherent and show a positive relationship with initial Hf isotopes (Nebel et al., 2013a), except for the magnetitites (Figure 3) at the base of the intrusion, which show a reverse trend between Fe and Hf isotopes.

DISCUSSION
Hafnium isotopes do not fractionate during igneous differentiation, thus any systematic isotope variations within the intrusion must result from two-component mixing of magmas with distinct Hf isotope values. Because we find that the whole-rock Fe isotopes largely follow a similar trend to the Hf isotopes, we conclude these variations must also relate to mixing of end-members with different Fe isotope values. The alternative is that the Fe isotope variation results from progressive crystal fractionation, which may have occurred contemporaneously with

Crystal Fractionation
Variations in Fe isotopes may reflect progressive crystal fractionation through removal of Fe 2+ -bearing phases (Weyer, 2008) within the Upper Zone magma batch during cooling. The trend observed in Figure 3 requires that minerals with an isotopically heavier composition accumulated at the roof of the intrusion, or that crystallization started in the centre of the Upper Zone and progressed toward the top and bottom with residual, isotopically heavier melts. In a closed magmatic system this indicates a bottom-to-top crystallization because light isotopes are incorporated first with the removal of Fe 2+ (and thus isotopically light iron) into olivine, notably with the exclusion of magnetitites. However, Hf isotopes do not fractionate during this process, and the accompanying trend in these radiogenic isotopes requires that during this processes of crystal settling, additional processes must have occurred that can explain the variation in radiogenic Hf isotopes. The only viable scenario to create a systematic Hf isotope gradient is involves a replenishing melt with a distinct Hf isotope composition, which would have progressively filled either the roof or the basement sections of an Upper Zone crystal mush area. In this scenario, the combined Hf-Fe isotopes demand that crystal fractionation (monitored through Fe) needs to progress contemporaneous to, and in concert with, melt injection (monitored through Hf isotopes). However, such a scenario, even though not impossible, seems thoroughly implausible. Effective and progressive crystal separation from liquid in a magma mush zone, as would be required to explain the combined trend in Fe-Hf isotopes through igneous differentiation alone, also ignores the repetitive modal layering of the intrusion. Double-diffusive convection (Huppert and Sparks, 1984) or small perturbation density horizons (Bons et al., 2015) would allow layers to develop individual Fe isotope batches within repetitive layering horizons. Yet this would need to occur progressively at the same rate that magma is being replenished into the system. The latter is required to satisfy the gradient in Hf isotopes and does not even account for potential differences in Fe isotopes between replenished and existing melts. We thus exclude igneous differentiation as the cause for the observed Fe-Hf isotope systematics.

Two-Component Mixing
Magma mixing between two components with distinct Hf isotope signatures requires one end-member to have a highly radiogenic Hf isotope signature at 2.8 Ga (i.e., much higher than the depleted mantle, ε Hf (t) ∼ +10 in whole rocks and >+15 in mineral separates), and a second component with a near primitive-mantle like composition with ε Hf (t)∼ +1 (Nebel et al., 2013a). For a mixing scenario the co-variations between Fe and Hf isotopes in Figure 4 for the silicate-rich horizons only (i.e., excluding magnetitites, see section "Discussion") also implies that the primitive mantle ε Hf (t) component is also primitive mantle-like in its Fe isotope composition with δ 57 Fe ∼0 (Figure 4).
For a layered intrusion whose magma is supposedly sourced from a mantle plume, it is plausible to assume that parental melts directly derived from within the plume are of basaltic to FIGURE 2 | Harker diagrams for the samples (Nebel et al., 2013b) indicating the modal mineralogy of the cumulate layers. Grouping of samples in these plots is defined by the dominant mineralogy with one group being dominated by olivine and pyroxene, one by feldspar and one by Fe-Ti oxides.
komatiitic composition (Campbell et al., 1989). Archean mantlederived liquids from komatiites show variations in their Fe isotope composition from δ 57 Fe −0.2 to +0.2 (Dauphas et al., 2009;Hibbert et al., 2012;Nebel et al., 2014). This range reflects to a large extent mineral accumulation and crystal fractionation during igneous differentiation of isotopically light olivine. The large degree of melting within a mantle plume (e.g., Sossi et al., 2016a) prevents Fe isotope fractionation during partial melting, such that melt isotope compositions directly relate to their source Fe isotope composition . One component in the Windimurra Upper Zone, being directly sourced from the mantle, can thus be constrained to a primitive mantle component (Sossi et al., 2016b).
The observed range in Fe isotopes requires that the second mixing component must have a heavy Fe isotope signature of ca +0.2 . Because tholeiitic layered intrusions fractionate to high Fe contents (>25 wt. FeO in whole rocks Vantongeren et al., 2010), and considering massive olivine-pyroxene fractionation with the preferential removal of isotopically light Fe 2+ during the formation of an underlying Middle Zone, residual liquids from this differentiation are expected to develop toward high Fe content with a relatively heavy Fe isotope composition. In addition, progressive crystallization of the Middle Zone will concentrate these residual, isotopically heavy liquids toward the top of the intrusion through compaction; compression of the cumulate mass is likely to trigger filter pressing (McBirney, 1995;Mathez et al., 1997). There, it forms interstitial liquid or a melt pool within a crystal mush. This liquid is defined here as liquidin-residence with a δ 57 Fe ∼ +0.2 . The origin of the extreme Hf isotope signature of this melt is uncertain, but is not unique and may be related to similar sources of extreme Hf isotopes observed in komatiites elsewhere (cf. Nebel et al., 2014).
Pulses of hot and buoyant replenishing melt from the mantle plume source will preferentially intrude into these existing crystal mush zones in an evolving Middle Zone magma chamber. Replenishing melt can thus be expected to fill areas where the crystal/melt ratio is at its lowest, thereby creating a new, larger melt-pool in the crystal mush. Mixing of liquids-in-residence was probably filter-pressed and pooled from underlying sequences in an evolving Middle Zone (heavy in Fe isotopes with a radiogenic Hf isotope signature) in a homogeneous Hf-Fe isotope melt pond, yet with various antecrysts and possibly also partly layered. Given the sheer volume of melt required to form the Upper Zone, represented here by ca. 500 m of drill core intersection, the newly formed melt pool would likely be dominated or at least strongly affected by the fresh melt batch. Inevitably, however, crystals from the previous mush zone would be present in this blend of melts. It remains unclear as to what extent a melt injection can disturb existing layering, but the perturbation within the melt pool is expected to be severe with respect to temperature and viscosity difference between new melt and liquid-in-residence. If magma cools below its liquidus temperature, rheologic properties change significantly in an interlocking crystal network . However, reheating through magmatic recharge can cause remobilization of crystals (Irvine et al., 1998) and magma mixing occurs while a melt network is still present , especially if both liquids are compositionally similar. Replenishing liquid will likely be much hotter and with lower viscosity, and it is assumed that within the Middle Zone magma ponds, crystallinity was still <25 , which is plausible for a spongy crystal mush at the timescales of LMI formation (Bergantz et al., 2015). We thus propose that the replenishment forged a new Upper Zone in the intrusion through a massive inflation in volume consequent to a melt surge.

Unmasking the Mixing of Similar Liquids
Despite the apparent difference in isotope compositions, both the replenished and residential liquids will be (although not isotopically but) compositionally similar, yet not identical. It is expected that the residual liquids have similar Si but higher Fe concentrations following a tholeiitic fractionation trend, whereas newly injected melts are higher in MgO with a possible komatiitic composition. The new melt blend must thus be in disequilibrium with existing crystals that formed from earlier liquids, causing some remelting. Overgrowth of crystals (and zonation) is likely common during subsequent cooling. This would lead to moderate chemical and strong isotope heterogeneity in individual minerals, a feature observed for Fe isotopes in some Liu et al., 2014), but not all (Bilenker et al., 2017) layered intrusions. In the mixing scenario envisaged here, Fe isotope systematics of mineral pairs or within minerals are thus not useful to decipher the history of LMI petrogenesis, or would require high resolution isotopic analyses that is currently not achievable. Strong and non-systematic variability is expected FIGURE 4 | Plot of stable Fe vs. age-corrected radiogenic Hf isotopes. Replenishing liquid (red star) with δ 57 Fe ∼ 0 and ε Hf (t)∼ + 1 is admixed to liquid-in-residence, forming a new liquid (orange star). The Fe isotopes are a mixture between both liquids, whereas Hf is initially controlled by the liquid in residence. The final whole rock compositions of the samples, however, represent a mixture between minerals that formed from the liquid-in-residence (purple star) or stage-1 crystals that formed from it and newly formed minerals from the mixed liquid (orange star). Dependent on respective proportions, the whole rock will plot along the mixing trend. Magnetitites have probably retained a predominant memory of stage-1, Middle Zone liquid but will probably also record crystal chemistry affects in Fe isotopes.
if melts are very different in Fe isotope compositions, as observed in the Beiba intrusion Liu et al., 2014). If they are identical, no variation is observed at all, as is the case for the Upper Zone of the Bushveld Complex (Bilenker et al., 2017).
The Windimurra Upper Zone, however, shows considerable heterogeneity in Hf isotopes that are systematic with depth and Lu/Hf Nebel et al. (2013a). Because blends of melts are expected to have a homogeneous isotope signature, the only way this can occur is through isolation of melt batches or through freezing of Hf isotopes in antecrysts from a previous melt batch. Isolating melts is difficult to reconcile with a turbulent melt injection, so we consider the latter scenario more plausible. Indeed, Nebel et al. (2013a) reported extreme Hf isotope signatures in pyroxene separates with ε Hf [t] = + 15 to + 20, which is in strong support for crystals that grew from an isotopically different melt batch. It is argued here that this melt batch is the former Middle Zone.
In the proposed scenario of melt injection into an existing crystal mush zone, the systematic increase in initial Hf isotopes toward the roof (Nebel et al., 2013a) requires that a substantial amount of crystal cargo in the Upper Zone is inherited from crystallization of Middle Zone melts with distinct isotope signatures, and that these are located at higher stratigraphic levels. This initial Middle Zone crystallization is referred to here as stage 1 crystallization and is illustrated in Figure 5. Stage 2 crystallization in the newly established Upper Zone will then follow with overgrowth on existing grains and growth of new minerals with a different Hf-Fe isotope composition. The ratio of stage 1/stage 2 minerals then defines the whole-rock isotope mixing trend observed in Figure 4 (cf. Nebel et al., 2013a). Scatter is expected in any trend resulting out of this mix, as interstitial liquid will have Hf isotope signatures of the new melt mix, whereas Fe isotopes can be further subject to fractionation through crystal formation, two independent factors captured in a single analysis.

Isotope Contraints on Magma Mush Dynamics
Adopting the mixing scenario outlined above for Fe isotopes implies that, within the silicate layers, each whole rock Hf-Fe isotope composition is a rough reflection of the stage 2/stage 1 crystal ratio (or replenished melt blend vs. inherited crystals), which appears to gradually decrease toward the roof of the intrusion. Crystal settling in the chamber must have occurred through convection and a density contrast within the magma chamber, with the majority of the crystal cargo accumulated at the roof of the intrusion. It remains unclear why melt with the lowest crystal fraction is concentrated toward the base of the intrusion, but possibly because of an injections into an already existing, dense melt pond. Notable is that the lower density of komatiite liquid compared to existing crystals requires a physical process stopping the liquid to progress further towards the roof. A second effect can be imperfect and local limited perturbation of existing crystals during replenishment toward the roof, such that the chamber is filled and expanded (like a balloon) through an existing melt pond at the base of the newly established Upper Zone. The roof section is less affected by the melt injections. Melt infiltration, however, must have been sufficient to infiltrate all layers, albeit to different degrees, and most likely followed by rapid establishment of rhythmic layering in individual layers (Bons et al., 2015).
A remaining puzzling aspect of the intrusion are the basal Fe-Ti-rich layers of the Upper Zone, which show among the lowest ε Hf (t) of the sequence but a reversal in the Fe isotope signature. The Hf isotopes demand that these layers have formed from a blend of melts, the replenished liquid and the liquid-in-residence. At least some Fe-Ti oxides must thus postdate the replenishment, and were only then accumulated close to the base of the newly formed Upper Zone. For this, it seems plausible that the replenished liquid, albeit lower in Fe content, was more oxidized than the liquid-in-residence, triggering magnetite saturation through elevated Fe 3+ / Fe levels. The lowermost magnetite horizons show the highest δ 57 Fe with incorporation of isotopically heavy Fe 3+ into these early during crystallisation. These either form instantly and with this crystal chemical control on Fe isotopes overrides the Fe isotope signatures of the liquids (incorporating preferentially heavy isotopes, Shahar et al., 2008). Alternatively, a redoxreaction forced instant magnetite saturation of the liquidin-residence with later magnetite formation being dominated by Fe of replenished liquid, and with this a lighter Fe isotope composition. An aspect that needs to be accounted FIGURE 5 | Schematic sketch of the temporal two-stage evolution of the Upper Zone magma-reservoir. Prior to melt replenishment, a melt pool, possibly with initial layering, was established in the roof zone of the (now) Middle Zone of the intrusion, which formed the top part of the intrusion. Residual liquids accumulated through filter pressing with a heavy Fe isotope composition of δ 57 Fe ∼ + 0.2 and extreme radiogenic Hf isotopes (purple star in Figure 4). During melt replenishment, a perturbation in the magmatic pool allowed melt and crystal re-organization. During this second stage, the new Upper Zone formed, which constitutes a mixture of newly injected melt and liquid-in-residence and its crystal cargo. Crystals that form during this process would be in disequilibrium with any hybrid melt pool. A single parental liquid did not exist.
for in our study is, however, that we report whole rock analyses. The Hf isotope composition in these rocks, composed of almost entirely magnetite, can then be dominated by interstitial liquid, which is compatible with regular cumulate behavior, whilst the Fe is driven by magnetite. The Hf-Fe isotope co-variation in the magnetitites is then only an artifact of two independent factors, combined in a whole rock analysis.

CONCLUSION
The scenario promoted here can explain a range of geochemical oddities in the Upper Zone of the Windimurra Igneous Complex and possibly other layered intrusions. Cryptic melt replenishment with incomplete mixing and variable melt-crystal ratios may be responsible for: (1) apparent chemical disequilibrium among crystals, as to overgrowth, remelting and formation of a second generation of crystals from a different magma batch; and (2), cryptic isotope layering among phases, which relate to magma replenishment of melts that appear similar in major elements but are distinct in their trace isotope systematics. Notable is that this is the case for major and trace elements. This sequence of events can explain reported Fe isotope disequilibrium among some mineral separates Liu et al., 2014). For a cyclic layering to develop, however, smaller zones must have subsequently developed into density horizons with possible internal small-scale convection. This process of layering remains enigmatic and is unrelated to the replenishment, but likely continues as proposed for other layered intrusions (e.g., Bons et al., 2015).
Petrologic investigations of individual mineral phases may, in general, fail in tracking the replenishment process, because of the mix of crystals and melts in disequilibrium, perturbation of crystals, coupled with subsequent convective reorganizations of different generations of crystals. This will likely result in large variations both chemically and isotopically on a mineral scale that, in addition to subsequent layering, create a pseudosystematic agglomeration of minerals with only some chemical memory of their turbulent past. Whilst mineral studies and whole rock analyses are often complementary, individual minerals would have been affected by various stages of the mixing process and to different degrees, so may only show a snap shot of what is a complex mechanism of recharge, mixing and physical sorting. Under these circumstances, and counterintuitive to most recent studies, whole rock composition possibly yield a more detailed picture of the history of the intrusion. Noteworthy is, though, that this requires a special sequence of intrusive events and mantle sources, as is the case for Windimurra. In the context of melt replenishment into magma chambers of any kind, a twocomponent mixing as is observed here may well, at least in parts, explain the widely observed Hf isotope heterogeneity (in zircons) in plutonic bodies.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary material.

AUTHOR CONTRIBUTIONS
ON and RA designed the study and the model. ON, PS, and NG wrote the manuscript. TI, AB, and RL contributed to the sampling and the data processing. All authors commented on the manuscript.

FUNDING
This study was supported by a ARC linkage grant to RA and an ARC Future Fellowship to ON (FT140101062).

ACKNOWLEDGMENTS
The Melbourne TIE team is acknowledged for support. TI publishes with permission from the Executive Director of the Geological Survey of Western Australia. Maximus Resources is kindly acknowledged for providing core material and significant contribution to an associated ARC linkage project. Reviewers are kindly acknowledged for their comments and MA for his editorial handling.