Exsolution by spinodal decomposition in multicomponent mineral solutions.

Phase separation in a three-component system that results from the uphill diffusion of chemical components is considered. The binary decomposition model of Cahn and Hilliard is generalized to account for the interdiffusion of several chemical components with considerably different diffusion constants. Thereafter the decomposition dynamics and the phase relations of the final system state are investigated by means of finite-element modeling. Examples from a hypothetical regular solution and from ternary feldspar are addressed. Special attention is given to situations in which different diffusivities affect decomposition dynamics and the final system states. Good qualitative agreement between our modeling and petrographic observations on exsolved feldspar is achieved. Our model explains systematic deviations from equilibrium element partitioning between the two phases exsolving from an initially homogeneous ternary feldspar during slow cooling.

Introduction

Uphill diffusion of chemical components and the resulting spinodal decomposition of phases in alloys, polymers and minerals are usually associated with the interdiffusion of atoms. Many mineral systems may be described as ionic crystals and interdiffusion occurs by a vacancy mechanism. Provided that the concentration of vacancies is constant over time, such interdiffusion implies that the fluxes of the substitutional components in a binary solid solution take identical values and have opposite directions in the material fixed frame of reference. Therefore both moving components can be characterized by a single mobility [1] or by a single diffusivity [2] parameter. Following these lines, spinodal decomposition in a multi-component system is usually considered under the assumption of one common mobility (diffusivity) constant [3,4,2,5,6].

The proposition of nearly equal diffusivities of components in a multicomponent solution may be a reasonable approximation for spinodal decomposition in alloys and in polymers, but it definitely fails for mineral solid solutions. The latter show diffusivity parameters differing up to several orders of magnitude [7]. Investigations of phase separation in multicomponent systems with different diffusivities are relatively scarce [8,9]. These studies typically conclude that neither the newborn phases nor the morphology of the final system state is affected by possibly different mobilities of the components. This is a natural conclusion because the stationary equilibrium state, say for given pressure and temperature, is determined by equilibration of chemical potentials. As there are no fluxes, the mobilities cannot affect the uniform phases with the exception of the transition regions between the different phases. The transition regions are coined by the interplay of particle mobilities and interfacial energy [10,11].

Mineral solid solutions that are originally formed at high temperatures that are typical of the Earth’s crust and mantle are often observed far from the equilibrium state. Such transient states may get “stranded” due to successively more sluggish kinetics during cooling. Therefore the whole evolution path may be of interest. This path is of course affected by the contrasting diffusivities of the components. Spinodal decomposition in mineral solid solutions is the main topic of the present study. We first revisit the standard derivation of the Cahn–Hilliard equation (see e.g. the review by Nauman and He [2]) and carefully account for the arbitrary diffusivities of the components. The key point is that the flux j of a component A is not proportional solely to −∇μ but results from additive contributions of gradients of all chemical potentials in the spirit of non-equilibrium thermodynamics [12]. Substitutional interdiffusion is ensured by special properties of the mobility matrix. In this paper we are primarily interested in new physical effects caused by the considerably different mobilities of the components. To illustrate these new phenomena as simply and cleanly as possible, we deliberately ignore all secondary effects, e.g. contributions of the induced stresses, possible changes in the number of vacancies and possible influence of the induced electric fields on diffusion of components in ionic crystals. After the multicomponent generalization of the Cahn–Hilliard equation is derived, we consider how the final equilibrium or a long-lived quasi-stationary state is approached by the system. To this end, the generalized Cahn–Hilliard equation is solved numerically in several typical situations for a hypothetical regular solution and for ternary feldspar. Special attention is given to considerably different diffusivities and to the intermediate long-living quasi-stationary system states that may differ from those predicted by equilibrium thermodynamics. Finally, our modeling results are discussed in the light of observations on naturally exsolved feldspar from slowly cooled rocks.

Derivation

Interdiffusion

In this paper we consider spinodal decomposition in a multicomponent system. The latter is characterized by specifying its uniform temperature T = const and an externally applied uniform pressure P = const. We presuppose that atoms do not leave the system, occupy identical volumes and only rearrange within some fixed crystalline sublattice. The remainder of the crystal structure remains unchanged, local stresses are not created and the system volume V is a constant. The difference between the Helmholtz free energy and the Gibbs free energy is then also a constant, and both quantities are equally suitable for describing the system. In what follows, the free energy refers to the Gibbs free energy. Moreover, there is no distinction between the laboratory reference frame and the lattice reference frame because both frames are associated with the same fixed crystalline sublattice.

