Chemical Heterogeneities in the Mantle: Progress Towards a General Quantitative Description

. Chemical equilibration between two different assemblages (peridotite-type and gabbro/eclogite-type)has been determined using basic thermodynamic principles and certain constraints and assumptions regarding mass and reactions exchange. When the whole system (deﬁned by the sum of the two sub-systems) is in chemical equilibrium the two assemblages will not be homogenized but they will preserve distinctive chemical and mineralogical differences. Furthermore, the mass transfer be-tween the two sub-systems deﬁnes two petrological assemblages that separately are also in local thermodynamic equilibrium. 5 In addition, when two assemblages previously equilibrated as a whole in a certain initial mass ratio are held together assuming a different proportion, no mass transfer occurs and the two sub-systems remain unmodiﬁed. By modeling the chemical equilibration results of several systems of variable initial size and different initial composition it is possible to provide a quantitative framework to determine the chemical and petrological evolution of two assemblages from an initial state, in which the two are separately in chemical equilibrium, to a state of equilibration of the whole system. Assum-10 ing that the local Gibbs energy variation follows a simple transport model with an energy source at the interface, a complete petrological description of the two systems can be determined over time and space. Since there are no data to constrain the kinetic of the processes involved, the temporal and spatial scale is arbitrary. The evolution model should be considered only a semi-empirical tool that shows how the initial assemblages evolve while preserving distinct chemical and petrological features. Nevertheless despite the necessary simpliﬁcation, a 1-D model illustrates how chemical equilibration is controlled by the size 15 of the two sub-systems. By increasing the initial size of the ﬁrst assemblage (peridotite-like), the compositional differences between the initial and the ﬁnal equilibrated stage become smaller, while on the eclogite-type side the differences tend to be

mined using basic thermodynamic principles and certain constraints and assumptions regarding mass and reactions exchange.When the whole system (defined by the sum of the two sub-systems) is in chemical equilibrium the two assemblages will not be homogenized but they will preserve distinctive chemical and mineralogical differences.Furthermore, the mass transfer between the two sub-systems defines two petrological assemblages that separately are also in local thermodynamic equilibrium.
In addition, when two assemblages previously equilibrated as a whole in a certain initial mass ratio are held together assuming a different proportion, no mass transfer occurs and the two sub-systems remain unmodified.
By modeling the chemical equilibration results of several systems of variable initial size and different initial composition it is possible to provide a quantitative framework to determine the chemical and petrological evolution of two assemblages from an initial state, in which the two are separately in chemical equilibrium, to a state of equilibration of the whole system.Assuming that the local Gibbs energy variation follows a simple transport model with an energy source at the interface, a complete petrological description of the two systems can be determined over time and space.Since there are no data to constrain the kinetic of the processes involved, the temporal and spatial scale is arbitrary.The evolution model should be considered only a semi-empirical tool that shows how the initial assemblages evolve while preserving distinct chemical and petrological features.
Nevertheless despite the necessary simplification, a 1-D model illustrates how chemical equilibration is controlled by the size of the two sub-systems.By increasing the initial size of the first assemblage (peridotite-like), the compositional differences between the initial and the final equilibrated stage become smaller, while on the eclogite-type side the differences tend to be larger.A simplified 2-D dynamic model in which one of the two sub-systems is allowed to move with a prescribed velocity, shows that after an initial transient state, the moving sub-system tends to preserve its original composition defined at the influx side.The composition of the static sub-system instead progressively diverges from the composition defining the starting assemblage.The observation appears to be consistent for various initial proportions of the two assemblages, which simplify somehow the development of potential tools for predicting the chemical equilibration process from real data and geodynamic applications.
Four animation files and the data files of three 1-D and two 2-D numerical models are available following the instructions in the supplementary material.

