Multicomponent diffusion in aluminosilicate garnet: coupling effects due to charge compensation

ABSTRACT Three kinds of multicomponent diffusion effects, arising from three distinct physical mechanisms, are evident in stranded diffusion profiles at the rims of partially resorbed garnets from the contact aureole of the Makhavinekh Lake Pluton, northern Labrador. Profiles that display a subtle maximum in Ca concentration are explained by the prevailing ideal mean-field theory of multicomponent diffusion, but models implementing that theory cannot replicate inverted profiles for Li and internal maxima for Nd, Sm, and Eu. The anomalous profiles are quantitatively reproduced, however, by numerical simulations employing a model based on coupled movement of charge-compensating groups during diffusional transport of yttrium and the rare-earth elements (Y+REEs). An alkali-type charge-compensation mechanism for the heterovalent substitution of Y+REEs on dodecahedral sites in garnet produces direct charge coupling between Li+ and (Y+REE)3+ and leads to co-diffusion of Li+-(Y+REE)3+ pairs, with the result that Li profiles closely mimic those for Y+REEs. A menzerite-type charge-compensation mechanism produces indirect charge coupling among all Y+REE components, with the result that the fluxes of low-abundance REEs become partly dependent on the fluxes of Y+REEs present in higher abundance. These findings have implications for the robustness of Li profiles in garnet as monitors of fluid–rock interaction, for geochronology based on the Sm–Nd and Lu–Hf systems, and for future experimental attempts to quantify rates of diffusion in garnet.


Introduction
Multicomponent diffusionin which the flux of a component is driven not only by its own concentration gradient, but also by factors related to other componentshas conventionally been approached in silicate minerals by application of an ideal mean-field theory (e.g. Lasaga 1979). In that approach, the flux of each component in a solid solution is linked, through a combination of electromagnetic and thermodynamic forces, to the gradients in concentration of all components, and the simplifying assumption is made that the diffusing components comprise a solution that can be regarded as thermodynamically ideal.
In many instances, the ideal mean-field theory correctly captures diffusional behaviour in silicates. But in other cases, multicomponent effects appear that are not encompassed by such a theory. A salient example is the diffusion in garnet of monovalent and trivalent ions, substituting for divalent ions, on the sublattice of dodecahedral sites: observations in natural systems reveal diffusional coupling among components that can be traced to the requirement of local charge compensation for ions involved in heterovalent atomic substitutions.
This article reviews evidence for diffusional coupling arising both directly and indirectly from local charge compensation, then derives a new empirical model for indirect coupling that quantitatively replicates previously unexplained observations of some multicomponent diffusion effects in garnet. light when dealing with components involved in heterovalent substitutions in minerals.

Ideal mean-field theory
The specific method most frequently employed for evaluation of multicomponent diffusion in silicates is the one promulgated by Lasaga (1979), which falls into the category of mean-field models. Multicomponent effects are treated by considering both electromagnetic forces arising from the movement of diffusing ions and thermodynamic forces arising from gradients in chemical potential. The electromagnetic forces are evaluated by neglecting forces arising from the instantaneous correlation of the motion of one ion and the motion of a nearby ion. This simplification is justified by the fact that the neglected forces act only over length scales much shorter than those of the forces arising from Coulomb interactions in an ionic crystal. In consequence, the electromagnetic force on an ion can be treated simply as the average force resulting from superposition of the electrostatic fields produced by all ions in the crystala 'mean-field' description of the electrical potential. The thermodynamic forces are the gradients in chemical potential, which are linked to concentration gradients through activity coefficients.
Using this mean-field approach, the total diffusional flux of a component in a multicomponent solution is described by combining the effects of both electromagnetic and thermodynamic contributions. Such a description leads directly to an expression for the multicomponent diffusion coefficients D ij that relate the flux J i of each component i to the concentration gradient ∇c j for each component j in the non-equilibrium thermodynamic relation in which n is the number of components in the solution.
The D ij are commonly grouped into a matrix so that the full set of diffusional fluxes, gathered into a vector J, can be obtained as the product of the D matrix and a vector of concentration gradients ∇c, via: A lengthy expression for the multicomponent diffusion coefficients D ij in a mean-field model appears as Equation (14b) of Lasaga (1979); that equation is applicable to non-ideal solutions. In common practice, however, the complexity of such expressions is markedly reduced by making the assumption thatat least for the purpose of assessing diffusional behaviourthe solution is thermodynamically ideal; this simplification also obviates the need for determining the activity coefficients for the solution, which may be unknown quantities. For ideal solutions, it is possible to express the multicomponent diffusion coefficients in a simple form, as for example in Equation (15) of Lasaga (1979): in which D 0 i is the self-diffusion coefficient for component i; δ ij is the Kronecker delta (δ ij = 0 if i ≠ j, or δ ij = 1 if i = j); z i is the charge on ion i; and c i is the concentration of component i in the solution. Equation (3) is the defining expression for the ideal mean-field theory of multicomponent diffusion, the approach most commonly employed to treat diffusion in silicates.

Application of ideal mean-field theory to diffusion in garnet
The ideal mean-field approach to multicomponent diffusion is widely used in many mineral systems, and it has been essential to understanding diffusion in garnet, which is a solution of (at minimum) four significant chemical components. The ideal mean-field theory has long found application as a means of exploiting experimental and empirical data to constrain rates of diffusion of divalent cations on the dodecahedral sublattice of garnet; prominent examples are Chakraborty and Ganguly (1992), Ganguly et al. (1998a), Carlson (2006), and Perchuk et al. (2009). Geological applications of this approach, used as a means of interpreting diffusional effects on chemical zoning in natural garnet, are exceedingly numerous; many examples may be found among the selection of studies highlighted in the review by Ganguly (2010).
Although application of the ideal mean-field model is the norm for studies of diffusion in garnet, two noteworthy exceptions deserve mention. In an experimental study of ternary Fe-Mg-Ca interdiffusion, Vielzeuf and Saúl (2011) produced profiles with a variety of features characteristic of multicomponent diffusion effects, and retrieved by inversion methods the four elements of the D matrix relevant to the specific compositions and physical conditions of each experiment; thus in each instance they acquired values for the two off-diagonal coefficients directly, obviating the need for any fundamental theory for the origin of the multicomponent interactions. Borinski et al. (2012) broke new ground by employing the full non-ideal mean-field equations of Lasaga (1979), incorporating previously determined activity coefficients for the quaternary Mg-Fe-Mn-Ca system. This permitted exploration of the range of compositions over which non-ideality plays a significant role in the diffusion of these divalent cations.