The local concentrations of atoms per unit volume are denoted by n(r, t), index α = A, B, C, … enumerates components. We consider substitutional interdiffusion, which is a special case of interdiffusion for which the total volume concentration of atoms does not change,In particular, the total flux of atoms is zero in a suitable frame of reference. To quantify the rearrangement dynamics, we have to introduce an expression for the local flux of each component j(r, t) such thatin accord with Eq. (1). The flux model will be derived by incorporating Onsager’s approach [12].

Flux model

Rearrangement of the atoms pushes the system towards an equilibrium state with the smallest possible free energy. The corresponding thermodynamic forces act on atoms and are opposite to the spatial gradients of the chemical potentials. The resulting atomic fluxes are given by the Onsager-type relationin which μ = μ(r, t) denotes the chemical potential of a component, −∇μ is the thermodynamic force per atom and L is the kinetic coefficient, also referred to as mobility. One naturally expects thatIn other words, the driving force −∇μ in Eq. (3) yields a flux of component β in the same direction as the force and fluxes of all the other components α ≠ β in the opposite direction. The matrix L is assumed symmetricin analogy with the reciprocal relation [12].

To address interdiffusion, we impose an additional restrictionNow the total flux (2) vanishes automatically because Eqs. (3) and (7) implyirrespective of the specific values of ∇μ.

The symmetry relation (6) and the interdiffusion condition (7) imply that Eq. (3) is equivalent toand the particle flux actually depends only on the difference of chemical potentials. In what follows we will see that, although the expressions for the chemical potentials may be quite complicated, their differences are considerably simpler. It is therefore more convenient to deal with Eq. (8) instead of the parent Eq. (3).

Another consequence of condition (7) is that only the off-diagonal elements of the matrix L need to be specified. For instance, for a binary system with α = A, B, we have an especially simple formin which only one parameter L < 0 counts. The fluxes are given by expression (8)accounting for the fact that the positive factor −L relates the flux j and the corresponding driving force −∇(μ − μ) for interdiffusion of unlike atoms moving in opposite directions. Cahn and Hilliard [1] considered Eq. (9) as a phenomenological definition of the diffusive mobility M = −L > 0.

For a three-component system, all nine elements of matrix L are determined by only three independent negative parameters:Particle fluxes are given by the relationsand the total flux automatically vanishes.

Diffusivities

Interdiffusion of unlike atoms is only possible if both types of atoms are present in a given volume element. For example, in a two-component system the mobility L in Eq. (9) must be concentration dependent and must vanish if either n = 0 or n = 0. It is then useful to introduce a positive coefficient D such thatwhere k is Boltzmann’s constant and is defined by Eq. (1). The parameter D is referred to as the interdiffusion coefficient (this notation is justified later on in Section 3.1). Here we stress that the mobility coefficient L is composition dependent whereas the interdiffusion coefficient D may be constant, as discussed in Ref. [2]. This is why it is more convenient to operate with D.

Following these lines, we return to the multicomponent case and define a set of positive parameters D such thatThe diagonal values L are uniquely deduced from Eq. (7). The symmetry relation (6) immediately indicates that D = D. Using definition (11), the flux model (8) can be rewritten aswhere α, β = A, B, C, … enumerate all components. Eq. (12) specifies our flux model for a multicomponent interdiffusion system.

Chemical potential

To apply Eq. (12), we require explicit expressions for the chemical potentials. These are derived from the system free energy , the latter being defined as a spatial integral of the volume density Gwhere n is a short notation for all particle densities and a dependence of G on ∇n is introduced in the spirit of the seminal work [1]. The integration in Eq. (13) is performed over the system volume V. In principle, the latter is a self-organized quantity determined by P, T and the total number of atoms of each component. However, V = const in our interdiffusion framework.