Introduction
Our understanding of the Earth and planetary interiors is based on the underlying assumption that thermodynamic equilibrium is effectively achieved on a certain level, which means that the system under consideration is in thermal, mechanical and chemical equilibrium within a certain spatial and temporal domain.Although this may appear to be just a formal definition, it affects the significance of geophysical, petrological and geochemical interpretations of the Earth's interior.While the assumption of thermodynamic equilibrium is not necessarily incorrect, the major uncertainty is the size of the domain on which the assumption is expected to be valid.
The Earth and planetary interior as a whole could be defined to be in mechanical equilibrium when the effect of the gravitational field is compensated, within a close limit, by a pressure gradient (for simplicity variations of viscous forces are neglected).Even when this is effectively the internal state (one example could be perhaps the interior of Mars), thermodynamic equilibrium most likely is not achieved because it requires also thermal equilibrium (i.e.uniform temperature) and chemical equilibrium (for possible definitions of chemical equilibrium see for example Prigogine and Defay, 1954;Denbigh, 1971;Smith and Missen, 1991;Kondepudi and Prigogine, 1998).On a smaller scale instead, local thermodynamic equilibrium could be a reasonable approximation.If the system is small enough, the effect of the gravitational field is negligible and a condition close to mechanical equilibrium is achieved by the near balance between the gravitational force and pressure (locally both density and pressure are effectively uniform and viscous forces are neglected for simplicity).Clearly a perfect balance will lead to static equilibrium.On the other end dynamic equilibrium makes harder for chemical and thermal equilibrium to be maintained.
In studies of planetary solid bodies it is often reasonable to assume dynamic equilibrium close to a quasi-static condition in which the forces balance is close but not exactly zero.At a smaller scale it is then easier to consider that the temperature is also nearly uniform.The main uncertainty remains the chemical equilibrium condition.On a planetary scale, whether the size of system under investigation is defined to be on the order of hundreds of meters or few kilometers, it has little effect on the variation of the gravity force and in most cases on the temperature gradient.But for chemical exchanges, the difference could lead to a significant variation of the extent of the equilibration process.For the Earth's mantle in particular this is the case because it is generally considered to be chemically heterogeneous.The topic has been debated for some time (Kellogg, 1992;Poirier, 2000;Schubert et al., 2001;van Keken et al., 2002;Helffrich, 2006) and large scale geodynamic models to study chemical heterogeneities in the Earth's mantle have been refined over the years (Gurnis and Davies, 1986;Ricard et al., 1993;Christensen and Hofmann, 1994;Walzer and Hendel, 1999;Tackley and Xie, 2002;Zhong, 2006;Huang and Davies, 2007;Brandenburg et al., 2008;Li and al., 2014;Ballmer et al., 2015Ballmer et al., , 2017)).Geochemical (van Keken and Ballentine, 1998;van Keken et al., 2002;Kogiso et al., 2004;Blusztajn et al., 2014;Iwamori and Nakamura, 2014;Mundl et al., 2017) and geophysical (van der Hilst et al., 1997;Trampert et al., 2004;Tommasi and Vauchez, 2015;Zhao et al., 2015;Tesoniero et al., 2016) data essentially support the idea that the mantle develops and preserves chemically heterogeneities through the Earth's history.Even though all the interpretations of the mantle structure are based on the assumption of local thermodynamic equilib-rium, the scale of chemical equilibration has never been investigated in much detail.An early study (Hofmann and Hart, 1978) suggested that chemical equilibrium cannot be achieved over a geological time, even for relatively small systems (kilometer scale), hence it must preserve chemical heterogeneities on the same scale.The conclusion was inferred based on volume diffusion data of Sr in olivine at 1000 o C. At that time the assessment was very reasonable, albeit the generalization was perhaps an oversimplification of a complex multiphase multicomponent problem.In any case, significant progress in the experimental methodology to acquire kinetic data and better understanding of the mechanisms involved suggest that the above conclusion should be at least reconsidered.Based on the aforementioned study, the only mechanism that was assumed to have some influence on partially homogenizing the mantle was mechanical thinning/mixing by viscous deformation (Kellogg and Turcotte, 1987).In addition very limited experimental data on specific chemical reactions relevant to mantle minerals (Rubie and Ross II, 1994;Milke et al., 2007;Ozawa et al., 2009;Gardés et al., 2011;Nishi et al., 2011;Dobson and Mariani, 2014) came short to set the groundwork for a general re-interpretation of chemical heterogeneities in the mantle.
Perhaps a common misconception is that chemical equilibrium between two lithologies implies chemical homogenization.In other words, if the mantle is heterogeneous, chemical equilibration must have not been effective.This is not necessarily true.A simple example may explain this point.If we considering for example the reaction between quarz and periclase to form variable amount of forsterite and enstatite: would require the formation of a bimineralic single layer made of a mixture of enstatite and forsterite crystals.However experimental studies (e.g.Gardés et al., 2011) have shown instead the formation of two separate monomineralic layers, one made of policrystalline enstatite and the other one made of forsterite.
In summary there are still unanswered questions regarding the chemical evolution of the Earth's mantle, for example, at what spatial and temporal scale we can reasonably assume that a petrological system is at least close to chemical equilibrium?and how does it evolve petrologically and mineralogically?
This study expands a previous contribution that aimed to provide an initial procedure to determine the chemical equilibration between two lithologies (Tirone et al., 2015).The problem was exemplified in a illustration (figure 1 in Tirone et al. (2015)).
Because certain assumptions need to be made, the heuristic solution, further developed here, is perhaps less rigorous than other approaches based on diffusion kinetics that were applied mainly for contact metamorphism problems (Fisher, 1973;Joesten, 1977;Nishiyama, 1983;Markl et al., 1998).However the advantage is that it is relatively easy to generalize, and it leads towards a possible integration with large scale geodynamic numerical models while still allowing for a comparison with real petrological data.At the same time it should be clear that to validate this model approach and to coinstrain the extent of the chemical equilibration process, experimental data should be acquired on the petrological systems investigated here and in the previous study.
The following section (section 2) outlines the revised procedure to determine the two petrological assemblages forming together a system in chemical equilibrium.The revision involves the method used to determine the composition of the two assemblages when they are in equilibrium together, the database of the thermodynamic properties involved and the number of oxides considered in the bulk composition.In addition since the solids are non-ideal solid mixtures (in the previous study all mixtures were ideal), the chemical equilibration requires that the chemical potential of the same components in the two assemblages must be the same (Prigogine and Defay, 1954;Denbigh, 1971).The method is still semi-general in the sense that a similar approach can be used for different initial lithologies with different compositions, however some assumptions and certain specific restrictions should be modified depending on the problem.The simplified system discussed in the following sections assumes on one side a peridotite-like assemblage, and a gabbro/eclogite on the other side.Both are considered at a fixed pressure and temperature (40 kbar and 1200 o C) and their composition is defined by nine oxides.The general idea is to conceptually describe the proxy for a generic section of the mantle and a portion of a subducting slab.A more general scheme that allows for variations of the pressure and temperature should be considered in future studies.The results of the equilibration method applied to 43 different systems are presented in section 2.1.The parameterization of the relevant information that can be used for various applications is discussed in section 2.2.Section 3 presents the first application of a 1-D numerical model applied to pairs of assemblages in variable initial proportions to determine the evolution over time towards a state of equilibration for the whole system.The next section (4) illustrates the results of few simple 2-D dynamic models that assume chemical and mass exchange when one side moves at a prescribed velocity while the other side remains fixed in space.These simple models only serve the purpose to illustrate how distinct mineralogical and petrological features are preserved after chemical equilibration has been reached.
All the necessary thermodynamic computations are performed in this study with the program AlphaMELTS (Smith and Asimow, 2005), which is based on the thermodynamic modelization of Ghiorso and Sack (1995); Ghiorso et al. (2002) for the melt phase, the mixture properties of the solid and certain end-member solids.The thermodynamic properties of most of the endmember solid phases are derived from an earlier work (Berman, 1988).Even though melt is not present at the (P,T,X) conditions considered in this study, and other thermodynamic models are also available (Saxena, 1996;Stixrude and Lithgow-Bertelloni, 2005;Piazzoni et al., 2007;de Capitani and Petrakakis, 2010;Holland and Powell, 2011;Duesterhoeft and de Capitani, 2013), AlphaMELTS proved to be a versatile tool to illustrate the method described in this work.It also allows for a seamless transition to potential future investigations in which it would be possible to study the melt products of two equilibrated, or partially equilibrated assemblages when the P,T conditions are varied.