Effects of incorporation of monovalent and trivalent ions into dodecahedral sites
The applications described above are restricted to instances in which the interacting components are all divalent ions, the predominant occupants of dodecahedral sites in garnet. But diffusion of monovalent and trivalent ions on those sites introduces a vital new factor. In such cases, one must consider the need for local compensation of electrostatic charge; such compensation requires pairing of cations with reciprocal excesses and deficiencies of charge, and maintenance of mutual proximity for ions in those pairs during their movement.
The identity and nature of these charge-compensating pairs in garnet is revealed by recent lattice-dynamics calculations that examined the relative energetic costs of different means of incorporation of monovalent and trivalent ions into aluminosilicate garnet . Two substitution mechanisms were found to be substantially lower in energetic cost than all known alternatives, and thus are regarded as significantly more probable in nature: these are the alkali-type and menzerite-type substitutions.
Both substitutions incorporate large trivalent ions like yttrium and the rare-earth elements (Y+REEs) into the dodecahedral sites of the garnet lattice, sites with eightfold coordination by oxygen. Because these sites normally house divalent ions, balancing the electrostatic charge requires compensating substitutions on adjacent sites. An alkali-type substitution mechanism ( VIII LI or VIII NA) provides charge balance for a (Y+REE) 3+ ion on a dodecahedral site by introducing Li + or Na + into a neighbouring dodecahedral site. A menzerite-type substitution mechanism (MNZ-MG or MNZ-FE) instead introduces Mg 2+ or Fe 2+ into a neighbouring octahedral site, in place of Al 3+ .
The evidence from natural occurrences reviewed in Carlson et al. (2014Carlson et al. ( , p. 1030Carlson et al. ( -1031 corroborates the prediction from lattice dynamics that the alkali-type and menzerite-type substitution mechanisms are significantly more prevalent in nature than any alternatives, specifically the IV AL (YAG), VI LI, or VIII  (vacancy) substitutions. In particular, for rare cases in which Y+REE contents of natural garnet achieve levels high enough for sufficiently accurate analytical determination of their charge-compensation mechanisms (e.g. Enami et al. 1995;Grew et al. 2010), the alkali-type and menzeritetype schemes are found to predominate strongly.
The linkage of ions into charge-compensating pairs is problematic when contrasted with the assumptions underlying ideal mean-field theory. Ion-pairing requires strong short-range interactions in which the motion of one ion in the pair directly influences the motion of the other, and those effects cannot be merely transient (as for unpaired ions), because they are required for the ion pair to remain associated as it moves through the crystal. Analogously, Lasaga (1979, p. 458) alludes to the inadequacy of the ideal mean-field theory for treatment of diffusion in aqueous solution of components like Na 2 SO 4 , whose incomplete dissociation generates appreciable concentrations of dissolved ion pairs. In crystalline solutions, moreover, the interactions of substituent monovalent and trivalent ions with divalent ions in the host are more likely to depart significantly from thermodynamic ideality than in the case of homovalent substitution, owing to the differences in charge and potentially large differences in size from ions in the host.
Because these key features of charge-coupled diffusion of monovalent and trivalent ions in garnet are not encompassed by the ideal mean-field theory, one might expect those ions to display unanticipated diffusional effects in natural multicomponent systems. Anomalous behaviourdiffusion that defies explanation in terms of the ideal mean-field modelhas in fact been observed for both Li and some REEs (Nd, Sm, and Eu), in a wellstudied example of partially resorbed garnet in an ordinary metapelite, as described in the following.

Previous work: multicomponent diffusion effects in garnet from the MLP aureole
The contact aureole of the Makhavinekh Lake Pluton (MLP), part of the Nain Plutonic Suite of northern Labrador, Canada, provides an extraordinary natural laboratory for the study of multicomponent diffusion in aluminosilicate garnet. Garnet crystals in metapelites initially achieved a high degree of compositional homogeneity during prolonged granulite-facies regional metamorphism; then, after substantial exhumation, they were partially resorbed during a contact-metamorphic episode resulting from intrusion of the pluton. In the course of resorption, differential partitioning of both major and trace elements between garnet and its reaction products drove diffusion that produced dramatic zoning in the outer rims of the relict garnet crystals. Analysis of the stranded diffusion profiles that comprise this zoning, together with numerical simulation of their evolution, generates novel insight into the fundamental nature and mechanisms of multicomponent diffusion in natural aluminosilicate garnet.

Evolution of stranded diffusion profiles in partially resorbed garnet
The Makhavinekh Lake Pluton and its surrounding contact aureole have been analysed extensively in a series of field, mineralogic, petrographic, thermobarometric, geochronologic, and thermal-modelling studies (e.g. Berg 1977aBerg , 1977bLee 1987;Ryan 1991;Bertrand et al. 1993;Van Kranendonk et al. 1993;McFarlane et al. 2003McFarlane et al. , 2005aMcFarlane et al. , 2005bPattison et al. 2003;McFarlane and Harrison 2006;Carlson 2010;Kelly et al. 2011). As a consequence, the evolution of stranded diffusion profiles (SDPs) in partially resorbed garnet crystals from the aureole is well understood, and the conditions and processes that produced the SDPs are well constrained. A brief synopsis of key points from the extensive literature cited above is provided here.
The Tasiuyak Gneisshost to the MLP intrusionis a geographically extensive, lithologically and geochemically homogeneous, layered migmatitic paragneiss that underwent granulite-facies metamorphism (~850°C at 0.6-0.9 GPa) and pervasive partial melting at 1860-1850 Ma. Exhumation to shallower depths occurred at ca. 1794-1786 Ma, so granulite-facies conditions persisted for 60-70 million years. This interval was sufficient to produce complete homogenization of major divalent cations (Mg, Fe, Mn, Ca) in garnet, and to very nearly homogenize concentrations of slower-diffusing cations (e.g. Li, Y+REEs, and Cr), leaving concentration profiles that slope extremely gently from core to rim.
At 1322 Ma, a 1000 km 2 composite pluton was intruded, cored by cumulate troctolite, anorthosite, and ferrodiorite, and sheathed in charnockitic quartz monzonite grading outward to granite. The intrusion produced a 5-6-km-wide, virtually anhydrous aureole during an isobaric contact-metamorphic episode at~0.53 GPa that reached peak temperatures ranging from~700°C tõ 900°C. Garnet crystals within the aureole were partially to completely resorbed during the contact-metamorphic event, in reactions that produced coronal symplectites of cordierite + orthopyroxene ± plagioclase ± K-feldspar. Outside the aureole, at distances greater than 5-6 km from the intrusion, unmodified garnet persists; within the aureole, the extent of garnet resorption increases with proximity to the contact, and little or no garnet survives within~450 m of the pluton.
The strongly anhydrous nature of both the anorthositic intrusion and its high-grade granulite host rules out appreciable convective transport of heat during the contact metamorphism, allowing confident use of conductive thermal models to determine T-t histories as a function of distance from the contact with the approximately cylindrical pluton. These models are tightly constrained by Al-in-orthopyroxene thermometry across the coronal reaction textures, which provides narrow limits on peak temperatures during the resorption event and on temperatures at the cessation of garnet resorption, both of which vary, of course, with distance from the intrusive contact.
During garnet resorption, the rim of a garnet crystal continually changes composition, tending towards local chemical equilibrium with product minerals in the surrounding reaction zone. Some elements partition preferentially into the relict garnet, building up high rim concentrations that drive inward diffusion; this effect is seen for Mn and Fe, and likewise for Y, Cr, Gd, Tb, Dy, Ho, Er, Tm, Yb, and Lu. Other elements partition preferentially into the product minerals, yielding depleted concentrations in relict garnet rims that drive outward diffusion; this effect is seen for Mg and Ca, and likewise for Nd, Sm, and Eu. The radial extent (width) and shape of the diffusionally generated concentration gradients is determined by the relative rates of resorption reaction and inward communication by diffusion of changes in rim composition. Thus the character of the concentration profiles varies dynamically over the course of the resorption reaction, reflecting the changing kinetics of reaction and diffusion during heating and cooling of the aureole. The resorption reaction ceased at temperatures high enough to permit continued diffusion within the relict garnet crystal during cooling, which introduces inflection points into the concentration profiles while producing considerable flattening of the gradients present at the end of resorption. An example of the modelled time evolution of a concentration profile for Y appears as Fig. 2 of Carlson (2012).
Samples selected from a suite collected previously along traverses outward from the intrusive contact ) yielded a diverse array of variably resorbed garnet crystals for analysis. In prior work (Carlson 2006(Carlson , 2012Cahalan et al. 2014), relict crystals were sectioned precisely through their morphological centres as located in 3D by high-resolution X-ray computed tomography (HRXCT), and SDPs for major elements and a variety of trace elements were measured by one or more techniques, depending on elemental abundance: techniques employed were electron probe microanalysis (EPMA), laser-ablation inductively coupled plasma mass spectrometry (LA-ICPMS), and secondary-ion mass spectrometry (SIMS). Details of analytical procedures and data-reduction methods can be found in the literature cited above.

Observed effects of multicomponent diffusion
Stranded diffusion profiles measured in prior studies evince three separate modes of multicomponent interaction in partially resorbed garnet in the MLP aureole. The observational evidence for each such mode of interaction is reviewed here, and for each, a distinct atomic-scale mechanism is proposed to explain the salient features of the SDPs. In the following section, each proposed mechanism will be tested against models for evolution of the measured SDPs.

Uphill diffusion of Ca
A subtle maximum in Ca concentration is a common, though not universal, feature of profiles measured on garnet crystals from the outer portion of the MLP aureole (Carlson 2006); an example is shown in Figure 1. The development of a maximum in an SDP, starting from an initially homogeneous composition, requires transport of material in directions of increasing concentration. This phenomenon, termed 'uphill diffusion', is an expected consequence of multicomponent interactions in the context of the ideal mean-field model, as discussed in detail by Vielzeuf and Saúl (2011). It results when one (or more) of the cross-coefficients in a matrix of multicomponent diffusion coefficientsthese are the off-diagonal terms that link the flux of one component to the concentration gradients of other componentsis negative, and thus opposite in sign to the self-diffusion coefficient. If that cross-coefficient and its corresponding concentration gradient are large enough, the resultant flux (the mathematical product of those two factors) can overbalance and thus counteract the flux of the component down its own concentration gradient.
If this mechanism is responsible for the origin of the Ca maxima observed in garnets from the MLP aureole, then a model based on the ideal mean-field theory should be capable of replicating it.

Reversed partitioning and inward diffusion of Li
Measurement of Li profiles in garnets from the inner portion of the MLP aureole (Cahalan et al. 2014) produced an unforeseen and initially puzzling result. Partition coefficients measured by Dutrow et al. (1986) and by Dutrow et al. (2011) indicate that for garnet with compositions typical of those found in pelitic assemblages, equilibrium partitioning of Li strongly favours cordierite over garnet (K D = 200). Thus the expectation before measurement was that in garnets of the MLP aureole, which are surrounded by coronal reaction zones dominated by cordierite, Li concentrations would fall off steeply towards the rims of the relict crystals. But SDPs for Li were instead found to be inverted, in the sense that they have highly elevated concentrations within the relict rims ( Figure 2); thus Li behaved like Mn or Y, elements that preferentially partition strongly into garnet rather than into its reaction products. In fact, Li profiles exhibit close similarities to Y and Yb profiles, and all three elements share precise congruence in the location of the steepest portion of  their SDPs, marked by the inflection point in each profile's curvature, which is a sensitive measure of each element's inward diffusive penetration. The inversion of Li profiles and their close correspondence to profiles for Y and Yb both suggest that Li transport is linked directly to the transport of Y+REEs. As described above, the VIII LI substitution mechanism is a low-energy means of providing charge compensation for some of the Y that partitions into the relict garnet during resorption. Consequently, enrichment of Y+REE in the relict rim stabilizes Li there as well, reversing the tendency towards its outward partitioning that is observed for garnet of low Y+REE content. Because local charge compensation must be maintained as Y +REEs diffuse inward down their concentration gradients, some of the Y+REEsand most if not all of the Li 1 must be moving as Li + -(Y+REE) 3+ ion pairs. If this mechanism is responsible for the inverted Li profiles, then it should be possible to replicate them using a model implementing the ideal mean-field theory, in which the self-diffusivity of Li matches that of Y+REEs.
Uphill diffusion of Nd, Sm, and Eu As described in detail in Carlson (2012Carlson ( , p. 1605Carlson ( -1607, a universal feature of the SDPs for Nd, Sm, and Eu measured by SIMS in relict garnets of the MLP aureole is the presence of a small but distinct maximum near the inward limit of penetration of diffusional effects for these and other trivalent species (e.g. Figure 3), contrasting with a steep drop in concentration from the maximum outward. (Hereafter, these three elements with similar behaviour are grouped together and for convenience are designated as L-REEs = 'lighter rare-earth elements'. This designation emphasizes that these elements the only REEs that partition preferentially out of garnet into the product assemblage during the MLP resorption reactionshave smaller atomic numbers, lower gramformula weights, and larger ionic radii than the other REEs studied, while recognizing that Eu is commonly considered to be a 'middle' rare-earth element rather than a 'light' rare-earth element sensu stricto.) The maxima for the L-REEs are systematically positioned near locations where the concentrations of Y and the heavier REEs begin to rise appreciably, suggesting that the maxima might be related in some way to multicomponent diffusion effects. But their origin was not definitively determined in prior work, partly because the ideal mean-field model used in Carlson (2012) to extract rates of diffusion for Y+REEs was limited to the four principal divalent components (Mg, Fe, Mn, Ca) and one additional component, either Y or a single one of the REEs. This meant that while the model could reproduce ideal mean-field effects involving the four divalent cations and any one of the three L-REEs, multicomponent effects due to interactions between or among Y +REEs would not be captured. None of the five-component models succeeded in generating profiles that contained maxima as observed, 2 but the efficacy of an ideal mean-field model that accounts for interactions between Y+REEs (e.g. four divalent cations + Y + one L-REE) remains untested.
The hypothesis advanced in Carlson (2012Carlson ( , p. 1613Carlson ( -1614 was that the L-REE maxima were produced by multicomponent effects that conjoined the motion of the L-REEs with the motions of Y and all REEs heavier than Eu. If such a linkage exists, the L-REEs are subject to two competing processes during resorption: on one hand, their preferential partitioning into the product assemblage reduces their rim concentration so that diffusion draws them outward; on the other hand, to the extent that their motion is coupled to that of components that diffuse inward because they partition preferentially into the garnet, the L-REEs would be drawn towards the centre of the relict crystal, moving along with Y and all REEs heavier than Eu. The maxima for the L-REEs were conceived to be the net result of these two competing processes, in much the same way that maxima are generated under the ideal mean-field theory when the net flux is the sum of competing partial fluxes with opposite sign due to off-diagonal terms in the diffusion matrix. Lengths of horizontal bars are equal to width of incident ion beam. Initial rise at~300 µm corresponds to location of furthest inward diffusive penetration of Y+REEs, and contrasts with abrupt outward drop in concentration, which results from preferential partitioning of Sm into reaction products during resorption. Traverse M04A1f, after Carlson (2012); used by permission.
The mechanism for coupling among Y+REEs was envisioned to be an indirect linkage that followed from a menzerite-type charge-compensation mechanism. By assuming that the Y+REEs in the garnet structure were incorporated primarily via a menzerite-type substitution (in which Mg 2+ or Fe 2+ replaces octahedral Al 3+ ), it was possible to explain some otherwise puzzling observations, namely the weak dependence of Y+REE diffusivity on ionic radius and on host-garnet composition, and the close equivalence of their diffusivity with that of Cr 3+ , which diffuses on the sublattice of octahedral sites (Carlson 2012(Carlson , p. 1613. The presumed predominance of a menzerite-type incorporation mechanism also implied that the diffusion of all Y+REEs should be coupled indirectly, by virtue of the fact that all Y+REEs are locally paired with VI (Mg, Fe) to maintain local charge balance. As the concentration of any one of the Y+REEs increases, the local concentration of VI (Mg, Fe) increases equally, favouring concomitant incorporation of all other Y+REE ions by simple homovalent exchange on dodecahedral sites.
As noted above, lattice-dynamics calculations subsequently confirmed that menzerite-type substitutions are most likely the dominant means of charge compensation for Y+REEs present in excess over Li + and Na + , lending credence to this hypothesis of indirect diffusional coupling among all Y+REEs. If this indirect-coupling mechanism is responsible for generating the L-REE maxima, then one would expect failure of a six-component model (four divalent cations + Y + one L-REE) based on the ideal mean-field theory, because that theory does not account for such effects. Instead, success might come from a six-component model in which L-REE fluxes are connected to the flux of the much more abundant Y.

Models of multicomponent diffusion in garnet from the MLP aureole
The mechanisms proposed above to explain the variety of observed multicomponent interactions constitute hypotheses that can be tested for consistency with numerical models of the evolution of the SDPs. After describing in brief how such simulations are conducted, this section evaluates models based on the proposed mechanisms by comparing their predictions with observed profiles.

Numerical simulations
The numerical simulations presented here for evolution of the SDPs are identical to those employed in prior work (Carlson 2012(Carlson , p. 1601(Carlson -1602, to which the reader is referred for full details), except for two differences. First, they treat a six-component solution: Mg, Fe, Mn, Ca, Y, and either Li or one L-REE (Nd, Sm, or Eu). Second, they allow, as an option, a means of coupling the flux of an L-REE component to the flux of Y by an additional mechanism that supplements the ideal mean-field theory interaction.
The simulations calculate, for each measured SDP, the evolution of chemical zoning within a spherical crystal undergoing simultaneous resorption and multicomponent intracrystalline diffusion, while traversing a specified temperature-pressure-oxygen-fugacity-time (T-P-f O 2 -t) path. The development of each SDP is fully determined by specification of four inputs: (1) an initial concentration profile; (2) the self-diffusion coefficients for each of the diffusing components, accounting for their dependence on variations in temperature, pressure, and oxygen fugacity during the resorption/diffusion episode; (3) the T-P-f O 2 -t history for the sample containing the modelled garnet; and (4) the activation energy for intergranular diffusion, which defines the dependence on temperature of the rate of the resorption reaction. In the present work, initial concentration profiles were those derived from measurements in Carlson (2006), Carlson (2012), and Cahalan et al. (2014), and intracrystalline diffusivities under the ideal mean-field theory were those determined by Carlson (2006) and Carlson (2012). In keeping with prior efforts, T-t histories as well as the activation energy for intergranular diffusion were those presented in Carlson (2010); P was regarded as constant at 0.53 GPa; and f O 2 determined from measured oxide-sulphide-silicate phase equilibria was extrapolated in temperature assuming a constant differential from the quartz-fayalite-magnetite buffer.
In the absence of charge-coupling effects, the model includes for each element only a single parameter that is not known a priori, namely the parameter that fixes the fraction of the element that is retained within the outermost rim of the relict crystal vs. the fraction lost to the reaction products as resorption progresses. This retention/loss factor is uniquely determined, however, by the material-balance constraint provided by the measured total amount of each element in the relict crystal at the end of the resorption episode. It must be evaluated by iterative trial-and-error adjustment over the course of multiple modelling attempts; iterative attempts continued until the material-balance constraint was satisfied within a narrow tolerance, typically 0.05% relative.
After the simulation reaches the temperature at which resorption ceased (as determined by Al-in-orthopyroxene thermometry for crystals in contact with the relict garnet rim), the garnet volume remains fixed and exchange with the reaction products ends, but additional intracrystalline diffusion occurs within the crystal until cooling along its T-t path brings it to temperatures low enough for further diffusion to be negligible.

Ca: mean-field coupling with other divalent cations
Because Ca diffusion requires only homovalent exchanges on the dodecahedral sites in the garnet structure, it conforms closely to the assumptions of the ideal mean-field theory. Figure 4 shows that a four-component numerical model implementing that theory quantitatively replicates the observed maximum in the Ca zoning profiles. (It also simultaneously replicates the spatially coincident profiles for Mg, Fe, and Mn, as seen in Fig. 4 of Carlson 2006). This result confirms that the ideal mean-field model can reliably capture multicomponent diffusion effects among the dominant divalent cations in garnet.
Further insight comes from the study of Borinski et al. (2012) mentioned earlier, which generated experimental data for Fe and Mg in almandine-pyrope-rich garnet, and compared the use of ideal and non-ideal mean-field models. Although non-ideal models yielded superior fits to some measured experimental profiles, differences between the two models were inappreciable except in the case of unusual Mg-Mn-Ca-rich compositions. In those compositions, the experimental profiles contained 'inflections/changes in curvature (including effects of uphill diffusion) . . . that cannot be produced using an ideal model' (Borinski et al. 2012, p. 575).
Li: direct charge coupling with Y+REEs As described above, the incorporation of Li + into a dodecahedral site in garnet proceeds dominantly via the VIII LI substitution mechanism, in which charge balance relies on Y+REEs occupying dodecahedral sites adjacent to Li. This is the fundamental control on the nature of the SDPs for Li in garnet of the MLP aureole: absent this effect, Li would partition preferentially into product cordierite, and the SDP for Li would exhibit strong depletion towards the garnet rim, the opposite of what is observed. The diffusional implication of the VIII LI substitution mechanism is that the motion of Li ions becomes directly coupled to the motion of Y+REE ions, so that from a macroscopic viewpoint, the diffusing component is a Li + -(Y+REE) 3+ ion pair.
Because the concentration of Y+REEs in MLP garnet is more than an order of magnitude larger than that of Li, the forces driving the fluxes of Y+REEs will dominate, and Li will diffuse down the Y+REE gradient, penetrating inward at the same rate as the Y+REEs. Thus the SDP for Li should be indistinguishable from one generated by the ideal mean-field theory with Li self-diffusivity equal to that of Y+REEs. Figure 5 provides an example of a model based on this premise, and documentation of the quantitative equivalence of diffusivities for Li, Y, and Yb in MLP garnet is found in Fig. 3 of Cahalan et al. (2014). The requirement of charge compensation produces a direct coupling between Li + and Y+REEs, making Li in garnet several orders of magnitude more resistant to diffusive redistribution than Li in other common minerals (Cahalan et al. 2014, their Fig. 4).

L-REEs: indirect charge coupling among all Y+REEs
The hypothesis advanced above for the origin of maxima in SDPs for L-REEs, like that shown for Sm in Figure 3, invokes charge coupling among all Y+REEs via menzerite-type substitution mechanisms, the dominant means of incorporation of Y+REEs at levels in excess of Li + and Na + contents. A key consequence of the menzerite-type substitution is that wherever total Y +REE concentration is high, the concentration of VI (Mg, Fe) is equally elevated, and interchange of one Y+REE atom for another is possible by simple homovalent substitution. Furthermore, local charge balance requires the co-diffusion of VI (Mg, Fe) 2+ and VIII (Y+REE) 3+ as an ion pair, so the diffusion of each Y+REE element becomes dependent on the diffusion of all others; Distance from Center (cm) Figure 4. Fit of numerical model (solid curve) to measured Ca profile of Figure 1 (circles). Model implements ideal mean-field theory for multicomponent diffusion in the quaternary system Mg-Fe-Mn-Ca. Dashed horizontal line marks the uniform concentration used as the initial pre-resorption condition for the model; its length indicates the original radius of the homogeneous crystal. The model computes small increases above the original concentration nearly to the core of the crystal, a distance corresponding to the limit of inward diffusive penetration of the compositionally dominant components (Fe and Mg). After Carlson (2006); used by permission.
they are linked indirectly through the movement of VI (Mg, Fe). This concept was tested in a series of three progressively more sophisticated attempts to replicate the Sm profile shown in Figure 3.

Mean-field approach
Because the mean-field model of Carlson (2012) was restricted to a five-component system (e.g. Mg-Fe-Mn-Ca-Y or Mg-Fe-Mn-Ca-Sm), it was inherently incapable of capturing multicomponent effects among multiple Y+REEs. Here, that shortcoming is addressed by employing a six-component ideal mean-field model (Mg-Fe-Mn-Ca-Y-Sm). As indicated by the dash-dot curve in Figure 6, the best-fit mean-field model fails to generate a perceptible maximum. In addition, in order to obtain a concentration matching the measurement at the outer edge of the relict garnet, it was necessary to introduce the arbitrary (though defensible) assumption that diffusive exchange at the garnet rim continued tõ 740°C, which is~50°C below the temperature at which resorption ceased. Otherwise, modelled final rim concentrations significantly exceed the measured value.
Indirect charge coupling: linear coupling function A closer approximation to the measured profile comes from an enhanced model for the six-component system Mg-Fe-Mn-Ca-Y-Sm, based upon a simple form of indirect charge coupling via the menzerite-type substitution. The essential concept is that in regions of high Y content, where VI (Mg, Fe) is also elevated, the entropy of mixing will drive simple homovalent exchange of Sm for Y. Thus any flux of Y will be accompanied by a corresponding flux of Sm, as a consequence of their shared requirement for charge compensation by VI (Mg, Fe).
The flux of Sm is then the result of two competing driving forces. The dominant force is the one captured by the mean-field theory; it depends principally upon the concentration gradient for Sm, adding in very minor effects from the cross-terms coupling to the gradients of the other five species in the model. The strong preferential partitioning of Sm into product phases in the reaction zone produces a sharp drop-off in concentration at the garnet rim; this gradient leads to an overall outward flux of Sm, and yields a monotonic rimward decrease in its concentration, as seen in the mean-field model just presented. Acting counter to this mean-field driving force is one that results from charge coupling due to the menzerite substitution, which results in a flux of Sm that mimics at lower levels the flux of Y, which partitions into the  Figure 6. Fits of three numerical models (curves) to measured Sm profile of Figure 3 (filled circles with error bars). All models were computed for the six-component system Mg-Fe-Mn-Ca-Y-Sm; all used the same initial concentration profile, flat in the crystal's interior, and with a very gentle outward decrease towards the rim of the original pre-resorption crystal. The dash-dot curve shows the best fit of a model implementing the ideal mean-field theory; it fails to produce an appreciable maximum along the profile. The dashed curve shows the best fit of a model implementing the linear coupling of the Sm flux to the Y flux expressed by Equation (4); it yields a maximum of the appropriate magnitude, but overestimates Sm concentrations rimward of that maximum. The solid curve shows the best fit of a model implementing the quadratic coupling of the Sm flux to the Y flux expressed by Equation (5); it matches measured concentrations along the length of the profile. garnet, creating an inward diffusional flow. The net flux of Sm is the sum of these two competing effects. This model initially computes all fluxes under the mean-field approximation, then adds to the Sm flux a term that is simply proportional to the computed meanfield Y flux. Thus the model posits that the additional Sm flux due to indirect charge coupling is a linear function of the Y flux, which introduces a parameter k 1 as the proportionality factor between the two: in which the superscript mf designates a quantity computed in the mean-field approximation. To maintain material balance, the mean-field Y flux is reduced by an amount equal to that added to the Sm flux; this expresses the fact that proposed mechanism involves homovalent exchange of Sm for Y. Because of the much higher concentrations of Y, however, this adjustment has a negligible effect on the evolution of the Y profile.
The best fit of this model to the target profile is shown by the dashed curve in Figure 6; note that as with the mean-field model, to match the measured rim concentration, it was necessary again to assume that diffusive exchange at the garnet rim continued to~740°C, which is~50°C below the temperature at which resorption ceased. The result shows that for a suitable choice of k 1 , a maximum of the correct magnitude is generated. However, when using this approach, the best fits to this and other profiles are all systematically in error: fits that match both the measured Sm concentrations along the inner rise towards the maximum and also the outermost rim concentration invariably overestimate Sm concentrations between the maximum and the rim, where Y concentrations are highest. This suggests that the coupling effect is not as simple as a direct proportionality between the additional Sm flux and the Y flux. Instead it appears that the coupling is diminished in regions of very high Y concentration, which then requires a more complex form for the coupling function.
Indirect charge coupling: quadratic coupling function A plausible explanation for the reduction in coupling in regions of highest Y concentration comes from recognition that incorporation of Ya trivalent ion larger than the divalent ion for which it substitutesinduces significant lattice strain. Ions like Nd, Sm, and Eu, which are the largest of the REEs studied, will therefore find greater barriers to their own incorporation as the total concentration of Y +REEs rises, making their homovalent interchange for Y +REEs energetically more costly and therefore less likely.
To approximate this effect on the flux of Sm in the present example, Equation (4) can be modified by including a term that reduces the coupling as the flux of Y increases.
The simplest such modification would add on a quadratic term, one proportional to the square of the Y flux, introducing a second adjustable parameter k 2 to describe the magnitude of the Sm-exclusion effect due to lattice strain. The equation for the Sm flux then takes the form in which the sign of the last term is chosen as needed to produce a flux in the direction opposite to the flux of Y. 3 The mean-field flux of Y is again adjusted accordingly, to maintain material balance. The best fit of this model to the target profile is shown by the solid curve in Figure 6, which satisfactorily replicates the measured concentrations across the entire profile. In this approach, the fit at the rim of the relict garnet does not rely upon an arbitrary assumption of further rim equilibration after resorption ceases; instead, the match to the measured concentration follows directly from proper specification of values for k 1 and k 2 . Using this same crystal as an example, the evolution through time of its diffusion profiles for Y and Sm is illustrated in Figs. S1 (syn-resorption phase) and S2 (postresorption phase), included as supplemental online material.
Fits of quadratic coupling model to measured L-REE profiles A multicomponent diffusion model implementing indirect charge coupling of fluxes in the manner specified by Equation (5) was fitted to SDPs for Nd, Sm, and Eu measured on garnet in samples M04A1f and M20D1c, collected at distances from the intrusive contact of 450 and 1100 m, respectively. Values for k 1 and k 2 were iteratively adjusted to obtain optimal fits to the measured profiles; optimal values are listed in Table 1. The resulting fits are shown in Figure 7, together with corresponding fits to the Y profiles. The excellent agreement with the measured concentrations validates the concept of indirect diffusive charge coupling among Y +REEs, incorporating reduction of coupling effects due to lattice strain as total Y+REE concentrations rise.
For all three L-REEs, values for k 1 are similar in magnitude to the ratios of their molar concentration to that of Y; this accords well with the concept of flux linkage via homovalent exchange of L-REEs for Y. Considering Nd and Sm, which have similar concentrations: one might expect the values for k 2 to be slightly larger for Sm, owing to its larger size; however, such effects cannot be clearly discerned amidst the scatter among the few data available in Table 1.

Form of the coupling function
The mathematical form of the coupling function that appears in Equation (5) while empirically adequateis nonetheless somewhat unsatisfying from a physical standpoint, insofar as it lacks a direct connection to  the proposed coupling mechanism. Conceptually, the proposed mechanism links the L-REE flux to the Y flux by a proportionality factor that is larger at low Y concentrations and smaller at high Y concentrations, due to lattice strain. One might therefore envision a form for the coupling function resembling that of Equation (4), but with the quantity k 1 cast as a function of Y+REE concentrationor even the overall concentration, considering that the strain effects of incorporating Y+REEs are much larger in structures with smaller unit-cell dimension (those rich in almandine and pyrope) than they are in more expanded structures (those rich in spessartine and especially grossular). Unfortunately, the available data are insufficient to constrain adequately such an approach, so the purely empirical relationship expressed by Equation (5) must suffice at present. In a similar vein, it is conceivable that one might implement a physically more satisfying model by approaching the problem from a different angle: if and when activity coefficients for Y+REE components in garnet solutions are determined over a range of compositions and P-T conditions, perhaps a mean-field model that accounts directly for the non-ideality of the solution could succeed. Again, however, the absence of required data makes it premature to take such an approach now.
Despite the empirical nature of the functional form for coupling employed here, the fundamental finding is robust: coupling due to indirect charge compensation produces detectable effects on diffusional flows of Y+REEs.

Universal coupling among Y+REEs: abundance effects
The mechanism responsible for the observed indirect charge coupling is universal in the sense that all Y+REEs are diffusionally linked: the flux of any one of them will induce a partial flux of all the others. In the particular case of the SDPs examined here, the coupling is made evident for the L-REEs by the fact that they partition during reaction in a sense opposite to that of Y and the other REEs. Such coupling, however, must also occur among all Y+REEs, although it will be largely obscured for elements whose individual concentration gradients are everywhere coincident with the direction of the overall gradient in total Y+REE concentration.
The concept of homovalent interchange among diffusing Y+REEs implies that elements present in the greatest abundance will dominate the overall flow of the Y+REEs. Elements present at low levels, with correspondingly small concentration gradients, will be strongly affected by the contribution to their overall flux that is induced by their homovalent interchange with diffusing elements present at high levels, with correspondingly large concentration gradients. The converse will not occur: small fluxes of elements with low concentrations will induce only a negligible flux of an element with high concentration.
In garnet, Y dominates the Y+REE budget, with appreciable contributions from Yb, Er, and Dy. For example, in garnet of the MLP aureole, Y amounts to~60 mol % of the total Y+REE budget, and Y + Yb + Er + Dy makes up 90 mol %. This is the result of the strong propensity for inclusion of the smaller HREEs rather than the larger LREEs, combined with the fundamentally higher elemental abundance of Y and of REEs of even atomic number. The net effect is that the diffusional redistribution of Y (± HREEs) will exert an influence on the diffusional motion of all REEs, and induced fluxes may be a prominent contributor to the redistribution of elements present at low levels.

Energetics of charge coupled diffusive groups
Vital to the notion of multicomponent diffusional coupling due to charge compensation is the idea that pairs of cations can move individually through a crystal structure, yet remain associated closely enough so that their reciprocal excess and deficiency of charge yields electroneutrality at the local scale. The probability is negligible that two vacancies will be positioned so that both ions in a pair can move synchronously from one set of adjacent polyhedra into a new set of adjacent polyhedra. Instead, the motions of the ions must be sequential but linked, which means that to move through the structure, they must be transiently separated from one another before recombining in their lowest-energy configuration in adjacent polyhedra.
Such motions are consistent with lattice-dynamics calculations for the relative energies of the ions at varying degrees of separation. Carlson et al. (2014, p. 1032-1033) evaluated the free energies at 1200 K, 2 GPa for ion pairs corresponding to the alkali-type and menzerite-type substitution mechanisms (along with others) in garnet end-members pyrope, almandine, spessartine and grossular; in all cases, the results were quite similar. Dissociated ions, those separated by a distance sufficient to produce only negligible interactions between them, were found to have free energies 200 kJ/mol higher than ion pairs occupying adjacent polyhedra, the minimum-energy configuration. This indicates that the ions will tend to remain strongly associated, producing local charge compensation, even at high temperature. Moreover, the minimum-energy configurations were found to be about 10-30 kJ/mol lower in energy than alternative arrangements in which the charge-compensating ions are progressively more distant from one another, that is, when occupying second-, third-or fourth-nearest-neighbour polyhedra. Considering that the thermal energy RT at 1200 K is 10 kJ/mol, an appreciable fraction of the ion pairs can experience transient occupancy of these slightly higherenergy configurations, as required for the coupled movement of charge-compensating groups during diffusional transport of Y+REEs.

Implications and conclusion
The ideal mean-field theory does not account for multicomponent diffusion effects arising from coupled motions of the ion pairs that provide charge compensation for heterovalent substitutions in crystalline solutions. Ion pairing can couple diffusional flows either directly, as exemplified here by the co-diffusion of Li + and (Y+REE) 3+ , or indirectly, as exemplified here by the induced fluxes of L-REEs generated by homovalent interchange with more abundant Y+REEs. Direct coupling effects can be modelled by matching the diffusivity of the subordinate ion to that of the dominant ion. Indirect coupling effects are more complex, but at least in some cases (as for L-REEs here) empirical approaches may be taken that modify the flux of the subordinate ion by adding terms that account for coupling via homovalent interchange, and for variation in the strength of that coupling as the flux (or concentration) of the dominant ion changes (e.g. Equation (5)).
The specific examples offered here of diffusional coupling due to charge compensation have clear petrologic implications. First, the direct coupling of Li and Y +REEs implies that the mobility of Li can be no greater than that of the Y+REEs, which diffuse in garnet much more slowly than the dominant divalent cations (Carlson 2012). Consequently, retention of Li zoning in garnet will persist to substantially higher temperatures than in most other minerals (Cahalan et al. 2014), making garnet an especially robust monitor of the behaviour of Li in deep crustal and upper-mantle systems. Because the concentrations and isotopic compositions of Li in minerals record fluid-rock interactions during petrogenesis (e.g. Tomascak et al. 2000;Magna et al. 2006;Ionov and Seitz 2008;Penniston-Dorland et al. 2010, 2012, garnet will have particular value in such studies, due to the relative diffusional stability of Li within it. Second, the indirect charge coupling among all Y+REEs could have important implications for geochronological applications involving garnet. Variable degrees of diffusional coupling of Y with Lu, Sm, and Nd may lead to differential redistribution of those elements during a diffusive episode. Especially when they are combined with partial resorption, these effects can modify apparent ages and may affect the variations in ages obtained by core-to-rim microsampling (cf. Pollington and Baxter 2010;Kelly et al. 2011;Bloch and Ganguly 2015). In addition, because REE systematics in garnet and co-existing phases can sometimes be used to infer fluid histories (e.g. Sorensen et al. 2006), such studies in high-temperature rocks should not overlook the potential for modification of Y+REE distributions by diffusion induced by multicomponent effects.
Recognition of multicomponent coupling effects should inform design of experimental attempts to measure rates of diffusion of monovalent and trivalent ions in garnet. Synthetic experimental systems in which nearend-member garnet is infused with a single impurity component (e.g. Ganguly et al. 1998b;Cherniak 2005;Tirone et al. 2005) will not capture the diffusional effects of cross-coupling among multiple components that are observed in natural systems. Moreover, experiments designed without accounting explicitly for the substitution mechanism(s) responsible for incorporation of the ions in heterovalent substitutions may yield conflicting or spurious results, as they will be subject to diffusional variations that arise from coupling via different types of substitution mechanisms operating at variable levels of concentration. The use in such experiments of carefully characterized natural materials of greater compositional complexity would provide an important control on the petrologic applicability of results from synthetic systems.
This study documents the operation in nature of two types of charge-coupled multicomponent diffusional effects in garnet, both arising from the requirement for local compensation of the electrostatic imbalance that accompanies any heterovalent substitution. Although this work offers conceptual explanations for the origin of these effects that are consistent with data from measured stranded diffusion profiles, its treatment of induced flows due to indirect coupling remains empirical, lacking full theoretical rigour. Development of a more comprehensive approach with broader quantitative predictive ability will likely require constraints from experiments directed at obtaining more precise data across a wider range of phenomena, perhaps combined with the casting of those data into a non-ideal extension of the conventionally applied ideal mean-field theory. Notes 1. Although incorporation mechanisms for alkali ions in garnet have not yet been systematically explored, the lattice-dynamics calculations of Carlson et al. (2014) and the analytical results of Enami et al. (1995) indicate that Li and Na are likely to be incorporated almost exclusively via mechanisms involving Y+REEs on neighbouring dodecahedral sites.
2. This result was anticipated, insofar as the ideal meanfield theory predicts that a single component present in tracer amounts always diffuses independently of other components (Lasaga 1979, p. 458). 3. An additional requirement of the numerical simulation is that for Sm, the computed flux out of any region can be no larger than that allowed by the abundance of Sm there, to prevent the physically unrealistic possibility of computing a negative concentration of Sm. Thus the Sm flux in any time step is limited to that which would deplete the region completely in that step.