To determine chemical potentials, one must consider a small variation of the particle concentrations n → n + δn. Note that the integration volume does not change, and contributions to the particle number due to fluxes through the system boundaries are ignored. In other words, δn represents some internal rearrangement of atoms (and in general their chemical interconversion). The corresponding variation of the system free energy can be transformed to the formwhere partial integration was performed to eliminate terms proportional to ∇δn anddenotes a variational derivative. System-interface related contributions to are set to zero. Now we identify the variational derivatives with the chemical potential of the componentssuch that variation of the system free energy in a unit volume takes the well-known form, .

To compute μ, all components n are considered as independent quantities. In other words, definition (15) is independent of the interdiffusion condition (1). It is the relation (7) for the mobility matrix that accounts for interdiffusion.

Dynamic model

We can now specify a dynamic equation for decomposition of phases in a multicomponent system. For a given free energy density G(T, n, ∇n), the chemical potentials are first calculated using Eq. (15). Second, they are inserted into Eq. (12) to derive fluxes. The latter are then inserted into a sourceless continuity equation for the particle densityEq. (16) is finally solved numerically to obtain n(r, t).

Our dynamic model can be written asand formally contains one equation for each component. However, these equations are not independent due to the interdiffusion condition (1), which is automatically supported by all solutions of Eq. (17). Specifically, one can always eliminate one component by reducing Eq. (17) to the equations for mole fractions such that, for example, a two-component system is described by a single equation.

Some useful system properties can be derived prior to numerical solution of Eq. (17). For instance, let us trace the evolution of the free energy . In accord with Eqs. (14) and (16), we obtainwhere partial integration was performed in a second step. For the flux model (3), we immediately obtainFurthermore, using Eq. (7), we transform Eq. (18) to the formand conclude that permanently decreases because of inequality (5). In particular, stable equilibrium system states can be calculated by constrained minimization of the free energy. The constraints are derived from the condition that is conserved for each component. Specific examples are given in Appendices A and B.

Model equations

Two components

A simple illustration of the model (17) is given by a two-component system with α = A, B. We define the mole fraction of the component Aand introduce the following expression for the free energy density [1]Here g(a) is just a free energy per particle in a uniform system. It accounts for both the configurational and the non-ideal contributions to the free energy. The gradient term with the parameter κ > 0 describes the interfacial energy associated with the phase boundaries. Using Eq. (15), one directly obtainswhere as expected, G = μn + μn. Although both chemical potentials look quite complicated, their difference is much more simple:so that, for α = A, Eq. (17) takes the formwhich is the well-known nonlinear Cahn–Hilliard equation for an isotropic binary system (for a recent review see Ref. [2] and references cited therein; anisotropic models for the surface free energy term are discussed in Refs. [13-17]). The dynamic equation for the second component, α = B, trivially follows from Eq. (21) and can be dropped. The quantityis referred to as an effective diffusion coefficient. If the κ term can be ignored, Eq. (21) is just a nonlinear diffusion equation:Moreover, for an ideal mixture withwe obtain Deff(a) = D = const, which is why parameter D is associated with the interdiffusion coefficient [2].

More complex behavior is possible in a non-ideal mixture where an unstable (spinodal) region of concentrations appears. This region is determined by the inequalityand physically corresponds to the uphill diffusion of components. A typical example is given by a non-ideal mixture in whichand the uphill diffusion and spinodal phase separation are possible if the parameter χ > 2. When uphill diffusion is present, the κ term in Eq. (21) cannot be ignored; this term finally stabilizes the system.

Three components

We now formulate a closed system of two equations for a three-component system with α = A, B, C. The system is completely described by two mole fractions,where . In full analogy with Eq. (19), the free energy density is presented aswhere the specific expression for f(a, b) is determined by the free energy per particle for the uniform system states. (In full analogy with g(a) from the previous section, a different notation is used to avoid confusion between two- and three-component systems.) The most simple f(a, b) is given by the expressionin which c = 1 − a − b. More complicated models for f(a, b) in real materials will be considered in the numerical section.

The gradient terms in Eq. (24) should give a positive contribution for arbitrary directions of ∇a and ∇b. This implies the restrictionsIn what follows, we use the short notationsThe chemical potentials of the components are calculated in accord with definition (15) and readwhere, as expected,One sees that the differences between chemical potentials are given byThe governing equations for the field variables (23) are obtained by inserting the latter differences into Eq. (17):where c = 1 − a − b. Similar to condition (22), the uphill diffusion of components is possible if either ∂2f/∂a2 < 0 or ∂2f/∂b2 < 0. More accurate conditions for spinodal decomposition are derived in Appendix B.