Modeling Chemical Equilibration Between Two Assemblages
This section describes in some details the procedure to determine the transformations of two assemblages after they are put in contact and the system as a whole reaches a condition of chemical equilibrium.The bulk composition is described by nine Retaining the input format of the AlphaMELTS program, the bulk composition is given in grams.Pressure and temperature are defined at the beginning of the process and they are kept constant.Water (thermodynamic phase) is not considered simply because the mobility of a fluid phase (or melt) cannot be easily quantified and incorporated in the model.Three independent equilibrium assemblages are retrieved using AlphaMELTS.These are standard equilibrium computations which consist of solving a constrained minimization of the Gibbs free energy (van Zeggeren and Storey, 1970;Ghiorso, 1985;Smith and Missen, 1991).The first two equilibrations involve the bulk compositions of the two assemblages separately.The third one is performed assuming a weighted average of the bulk composition of the two assemblages in a predefined proportion, for example 1:1, 5:1 or 100:1, also expressed as f:1 where f=1,5,100 (peridotite : gabbro/eclogite).This third computation applies to a whole system in which the two assemblages are now considered sub-systems.The variable proportion essentially allows to put increasingly larger portions of the subsystem mantle in contact with the sub-system gabbro/eclogite using the factor f to indicate the relative "size" or mass of material involved.By using AlphaMELTS the mineralogical abundance and composition in moles is retrieved from the file phase_main_tbl.txt,while the chemical potential for each mineral component in the solid mixture is retrieved from the thermodynamic output file (option 15 in the AlphaMELTS program).Knowing all the minerals components involved, an independent set of chemical reactions can be easily found (Smith and Missen, 1991).For the problem in hand, the list of minerals and abbreviations are reported in table 1, and the set of independent reactions are listed in table 2.
Given the above information, the next step is to determine the bulk composition and the mineralogical assemblages of the two sub-systems after they have been put together and equilibration of the whole system has been reached.For this problem the initial amount of moles n of mineral components i in the two assemblages is allowed to vary (∆n i ), provided that certain constraints are met.The set of constraints can be broadly defined in two categories.The first group consist of relations that are based on general mass, chemical or thermodynamic principles.The second set of constraints are based on certain reasonable assumptions that should be verified by future experimental studies.
The first and most straightforward set of constraints requires that the sum of the moles in the two assemblages should be equal to the moles of the whole system: where n i (A 0 ) represents the initial number of moles of the mineral component in the first assemblage (A) in equilibrium before it is put in contact with the second assemblage (B).A similar definition applies to n i (B 0 ).∆n i (A) and ∆n i (B) are the variations of the number of moles after the two assemblages are held together and n i (W ) is the number of moles of the component in the whole assemblage (A + B).The size of the whole assemblage is defined by f + 1 where f refers to the size of the first assemblage.
Another set of constraints imposes the condition of local chemical equilibrium (Prigogine and Defay, 1954;Denbigh, 1971;Kondepudi and Prigogine, 1998) by requiring that the chemical potentials of the mineral components in the two sub-systems cannot differ from the chemical potentials found from the equilibrium computation for the whole assemblage (W ): where µ i (A) is the chemical potential of the mineral component in the assemblage A whose number of moles is n i (A) = n i (A 0 ) + ∆n i (A), and a similar expression for the second assemblage B.
Another constraint is given by the sum of the Gibbs free energy of the two sub-systems that should be equal to the total Gibbs free energy of the whole system: where G(A) = i n i (A)µ i (A) and similar equations for B and W .
The list of reactions in table 2 allows to define a new set of equations which relates the extent of the reaction ξ r with the changes of the moles of the mineral components (Prigogine and Defay, 1954;Kondepudi and Prigogine, 1998) where all the extent of the reactions are considered to be potential new variables.However not necessarily all the ξ r should be treated as unknowns.This can be explained by inspecting for example table 3, which provides the input data and the results of the equilibrium modeling of on of the study cases, in particular the one that assumes an initial proportion 1:1 (f=1 which is uniquely coupled to ∆n OEss (A), furthermore ξ (T −4) coupled to ∆n OHd (A), also ξ (T −11) coupled to −∆n OJd (A), and finally ξ (T −17) fixed by ∆n Coe (B).
For the problem in hand the above set of relations does not allow to uniquely define the changes of the moles of the mineral components in the two sub-systems.Therefore additional relations based on some reasonable assumptions have been added to the solution method.Future experimental studies will need to verify the level of accuracy of such assumptions.Certain constraints on the mass exchange can be imposed by comparing the equilibrium mineral assemblage of the whole system (W ) with the initial equilibrium assemblages in A 0 and B 0 .For example table 3 shows that olivine is present in the whole assemblage W . However initially olivine is only located in sub-system A 0 .Therefore rather than forming a complete new mineral in B, the assumption is that the moles of fayalite (Fa), monticellite (Mtc) and forsterite (Fo) will change only in sub-system A to comply with the composition found for the whole assemblage W .Following this reasoning the changes in the two sub-systems could be set as: ∆n F a (A) = 0.0008090, ∆n Mtc (A) = −0.0000555and ∆n F o (A) = −0.0726300and In this particular case the same assumption is also applicable to the orthopyroxene components.It is clear that starting with different bulk compositions or proportions or (T,P) conditions, alternative assemblages may be formed, therefore different conditions may apply, but the argument on which the assumption is based should be similar.
Additional constraints based on further assumptions can be considered.For example, garnet appears on both sides A 0 and B 0 .
The components pyrope (Prp) and grossular (Grs) contribute only to two reactions, (T-1) and (T-12), and in both cases the reactions involve only olivine components which have been fixed in sub-system A, as previously discussed.The assumption that is made here is that the change of the moles of the garnet components in sub-system B will be minimal because no olivine is available in this sub-system.Therefore the following relation is applied: and similar relations can be also imposed to the other garnet components, Alm and Grs.The same argument can be applied to the clinopyroxene and spinel components.For example the spinel component hercynite (Hc) appears only in reaction (T-13), which involves olivine and orthopyroxene components (Fa, ODi) located in sub-system A, and the garnet component Alm which has been already defined by the previous assumption.
The overall procedure is implemented with the use of Minuit (James, 1994), a program that is capable of performing a minimization of multi-parameter functions.Convergence is obtained making several calls of the Simplex and Migrad minimizers (James, 1994).The procedure is repeated with different initial values for the parameters ∆n i (A), ∆n i (B) and ξ r to confirm that a unique global minimum has been found.

Results of the Chemical Equilibrium Model Between Two Assemblages
This procedure described in the previous section has been applied to 43 different cases, varying the proportion of the two subsystems from 1:1 to 1000:1 and considering different, but related, initial compositions.The initial bulk composition and the proportion factor f of the two sub-systems for all the 43 cases are included in a table available in the supplementary material.
For example the initial compositions for A 0 and B 0 applied to case #11 are taken from table 4 (column A * ) and from table 3 (column B 0 ), both tables discussed in this section.Tables 3-7 report the results of the procedure discussed in the previous section for few cases.Table 3 was briefly introduced earlier to show the initial bulk composition of the two sub-systems (upper portion of the table), the initial equilibrium assemblages and the mole changes after the chemical equilibration (lower part of the table ).The table also includes the bulk composition in the two sub-systems after the chemical equilibration procedure is completed (upper part, column 5 and 6).These bulk compositions are calculated from the mole abundance of the mineral components shown in the lower part (columns 4 and 7).The total mass of the sub-systems is reported as well.Note that negative abundance of certain mineral components is permissible according to the thermodynamic model developed by Ghiorso (Ghiorso and Carmichael, 1980;Ghiorso, 2013) as long as the related oxides bulk abundance is greater than zero.
In the example shown in table 3 there is a significant mass transfer from B to A: mass(A 0 )=100, mass(A)=146.36and mass(B 0 )=100, mass(B)=53.64 (grams).The table also includes the total Gibbs energy for the sub-systems, before and after the equilibration of the whole system which are computed from the output of the program AlphaMELTS after combining the moles of the components and the relative chemical potentials.The total Gibbs free energy is relevant for the parameterization discussed in the next section.Table 4 is a summary of a further analysis aiming to investigate whether there is any pattern in the compositions of the two sub-systems.The bulk compositions in the upper portion of the table (A * , B * ) are obtained by normalizing the oxides in A and B (upper part, column 5 and 6 of table 3) to a total mass of 100 grams.
For example SiO 2 in A * from table 4 (47.434) is 100×(SiO 2 in A)/(sum of oxides in A) from table 3, which is equal to 100×69.428/146.367.The normalized oxides (A * , B * ) represent the mass of the components in grams when the total mass is 100 grams, which is obviously also equivalent to the weight % of the components.These bulk compositions can be used for two new Gibbs free energy minimizations, one for each of the two sub-systems, to retrieve the correspondent equilibrium assemblages separately.The interesting observation that can be made following the summary in the lower part of table 4, is that the abundance of the mineral components remains unmodified after scaling the results for the total mass of the system.For example using the data from table 3, the proportion relation: n alm (A) : 146.347 = n alm (A * ) : 100 gives n alm (A * ) = n alm (A) × 100/146.347= 0.01453 × 0.6833 = 0.009928 which is remarkably close to the moles of almandine found from the separate equilibration calculation reported in table 4, n alm (A * ) = 0.0099353.In other words the scaling factor used to define the input oxide bulk composition can be also applied to the equilibrium mineral assemblage.
Based on this observation, some equilibration models have been carried out considering at least one of the initial composition from a previous model (e.g.A * from a previous equilibration model ⇒ input for a new model A 0 or alternatively B * ⇒ B 0 ), while for the other sub-system the initial bulk composition from table 3 is used again.A special case is the one shown in table 5 in which both A 0 and B 0 are taken from the equilibrated and normalized data of the previous model, A * and B * , reported in table 4. If the proportion in the new model remains the same, 1:1, then clearly no compositional changes are expected since the whole system is already in equilibrium.If the proportion is changed, for example to 5:1 (f = 5), the bulk composition of the whole system is different from the bulk composition of the whole system with 1:1 proportion and the assemblages in the two sub-systems may not remain unmodified after equilibration.However this does not appear to be the case, as shown in table 5, where ∆n i (A) and ∆n i (B) are very small.The results suggest that the moles of the mineral components remain unchanged.
A more general case with f = 5 is presented in table 6.The model is essentially the same shown in table 3, but with proportion of the two initial sub-systems set to 5:1.As expected the results of the equilibration process are different from the results starting with an initial proportion 1:1 (table 3).For example with 1:1, n alm (A) = 0.01453, while with 5:1, n alm (A)/5 = 0.00737.
The question is whether the observation made for the first studied case with proportion 1:1 can be generalized.In particular the observation that the minerals abundance in the two sub-systems from the equilibration procedure of the whole system is equivalent to the one that is obtained from two separate equilibration computations using the normalized bulk compositions A * and B * .Indeed it appears that the same conclusion can be made for the model with 5:1 initial proportion (table 7).The number of moles of the almandine component is (n alm (A)/5) × 100/110.064= 0.006698 (table 6) which can be compared with n alm (A * ) = 0.006695 from table 7. The similarity has been also observed for all the other models with f ranging from 1 to 1000.

Parameterization of the Equilibrium Model Results for Applications
While interesting observations have been made about the mineralogical assemblages in the two sub-systems after chemical equilibration, it is still unclear how this type of model can be applied for studies on the chemical evolution of the mantle.
Figure 1 summarizes the relevant data that allows to determine the bulk composition and the mineralogical assemblage in the two sub-systems after the chemical equilibration process is completed.
The key quantity is the normalized Gibbs energy of the two sub-systems after they have been equilibrated, G(A * ) and G(B * ).
The normalized Gibbs energy for an unspecified sub-system (either A * or B * ) is defined by the symbol G( * ).The quantity can be computed from the AlphaMELTS output after the Gibbs free energy minimization is applied to A * or B * , or it can be simply obtained by scaling G(A) or G(B).Panel 1-A) shows the relation between the ratio G(A * )/G(B * ) and G(B * ) which will be used later to define G( * ) at the interface between the two assemblages.The data in the figure for the 43 models have been fitted using a Chebyshev polynomials (Press et al., 1997).By knowing G( * ), it is possible to retrieve the abundance of all the oxides defining the bulk composition normalized to 100 grams.An example is shown in panels 1-B) and 1-C) which illustrate the data points for M gO in (A * ) and (B * ) in the 43 study models and the fitting of the points using Chebyshev polynomials.
The mass transfer between the two sub-systems can be related to the total Gibbs free energy variation in each of the two subsystems G(A) and G(B).The two relations are almost linear, as shown in panel 1-D).For practical applications, once a relation is found between G and the normalized G( * ), then the mass transfer can be quantified.Panel 1-E) of figure 1 shows the data points and the data fitting with the Chebyshev polynomial of the function More details on the use of the fitting polynomial functions are provided in the next section.