System (27) is the starting point for describing phase decomposition in a three-component system. It can be solved numerically after an expression for the free energy (24) and diffusion coefficients D are specified. An equation for the third component is derived in a similar manner,but does not provide new information.

Special cases

In this section we describe several special cases for the general system (27). First of all, a connection between the three- and two-component systems can be found, for example, by setting b → 0. Thereafter, the second equation in Eq. (27) is satisfied automatically, whereas the first one is simplified to the formwhich is equivalent to Eq. (21).

Another special case of Eq. (27) is that, for identical diffusion coefficients, i.e. replacing D, D and D with one common D, we obtainSpinodal decomposition in such systems was examined in Ref. [4], and possible morphologies are listed in Ref. [3].

In this paper we are primarily interested in systems with considerably different mobilities of components. For instance, in minerals, as opposed to alloys, such a difference may be several orders of magnitude. An extreme situation occurs when e.g. D = 0, so that direct interdiffusion of the components A and B is not possible. Both components are transported due to binary exchanges with the third component. The equations in Eq. (27) are then reduced to the systemInformally, one may say that Eq. (29) represents an effective two-component system in which in general D ≠ D.

Equilibrium states

If the chemical potentials of all components are constant throughout the system, the fluxes predicted by Eq. (12) vanish, and the system is in equilibrium. The equilibrium state does not depend on diffusion coefficients, the latter only determining the path to a final state. The main properties of the equilibrium states for two- and three-component systems are described in Appendices A and B. Here we stress that sometimes it is important to account for the differences in the diffusion coefficients even near to the equilibrium. This happens, for example, when the system is cooled down, such that some intermediate non-equilibrium system state becomes “frozen”. It is then critically important to reconstruct the path to which the intermediate state belongs. As we will see later, considerably different diffusion coefficients may force the system to evolve along a “wrong” conode, returning to a “correct” one only in close vicinity to the final state. Such situations will be considered later on for ternary feldspar.

Modeling

In this section we provide exemplary numerical solutions of the basic system (27). All solutions are obtained with the Multiphysics Modeling and Simulation Software, COMSOL [18].

Decomposition into two phases

As a first example, we consider a non-ideal three-component mixture where the free energy is given by the expression (25), in which χ = 2.6, χ = 2.5 and χ = 2.7. To begin with, we calculated the corresponding phase diagram. Barycentric coordinates are used, such that each point within the Gibbs triangle in Fig. 1 simultaneously represents a, b and c. Red regions indicate stable and metastable uniform system states. These regions are bounded by spinodal lines, which are calculated using Eq. (B.5). Straight black lines (conodes) indicate the directions along which the decomposition tends to occur; these lines are bounded by the binodal points. More precisely, if the initial (uniform) system state belongs to a conode, the final (uniform) phases are determined by the binodal points at the ends of the same conode. However, spinodal decomposition leads to highly non-uniform intermediate states, and compositions of the newborn phases may then temporarily leave the conode and finally return to it. Examples are given in Sections 5.3 and 5.4.

Fig. 1

Phase diagram for a three-component system with the free energy (25) and χ = 2.6, χ = 2.5 and χ = 2.7. The initial system states for two exemplary calculations of phase decomposition are indicated by stars and the final decomposed phases are indicated by points.

Returning to Fig. 1, we observe three smooth binodal curves within the red regions. These curves separate the stable and metastable system states. To construct the binodal curves, we numerically build a convex hull over the free energy surface. The line at which the convex hull intersects the free energy surface is the binodal curve. The facets of “tangent planes” give the conodes.

Two exemplary initial concentrations are chosen, as shown by the black star (a0, b0) = (0.25, 0.5) and blue star (a0, b0) = (0.5, 0.1) in Fig. 1. The finally decomposed phases are shown by the black and blue points on the binodal curves on the same phase diagram. As expected, the initial stars and the final binodal points belong to the same conode.

The distributions a = a0 and b = b0 with some small additive random perturbation are used to trigger decomposition. Periodic boundary conditions are applied. After the quick initial separation has taken place, the system evolves into a nearly stationary state. Exemplary nearly stationary spatial distributions of the mole fractions are shown in Fig. 2. Two phases are clearly observed. The corresponding values of the mole fractions agree with those taken from the phase diagram. The further evolution of such a state is extremely slow.

Fig. 2

Nearly stationary 1-D distributions of mole fractions for Eq. (27). Free energy parameters and initial concentrations correspond to those in Fig. 1. Top: (a0, b0) is given by the black star in Fig. 1. In the final state, a(x) ≈ c(x); therefore only a(x) and b(x) are shown. Bottom: (a0, b0) is given by the blue star in Fig. 1. Here all three mole fractions are different.

In a next step we consider decomposition of phases in the same system as above, but in two spatial dimensions. Exemplary snapshots of the density distribution for a(x, y, t) are shown in Fig. 3. The initial condition is given by the blue star in Fig. 1. The snapshots demonstrate all typical features of spinodal decomposition: a quick initial separation and a much slower further coarsening. The morphology is richer than in the one-dimensional (1-D) case, e.g. labyrinths are observed for our initial condition.

Fig. 3

Snapshots of a(x, y, t) for numerical solution of Eq. (27) in two spatial dimensions. The initial composition is given by the blue star in Fig. 1. The snapshots are taken at times 500, 1500, 3000, 5000, 10,000 and 18,000. The simulation domain is 80 × 80.

Decomposition into three phases

In this section, spinodal decomposition into three phases is considered. The system free energy is again described by Eq. (25), in which χ = 2.3, χ = 3.0 and χ = 2.4. The corresponding phase diagram is calculated as in the previous section. The diagram is shown in Fig. 4 and has a new feature as compared to Fig. 1, namely in the white (unstable) domain there is a small region free of conodes. If the initial concentration lies in this region, spinodal decomposition into three phases will take place. Let us consider one such example.

Fig. 4

Phase diagram for a three-component system with the free energy (25) and χ = 2.3, χ = 3.0 and χ = 2.4. The initial system state for an exemplary calculation of phase decomposition is indicated by the star, the final decomposed phases by three points.

Eq. (27) has been solved twice: in one and two spatial dimensions. The initial composition (a0, b0) = (0.27, 0.3) is shown by a blue star in Fig. 4. A final quasi-stationary state for one spatial dimension is shown in Fig. 5. In the latter, one can clearly distinguish three phases. The final mole fractions (a, b), i = 1, 2, 3 agree well with those obtained from the phase diagram.

Fig. 5

Nearly stationary spatial distributions of the mole fractions a(x), b(x) and c(x) for the three-component system (27) in one spatial dimension. The initial uniform concentrations correspond to system decomposition into three phases, as shown by the blue star in Fig. 4. The free energy parameters read χ = 2.3, χ = 3.0 and χ = 2.4.

Exemplary snapshots of the mole fraction distribution a(x, y, t) for the numerical solution of Eq. (27) in two spatial dimensions are shown in Fig. 6. Typical features of the spinodal regime of phase separation are clearly observable: at first, several individuals of the new phases appear. These grow in size and reduce in number, and the system develops into a quasi-equilibrium state. The subsequent evolution is extremely slow. Three different phases are now involved in the decomposition.

Fig. 6

Snapshots of a(x, y, t) for a numerical solution of Eq. (27) in two spatial dimensions. The initial composition is given by the blue star in Fig. 4. The snapshots are taken at times 7000, 16000, 48000, 88000, 350000 and 8e5. The simulation domain is 80 × 80.

Ternary feldspar