Application to the evolution of a 1-D Static Model with Variable Extension
The chemical and petrological evolution of two assemblages can be investigated with a 1-D numerical model, assuming that the two sub-systems remain always in contact and they are not mobile.The problem is assumed to follow a simple conduction/diffusion couple-type model with variable size for the local variation of G( * ) which can be expressed by the following equation for each sub-system: where S( * ) is a scaling factor and G( * ) and S( * ) refers to either A * or B * .Time t, distance d x ( * ) and the scaling factor S( * ) have no specific units since we have no knowledge of the kinetic of the processes involved.At the moment these quantities are set according to arbitrary units, S(A*) and S(B*) are set to 1, while t, d x (A * ) and d x (B * ) have different values depending on the numerical simulation.It should be clear that the dynamic model provides only a semi-empirical quantitative description of a complex process.The main purpose is to illustrate the general concept and to show that the two assemblages could develop distinct regions evolving towards the condition of chemical equilibrium, while far from the interface area the initial compositions can be preserved for a certain amount of time.The detailed description on how the two sub-systems will eventually reach chemical equilibration is beyond the scope of this study.
The numerical solution with grid spacing ∆d x ( * ), uniform on both sides, is obtained using the well-known Crank-Nichols method (Tannehill et al., 1997).At the interface (defined by the symbol if ) the polynomial of the function shown in panel 1-A) of figure 1 is used together with the flux conservation equation: G(A * ) l and G(B * ) l define the outside boundary limits of the whole system which represent either the closed-end of the system or the middle point of two mirrored images.
To determine the mass transfer and how it affects the length of the two sub-systems, the following steps are applied.The polynomial of the relation shown in panel 1-E) of figure 1 is used at the interface point to find G(B) if (from the relation with ), the length of sub-system B at complete equilibrium would be D x,eq (B * ) = D x (B 0 ) + D x (B 0 )∆G, where D x (B 0 ) is the total length of the sub-system at the initial time.The spatial average of G(B * ), defined as G(B * ) av can be easily computed.The quantity G(B * ) av is needed in the following relation to find the current total length of the sub-system at a particular time: The same change of length is applied with opposite sign on the other sub-system.The new dimensions D x,t (A * ) and D x,t (B * ) define also new constant grid step sizes, ∆ x (A * ) and ∆ x (B * ).The final operation is to re-mesh the values of G( * ) at the previous time step onto the new uniform spatial grid.
It is worth to mention that in the procedure outlined above here, converting the change of G to the change of the total length of the sub-system is a two steps process.The first step makes use of the relation between the change of G and the change of the total mass, which was illustrated in panel 1-D) of figure 1.In the next step the assumption is that the change of mass (and G) is proportional to the change of the total length of the sub-system.
To summarize the numerical procedure, at every time step the complete solution on both sides is obtained by solving equation 6 for G(A * ) and G(B * ) with the boundary conditions imposed for the limits of the whole system and preliminary values for the interface points.Then the interface points are updated using the polynomial function and equation 7. The total length is then rescaled to account for the mass transfer and the numerical grid size is updated.This procedure is iterated until the variation between two iterations becomes negligible (typically convergence is set by: where the labels # 1 and # 2 refer to two iterative steps).
Once convergence has been reached, the oxide abundance can be found easily using the Chebyshev polynomial parameterization in which each oxide is related to a function of G(A * ) or G(B * ) (e.g. for M gO see panel 1-B) and 1-C) of figure 1).
For convenience the composition is identified in wt% since the normalized oxides (*) represent the grams of the components with respect to a total mass of 100 grams.Finally, knowing temperature, pressure and the variation of the bulk oxides composition in space and time, a thermodynamic equilibrium calculation can be performed at every grid point using the program AlphaMELTS to determine the local mineralogical assemblage.Several 1-D numerical simulations have been carried out with initial proportion ranging from 1:1 to 100:1.Some results from a test case with proportion 1:1 are shown in figure 2. Initial total length on both side is set to D x (A 0 ) = D x (B 0 ) = 100 (arbitrary units), the initial spatial grid step is ∆d x (A 0 ) = ∆d x (B 0 ) = 1.Time step is set to 4 (arbitrary units) and S(A*)=S(B*)=1.
The initial bulk composition of the two assemblages, that separately are in complete thermodynamic equilibrium, is the same Note the increase of the length on the A side and decrease on the B side.Bulk oxides abundance is also computed at every grid point.The bulk M gO (wt%) is reported on panel 2-B), which shows the progressive decrease on the A side while MgO increases on the B side.The bulk composition can be used with the program AlphaMELTS to determine the local equilibrium assemblage.Panels 2-C) -2-H) show the amount of the various minerals in wt% (solid lines) and the M gO content in each mineral in wt% (dotted lines), with the exception of coesite in panel 2-H) (SiO 2 ).The complex mineralogical evolution during the chemical equilibration process can be studied in some detail.For example one can observe the progressive disappearance of orthopyroxene on the peridotite side and the exhaustion of coesite on the gabbro/eclogite side.
Similar results are shown in figure 3 and 4 for models with initial proportion set to 5:1 and 50:1, respectively.Differences in the numerical setup of the new test cases can be summarized as follow.For the 5:1 case: Few observations can be made by comparing the three simulations.For example, orthopyroxene on the peridotite side becomes more resilient and the total amount of Opx increases with the size of the initial sub-system.On the other side it appears that the M gO content in garnet (pyrope component) is greater for the model with starting proportion 5:1, compared to the 1:1 case.
However with initial proportion 50:1, the M gO content does not seem to change any further.
The supplementary material provides a link to access the raw data (all nine oxides) for the three test cases with initial proportion 1:1, 5:1 and 50:1.In addition two animations (1:1 and 5:1 cases) should help to visualize the evolution of the numerical models over time.