We now turn to a real mineralogical example and consider the spinodal decomposition of ternary feldspar. The latter is a solid solution consisting of the three components, Or, An and Ab (orthoclase, anorthite and albite), with compositions KAlSi3O8, CaAl2Si2O8 and NaAlSi3O8, respectively. In the following their mole fractions are denoted by a, b, and c such that (a, b, c) ↔ (XOr, XAn, XAb). The feldspar structure consists of a 3-D framework of interlinked [SiO4]4− and [AlO4]5− tetrahedra. The Na+, K+, and Ca2+ ions occupy large irregular cavities in the tetrahedral framework. It is important to note that exchange of K+ and Na+ ions is considerably easier than exchange of either K+ or Na+ ions with Ca2+ because, due to charge balance constraints, such an exchange requires simultaneous transport of Al3+ and Si4+ in the tetrahedrally coordinated sites, which is exceedingly difficult. Therefore, our system consists of two mobile components (orthoclase and albite) and one less mobile component (anorthite). The dynamic model (27) accounts for different diffusivities and is thus suited for tracing decomposition dynamics. We use a realistic free energy function for the uniform states [19]To specify the coefficients in the last equation, it is convenient to introduce a new index,Each coefficient w★ is related to the three Margules parameters, W, W and W (excess internal energy, excess entropy, and excess volume parameters respectively), by virtue of equationwhere R is the gas constant. All involved Margules parameters are listed in Table 1 in accord with Ref. [19]. Exemplary phase diagrams of ternary feldspar for two temperatures are shown in Fig. 7. The straight lines indicate conodes joining points on the binodal curves that represent compositions of feldspars that coexist stably at the temperatures and pressures for which the diagrams have been constructed. The red dotted areas indicate the stable (outside the binodal regions) and metastable (inside the binodal regions) domains. The white area represents the unstable region. The boundary between the metastable and unstable regions is calculated from condition (B.5). Note that relatively orthoclase-poor feldspars with compositions close to the albite–anorthite join coexist with relatively anorthite-poor feldspars with compositions close to the albite–orthoclase join. The former pertain to the plagioclase series, the latter constitute the alkali feldspars. As expected, the difference in the compositions of the coexisting feldspars increases with decreasing temperature (compare the top and bottom diagrams in Fig. 7). In the following we consider the decomposition of ternary feldspar in a closed system, i.e. we ignore any interaction with a fluid or melt that could change the bulk composition. Given this constraint, it is noted that when going from high to low temperature alkali feldspar becomes successively more orthoclas-rich and anorthite-poor, whereas the plagioclase becomes successively more orthoclase-poor and anorthite-rich.

Table 1

Margules parameters for feldspar [19].

★	WE (J/mol)	WS (J/(mol K))	WV (J/(mol bar))
AbOr	18810	10.3	0.394
OrAb	27320	10.3	0.394
AbAn	28226	0	0
AnAb	8471	0	0
AnOr	52468	0	−0.12
OrAn	47396	0	0
OrAbAn	8700	0	−1.094

Fig. 7

Exemplary phase diagrams of ternary feldspar for two temperatures. The binodal curves separate stable and metastable initial states (both red). The unstable (spinodal) white region is covered by conodes, which also cover metastable regions and are bounded by binodal points. The spinodal curves are calculated in accord with Eq. (B.5). The free energy model (30) follows Ref. [19].

We now consider the decomposition dynamics, taking into account that D ≃ D ≪ D. As suggested by our numerical solutions, the flux of the b-component is initially ignorable. The main results are as follows. The system first behaves like a two-component system in which only one diffusivity D is present. That is, the newborn phases evolve not along the conodes shown in Fig. 7 but along lines that are approximately parallel to the base of the Gibbs triangle, such that b ≈ const. Only when the separation of the albite and orthoclase components is nearly completed does exchange of the anorthite component become effective and the system slowly returns to a correct conode.

This effect is illustrated in Fig. 8. The initial condition corresponds to T = 973 K, P = 8000 bar and mole fractions (a0, b0) = (0.4, 0.07), as indicated by a green star. If all three diffusivities are similar, the system predictably evolves along the conode (Fig. 8, top). In contrast, the above-described two-stage decomposition regime is observed for realistic relations among the diffusivities (Fig. 8, bottom).

Fig. 8

Phase decomposition of a ternary feldspar shown using paths on the Gibbs triangles. Top: diffusivities of components are artificially taken to be identical. Bottom: realistic diffusivities are taken. In the second case, the system first leaves the conode and then returns to it.

To trace the evolution of the intermediate highly non-uniform system states using the Gibbs triangle, one must permanently evaluate the instant mole fractions of the phases. For a given t, these mole fractions are obtained from a histogram indicating, for example, the probability of different values of a(t,x,y). An example is shown in Fig. 9. We observe two pronounced maxima, which are attributed to the values of a for the developing phases. The mole fractions are determined through a normal fit and are used to plot the evolution paths in Fig. 8.

Fig. 9

An exemplary histogram of different values of a(t, x, y) observed throughout the decomposing system for a fixed t. The two maxima are associated with the compositions of the newborn phases. These compositions are used to plot the paths in Fig. 8.