Application to the Evolution of a 2-D Model with One Dynamic Assemblage and Variable Extension
A 2-D numerical model makes possible to study cases in which at least one of the two assemblages becomes mobile.The simplest design explored in this section, considers a rectangular box with a vertical interface dividing the two sub-systems.
The dynamic condition is simply enforced in the model by assuming that one of the two assemblages moves downwards with a certain velocity, replaced by new material entering from the top side, while the other assemblage remains fixed in the initial spatial frame.The whole system evolves over time following the same principles introduced in the previous section.The numerical solution of the 2-D model is approached at every time step in two stages.In the first stage the following equation is applied to both sub-systems: where d x ( * ) is the general spacing in the x-direction representing either d x (A * ) or d x (B * ) and the vertical spacing d y is assumed to be the same on both sides.This equation is solved numerically using the alternating-direction implicit method (ADI) (Peaceman and Rachford, 1955;Douglas, Jr., 1955) which is unconditionally stable with a truncation error O(∆t 2 , ∆d 2 x , ∆d 2 y ) (Tannehill et al., 1997).Similar to implicit methods applied for 1-D problems, the ADI method requires only the solution of a tridiagonal matrix.
The numerical procedure described in section 3 to determine G( * ) at the interface is also applied here to the 2-D model.
The limits of the whole system opposite to the interface (left/right) are also treated similarly, assuming either a closedtype or symmetric-type boundary.For the other two boundaries (top,bottom) the zero flux condition is imposed, G(A * ) t,b l = G(A * ) y=1,ny and G(B * ) t,b l = G(B * ) y=1,ny where n y is the total number of grid points in the y direction (excluding the boundaries).
In the previous section a procedure was developed to account for the mass transfer between the two sub-systems.The same method is applied for the 2-D problem.The conceptual difference is that in a 2-D problem the mass change in principle should affect the area defined around a grid point.For practical purposes however in this study it only affects the length in the horizontal x-direction, hence re-meshing due to the change of mass is applied only to determine D x,t (A * ) and D x,t (B * ) and the two uniform grid step sizes in the x-direction, ∆d x (A * ) and ∆d x (B * ).
Up to this point the evolution of the system is not different than what was described for the 1-D case.The dynamic component is included at every time step in the second stage of the procedure.It is activated at a certain time assuming that the chosen sub-system moves downwards with a fixed pre-defined vertical velocity (y-component).The material introduced from the top side is assumed to have the same composition of the initial assemblage as defined for the 1-D models, table 1 (and the same G(A 0 ) and G(B 0 ) values).This is accomplished by assigning G(A 0 ) or G(B 0 ) at a location near the interface which is defined by the imposed velocity.Then the G( * ) points are also shifted according to the prescribed velocity.Values of G( * ) y on the original orthogonal grid are obtained by linear interpolation of the shifted G( * ) points.
Oxides bulk composition is then retrieved at each grid point over time using the same polynomial functions applied for the 1-D problem.The complete mineralogical assemblage can be also computed using AlphaMELTS as part of a post-process step after the numerical simulation is completed.
Only few 2-D simulations have been performed, specifically considering the initial proportion 1:1, 5:1 and 50:1, assuming either one of the two assemblages moving downward.Figure 5 summarizes some of the results for the case 5:1(A), i.e. with moving sub-system A. Initial grid specifications are: arbitrary units).Time step is set to 16 (arbitrary units).The scaling coefficients S x ( * ) and S y ( * ) are set to 0.01 (arbitrary units).The dynamic component is activated at time=100000 with vertical velocity set to 0.00625 (arbitrary units).The figure is a snapshot of the whole system soon after sub-system A has been activated downwards (time=102400).Panel 5-A) shows the variation of G( * ), while panel 5-B) illustrates the bulk M gO distribution (wt%).The other panels, 5-C) -5-H), present an overview of the mineralogical distribution (flood contour-type) and the M gO content in each mineral phase (line contour-type), with the exception of panel 5-H) for coesite (SiO 2 ).The panels clearly illustrate the variations introduced by the mobile sub-system A. On the other side there is apparently no immediate effect on the assemblage B, however the long term effect is significant and becomes visible in a later figure (figure 7).
Figure 6 provides a similar overview for the case assuming 5:1(B) with sub-system B moving downward.The same numerical conditions described for the previous case apply for this case as well.This figure, which shows only one time-frame soon after the sub-system is mobilized, does not appear to reveal new remarkable features.However advancing the simulation, a clear effect becomes more evident near the interface.In particular changes of the chemical and mineralogical properties moving away from the top entry side are quite significant.An animation related to figure 6 is best suited to illustrate this point.This movie file and another file for the animation related to figure 5 can be downloaded following the link provided in the supplementary material.The raw data files which include all nine oxides for both simulations are also available online.