Cooled feldspar

Another interesting phenomenon resulting from the different diffusion coefficients occurs in a slowly cooled feldspar. The point is that diffusivities are very sensitive to temperature, e.g. nearly identical diffusivities of components for high temperatures may be replaced by considerably different ones in the course of cooling. Thereby the decomposition dynamics changes considerably. In particular, even the final quasi-stationary system state may lay off the binodal curve.

To illustrate the effect of cooling numerically, we have chosen a ternary feldspar at P = 8000 bar. The initial mole fractions are (a0,b0) = (0.4, 0.2). Cooling is performed in two steps. First, the ternary system with nearly equal diffusivities of the components decomposes at T = 1273 K until two nearly stationary equilibrium phases are formed. These phases are shown by the red points in Fig. 10 (top). Note that the initial decomposition occurs along the conodes.

Fig. 10

Different diffusion effects on cooling feldspar.

In a second step, the system is quickly cooled down to the temperature T = 973 K. The formerly created phases (red stars in Fig. 10, bottom) are used as a new initial condition. The less mobile anorthite component is now assumed to be “frozen”, unlike the others. Therefore D and D are set to zero, and further system evolution is determined by D only.

The system decomposes further, the final phases being shown by the green points in Fig. 10 (bottom). As b = const for the second stage, the compositions of both phases evolve along lines which are parallel to the base of the Gibbs triangle. Notably, the final compositions do not reach the binodal curves at all. The reason for this is that in the second stage we effectively have a two-component system, which does not behave in accord with the ternary phase diagram.

For comparison, we also calculated decomposition of the system at T = 973 K, taking identical diffusivities for all three components. In this case, the final phases (red points in Fig. 10, bottom) are different and predictably correspond to the conodes on the Gibbs triangle. We see that accounting for different diffusivities may lead to system states that deviate considerably from what is predicted for thermodynamic equilibrium.

It is a common phenomenon that the compositions of the two feldspars exsolved from a homogeneous ternary feldspar precursor do not lay on a binodal curve. This has been attributed to the ease of diffusion of the mono- and divalent cations in the feldspar differing. It has been argued [20,21] that prolonged exchange of Na+/K+ at cooling stages when Na+/Ca2+ and K+/Ca2+ interdiffusion has already essentially ceased leads to further separation of Na+ and K+, whereas the Ca2+ contents remain constant. In such a case, the compositions of the two feldspars produced from exsolution are shifted away from the common binodal curve and, if applied in two feldspar thermometry, give systematic errors. Voll et al. [20,21] suggested that this effect could be corrected for by drawing lines parallel to the albite–orthoclase join of the Gibbs triangle through the observed compositions and then identifying the equilibration temperature by finding the binodal curve for which both phases lay on the same conode. Although slightly inaccurate, this procedure in essence reflects the effect that we obtain from our numerical simulation. This simplistic correction cannot, however, account for the complex microstructural and compositional evolution during the decomposition of a homogeneous ternary feldspar. The information from both the observable microstructures and the mineral compositions can be fully exploited based on our kinetic model. An example of such an analysis was given for exsolution of a binary alkali feldspar [22]. Such a quantitative analysis for ternary feldspar is hampered because the mobility of the anorthite component is not yet well enough constrained. The knowledge of relative component mobilities in ternary feldspar could, in principle, be better constrained from confronting mineral chemical data from exsolved feldspars with our numerical simulations, but this is beyond the purpose of this communication.

Conclusions

In conclusion, we investigated spinodal decomposition in a three-component system with considerably different diffusion coefficients. After deriving a suitable generalization of the Cahn–Hilliard equation, we addressed phase separation for both an idealized regular solution and ternary feldspar. In principle, the system in question must finally reach the thermodynamically stable equilibrium state, in which there are no fluxes and the differences in diffusivities are unimportant. In actual fact, the system may become frozen in a long-lived quasi-equilibrium state, which is strongly affected by the specific values of the diffusivities. This is especially true when the system is cooled down. In particular, the system evolution may occur along paths that deviate markedly from conodes. Such behavior has frequently been observed in exsolved felspar from slowly cooled rocks. Our approach reveals the intrinsic reasons for the specific phase separation path and rigorously describes it by direct numerical solution of the generalized Cahn–Hilliard equation.