Summary of the 1-D and 2-D Models Approaching Chemical Equilibration
Figure 7 summarizes the results of all the 1-D and 2-D numerical test models when the whole system approaches or is close to chemical equilibration.In the static scenario, exemplified by the 1-D models (solid lines), by increasing the initial size of sub-system A, the mineralogical and compositional variations tend to be smaller (see panels 7-C) -7-H) and enlarged view around the interface, panels 7-C2) -7-H2)).It is the expected behavior since any change is distributed over a larger space of the sub-system.The variations of the minerals abundance in assemblage B (gabbro/eclogite-type) instead remain quite independent of the initial size of sub-system A. However the abundance of the minerals not necessarily is the same found in the initial assemblage.In particular the amount of garnet, clinopyroxene and coesite is quite different from the amount of these minerals in the initial assemblage.This difference is rather unaffected by the initial proportion of the two assemblages, which has been varied from 1:1 (f=1) to 100:1 (f=100).
The composition of the minerals in assemblage A (e.g.M gO illustrated in panels 7-CC) -7-HH)) follows a pattern similar to the minerals abundance.As the size of the initial sub-system increases, M gO tends to approach the oxide amount in the initial composition.A different result is observed for the composition of the minerals in assemblage B. Regardless whether the mineral abundance changes or remains close to the initial amount, the oxide composition varies quite significantly and in most minerals the difference is larger when f is set to higher values.
When one of the sub-systems is allowed to move (2-D models), the general observation on the long run is that the dynamic sub-system tends to preserve the assemblage that enters in the model.In this study this assemblage is set to be equal to the initial assemblage.Note that the 2-D data plotted in figure 7 refer to an horizontal section of extracted points at the middle vertical distance D y /2.When sub-system A is mobile (dotted lines), the behavior of assemblage B is similar to the static case, with some minerals changing their initial abundance, garnet, clinopyroxene, coesite and in part spinel.In the reverse case, with B set as the dynamic sub-system, the mineralogical abundance of A differs from the initial assemblage (dashed lines).But unlike the static cases, no significant variations can be noted with the increase of the initial proportion.
In terms of minerals composition (e.g.M gO, panels 7-CC) -7-HH) in figure 7), the dynamic sub-system preserves the composition of the entering assemblage.The immobile assemblage instead, shows a compositional variation that is larger than any change observed for the static cases.This variation remains somehow still independent of the initial proportion of the two assemblages, at least with f = 1, 5, 50.
Complete data for the bulk composition, which includes all nine oxides, is available for three 1-D models and two 2-D simulations following the instructions in the supplementary material.

Conclusions
The main objective of this work was to show that a chemical heterogeneous mantle does not necessarily mean that different lithologies are in chemical disequilibrium (at least not entirely).
Often geochemical and petrological interpretations of the Earth interior rely on the achievement of thermodynamic equilibrium on a certain scale.The use of phase equilibrium data and partition coefficients, for example, does imply that chemical equilibrium has been achieved and it is maintained.Curiously, while this assumption is tacitly imposed on the most convenient dimension to interpret observed data, chemical equilibration is ignored when it comes to discuss the presence or the extent of chemical heterogeneities (i.e. chemical equilibration, in this regard, is considered ineffective) (e.g.Morgan, 2001;Ito and Mahoney, 2005a, b;Strake and Bourdon, 2009;Brown and Lesher, 2014).
Geophyisical interpretations usually require to specify certain properties, such as the density for the Earth materials under consideration.For example when the density is considered representative of real rock assemblages, the system has to be sufficiently small that the gravitational force is almost completely balanced by the pressure effect (viscous forces are ignored for simplicity), effectively establishing a quasi-static or static condition.Under this condition then, thermodynamic equilibrium can be achieved when the system is also equilibrated chemically, so that petrological constraints can be applied to determine the density of the assemblage.When different lithologies are considered in geophysical applications, it is assumed that chemical equilibrium is never achieved among them, regardless of the size of the system or the temporal scale.For studies whose conclusions are based on geological processes lasting for hundreds or billion of years, such assumption should be carefully verified considering that chemical and mass exchange are always effective to a certain extent.
The results from 43 study models (section 2.1) suggest that the imposed condition of thermodynamic equilibrium for the whole system (sum of two sub-systems) defines two new assemblages that are not homogenized compositionally or mineralogically, and their equilibrated compositions are different from those in the two initial assemblages.The two new assemblages not only define a condition of chemical equilibrium for the whole system but they also represent the equilibration within each separate sub-system.In addition, mass exchange between these equilibrated assemblages does not progress any further when the initial mass proportion of the two is varied and a new equilibration model is imposed to the newly defined whole system.
The results of the study models have been condensed in a series of parameterized functions that can be used for various applications (section 2.2).
A semi-empirical quantitative forward model was also developed to describe the evolution of the chemical equilibration process in the mantle.The model has been restricted to one set of values for the pressure and temperature and one pair of bulk compositions indicative of a peridotite-type and a gabbro/eclogite-type.The gabbro/eclogite-type can be interpreted as a portion of a subduction slab.Ignoring a thin sedimentary layer, that possibly could peel off during subduction, a large portion of the slab consists also of a depleted peridotite.Three lithologies (mantle peridotite, gabbro, depleted slab peridotite) probably can be also approached with a chemical equilibration model similar to the one presented here.However it remains to be seen whether the difference in composition with respect to the generic peridotite assumed in this study would lead to significant new results that would justify the additional modeling effort.
A priority was given here to understand the influence on the final assemblages of various initial proportions of the two subsystems and, to a limited extent, the effect of the initial compositions.The spatial and temporal evolution necessarily assumes arbitrary units.The reason behind it is that a comprehensive approach to study chemical heterogeneities that would include time-dependent experiments and suitable models for the interpretation of the experimental results has not been developed yet.
Experimental data are also necessary to validate certain assumptions that were made to model the composition of the two equilibrated assemblages (section 2).
The choice made to describe the variation of G( * ) using the transport model presented in section 3 and 4 may seem rather arbitrary.While details of the transition towards chemical equilibration should be investigated by experimental studies, the main point of the models in section 3 and 4 (and of this study) is to show that different lithologies can evolve while preserving distinct chemical and mineralogical features.The idea of using the concept of local Gibbs free energy variations over time and space (Kondepudi and Prigogine, 1998) to describe the chemical changes is a practical mean to simplify a problem that otherwise becomes intractable for complex systems.The choice is not a complete abstraction, it is approximately based on the consideration that the mass exchange is not governed by the compositional gradient but by the differences in the chemical potential of the various components in the various phases (e.g.Denbigh, 1971).Ultimately only extensive experimental studies could determine whether the simple evolution model for G(*) applied in this work to an heterogeneous system can be considered a reasonable approximation for describing the chemical evolution in practical geodynamic mantle models.Two aspects of the numerical applications presented in the previous sections deserve perhaps a further consideration.The assumption made for the composition of the entering assemblage in the 2-D models perhaps should be reconsidered in future studies.The other consideration concerns the boundary condition imposed on the opposite side of the interface between the two assemblages.The assumption is that the whole system is either close to mass exchange or mirror images exist outside the boundary limits.From a geological perspective the first scenario is probably the more difficult to realize.On the other hand the possibility that periodic repetitions of the same model structure are replicated over a large portion of the mantle, if not the entire mantle, seems more reasonable.Assuming that the time scale is somehow constrained, an investigation of the temporal evolution would still require some kind of assessment of the periodic distribution of the thermodynamic system as a whole.
The 2-D simulations in which one of the assemblages is allowed to move, have shown that on the long run the mineralogical abundance and compositional variations are approximately independent of the size of the two sub-systems.This observations suggests the possibility of implementing large geodynamic models with evolving petrological systems, once the temporal and spatial scale of the chemical changes have been constrained.
At the moment the spatial and temporal variations are arbitrarily defined, but this study shows that the petrological and mineralogical changes may still be approximately quantified, at least at the (P,T) conditions that have been considered.It would be useful for example to select few bulk compositions for the two sub-systems and apply them to the dynamic equilibrium melting (DEM) and dynamic fractional melting (DFM) models that have been developed combining 1-D multiphase flow with AlphaMELTS (Tirone and Sessing, 2017;Tirone, 2018).Perhaps even a simplified model for non-equilibrium fractional crystallization could be applied to try to reproduce observed 3-D chemical zoning in minerals and multicomponent chemical zoning in melts (Tirone et al., 2016).More in general the results should be compared with existing data on melt products and residual solids observed in various geological settings to investigate indirectly, but from a quantitative perspective, the presence of chemical heterogeneities in the mantle.It becomes also possible to determine the variation of physical properties, such as bulk density, and relate them to certain observables, such as seismic velocities.At least on a relative scale, the effect of the compositional variations could be associated to seismic velocity variations, providing in this way another indirect evidence of heterogeneities in the mantle based on a quantitative forward description.In this model it is assumed that at time 100000 a new assemblage B enters from the top with velocity 0.00625 (arbitrary units).The description of the panels follows the caption of figure 5.     4).The initial proportion of the whole system is f:1 (f=5).The description of the results follow the outline of the caption of table 3. ------------------mol --------------   ------------------mol -------------- . Consider for example the garnet component almandine (Alm) which appears in reaction (T-1), (T-3), (T-10), (T-12), (T-13), (T-14), (T-15) and (T-16), the following relation can be established: 7) to retrieve G(A * ) if and G(B * ) if assuming that S(A * ) = S(B * ).The external boundaries defining the limits of the whole system (symbol l) are assumed to be of closed-type or symmetric-type.Both are obtained by the condition G(A * ) l = G(A * ) nA−1 and G(B * ) l = G(B * ) nB −1 , where n A and n B are the total number of grid points on each side (excluding the boundary points).

Figure 1 .Figure 2 .Figure 3 .Figure 4 .
Figure 1.Data and relative fitting of 43 study cases that are used to develop the chemical equilibration model.Panel 1-A) relation between the ratio G(A * )/G(B * ) and G(B * ) which is applied to constrain G(A * ) and G(B * ) at the interface.Panel 1-B) and 1-C) illustrate the relation between G(A * ) and G(B * ) with M gO bulk abundance.Similar relations are applied for all nine oxides defining the bulk composition.The normalized bulk abundance is intended as grams with respect to a total mass of 100 grams which is equivalent to wt%.Knowing G(B), the total size of the assemblage at equilibrium can be found assuming that a) a relation between the mass change and the change of G(B) is established (Panel 1-D), b) the extension of the assemblage is proportional to the mass change and it takes place along a direction perpendicular to the interface.The total length at equilibrium is then adjusted in accordance with the difference between the spatial average G(B * ) of the assemblage and G(B * ) at the interface (see the main text for a detailed explanation).The change of size of the second assemblage is also applied on the first assemblage but with opposite sign.Panel 1-E) allows to determine G(B) from the relation with G(B * ) at the interface.

Figure 5 .
Figure 5. Solution of a 2-D model simulation at time 102400 (arbitrary units).The starting proportion of the two assemblages is 5:1 (f = 5).In the initial setup the 2 assemblages are separately in chemical equilibrium.At time 100000 a new assemblage A enters from the top side with velocity 0.00625 (arbitrary units).The new assemblage is assumed to have been equilibrated but never previously in contact with assemblage B (the composition of the new assemblage is the same of the assemblage in the initial setup).Panel 5-A) spatial variation of G( * ).Panel 5-B) local distribution of MgO in the bulk assemblage.Similar results are obtained for all the other oxides defining the bulk composition.An animation file and raw data for all nine oxides are available online following the instructions provided in the supplementary material.Panels 5-C) -G) local minerals distribution (color map) and few contour lines for the abundance of M gO in the associate minerals.Panel 5-H) spatial distribution of coesite.Time and distance in arbitrary units.Pressure and temperature are fixed at 40 kbar and 1200 o C.The rest of the parameters for the numerical model are defined in the main text.

Figure 7 .
Figure 7. Summary of the results for all the 1-D and 2-D numerical models at conditions close to chemical equilibrium for the whole system.The models consider different initial proportions of the two assemblages.In addition for the 2-D models it is assumed that either assemblage A or B enters from the top side at time 100000 (arbitrary units) with velocity 0.00625 (arbitrary units).For the 2-D models the profiles represent an horizontal section at the vertical middle point (Dy/2).Panel 7-A) spatial variation of G( * ).For clarity, plot of the 2-D model with 50:1(B) is truncated at x ∼ 500.Panel 7-A2) enlarged view of G( * ) near the interface.Panel 7-B) variation of bulk M gO (wt%).Panel 7-B2) enlarged view of bulk M gO near the interface.Panels 7-C) -G) spatial variation of minerals abundance.Panels 7-C2) -G2) minerals abundance zoomed near the interface.Panels 7-CC) -GG) M gO content in the associated minerals near the interface.Panels 7-H) and 7-H2) distribution of coesite (SiO2).

Opx
).The second and third column on the upper side of the table report the input bulk composition on the two sides.The second and fifth column on the lower part of the table show the results of the thermodynamic equilibrium calculation applied separately to the two sub-systems.The last column shows the results for the whole system W .This last column indicates for example that orthopyroxene is not present at equilibrium in the whole assemblage.Considering the reactions in table 2 and the data in table3, the En component in orthopyroxene appears only in reaction T-2, and since no OEn is present on the B side, the mole change in A can be locked (∆n OEn (A) = −0.0700777).Therefore ξ (T −2) is fixed to -0.0700777.The same is also true for ξ (T −3)

Table 1 .
List of minerals and mineral components relevant for this study with chemical formulas and abbreviations.

Table 2 .
Set of independent reactions for the list of mineral components in table 1.

Table 4 .
Normalized bulk composition (A * ) and (B * ) in the two sub-systems taken from the results of the model in table 3, (A) and (B).The lower part of the table shows the equilibrium mineral composition computed with the program AlphaMELTS for each sub-system separately.

Table 5 .
Summary of the results of a chemical equilibration procedure in which the initial composition of the two-sub-systems (A0) and (B0) is taken from the outcome of the previous model (A * and B * from table

Table 6 .
Results from a chemical equilibration model with initial composition of the two sub-systems (A0) and (B0) analogous to the one presented in table 3. The only difference is that the initial proportion of the whole system is f:1 (f=5).

Table 7 .
Normalized bulk composition (A * ) and (B * ) of the two sub-systems taken from the results of the model in table 6.The lower part of the table shows the equilibrium mineral composition computed with the program AlphaMELTS for each sub-system separately.