Chemical mass transfer in magmatic processes

Contributions to Mineralogy and Petrology

Contrib Mineral Petrol (1987) 96:291-313

9 Springer-Verlag 1987

Chemical mass transfer in magmatic processes III. Crystal growth, chemical diffusion and thermal diffusion in multicomponent silicate melts Mark S. Ghiorso Department of Geological Sciences, University of Washington, Seattle, WA 98195, USA

Abstract. Lasaga's (1982) Master Equation for crystal growth is solved for multicomponent systems in situations which allow for coupled diffusion of melt species. The structure of the solution is explored in some detail for the case of a constant diffusion coefficient matrix. Incorporating these results, the growth of plagioclase is modeled in undercooled tholeiitic melts by approximating interface growth rates with (1) a reduced growth rate function and with (2) calculated solid-liquid solution properties obtained from the silicate liquid solution model of Ghiorso et al. (1983; appendix of Ghiorso 1985). For this purpose algorithms are provided for estimating the liquidus temperature or the chemical affinity of a multicomponent solid solution precipitating from a complex melt of specified bulk composition. Compositional trends in initial solids produced by successive degrees of undercooling are opposite to those predicted in the binary system NaA1Si3Os-CaA12Si2Os. Calculations suggest that the solid phase and interface melt compositions rapidly approach a "steady state" for a given degree of undercooling. Consequently, the overall isothermal growth rate of plagioclase forming from tholeiitic melts appears to be entirely diffusion controlled. In magmatic systems the multicomponent growth equations allow for the formation of oscillatory zoned crystals as a consequence of the "coupling" between interface reaction kinetics and melt diffusion. The magnitude of this effect is largely dependent upon the asymmetry of the diffusion coefficient matrix. Methods are described to facilitate the calibration of diffusion matrices from experimental data on multicomponent penetration curves. Experimental results (Lesher and Walker 1986) on steady state Soret concentration profiles resulting from thermal diffusion in MORB and andesitic liquids are analyzed using the theory of multicomponentlinear irreversible thermodynamics. Under conditions where the entropy production is minimized, a linear relationship is derived between liquid chemical potentials and temperature. This relationship is utilized to evaluate the validity of the solution model of Ghiorso et al. (1983) in melts up to 300 ~ C above their liquidus. The results indicate that configurational entropies are accurately modeled for MORB and andesite bulk compositions. The modeling fails in two four-component systems tested. Equations are derived which allow the calibration of multicomponent regular solution parameters from steady state Soret arrays. An algorithm is demonstrated which permits the calculation of steady state Soret concentration profiles, given an overall bulk melt composition and

temperature gradient. This algorithm uses the liquid solution properties of Ghiorso et al. (1983) and constants obtained from the experimental measurements of Lesher and Walker (1986).

Introduction

In the first two papers of this series (Ghiorso 1985; Ghiorso and Carmichael 1985; hereafter parts I and II) we presented a theoretical analysis of various thermodynamically based numerical algorithms suitable for the calculation of stable solid-liquid phase assemblages in magmatic systems. Our examples (Part II) clearly demonstrate the utility of equilibrium thermodynamic solution models for calculating the consequences of irreversible (disequilibrium) magmatic processes such as crystal fractionation and solid phase assimilation. The algorithms presented in Parts I and II approximate these irreversible phenomena as a series of discrete steps in system bulk composition, temperature or pressure, with each computed node in reaction progress characterized by thermodynamic equilibrium. Although generally valid, this procedure ignores the processes which detail the redistribution of mass between solids and liquid accompanying each irreversible step. These disequilibrium processes are the subject of the present paper. Our understanding of the mechanisms of mass transport in magmas, which include chemical diffusion and crystal growth, is greatly simplified by viewing them referenced to the equilibrium state of the system. In particular, if deviations from equilibrium are sufficiently small (Prigogine 1967; De Groot 1951; De Groot and Mazur 1984), then the fluxes of matter and energy which attempt to restore the system to an equilibrium state can be approximated as linear functions of spatial gradients of appropriate intensive variables (like T, P and chemical potentials) 1. Similarly, rates of chemical reactions, sufficiently close to equilibrium, become essentially proportional to the chemical affinity associated with the reaction; and the affinity itself is a direct measure of departures from the equilibrium state. By way 1 For example the well known applicability of (1) Fourier's law to problems of heat conduction, where the spatial gradient of temperature is proportional to the heat flux, and (2) the phenomenological equivalent of Fick's law to problems of chemical diffusion, where spatial gradients in chemical potential are linearly related to the mass flux of a component

292 of example, affinities associated with the crystallization of metastable magmatic liquids supercooled more than 100~ C below their liquidus are calculated in Part II. Their small magnitudes (5 to 10 cal/g of crystallized solid) when compared to R. T (the gas constant, R, multiplied by the absolute temperature, T) emphasizes that the driving force for crystallization in magmatic systems is small and hints that the rate of crystal growth will be approximately a linear function of chemical affinity and thereby a fairly sensitive function of melt composition. Viewed in this fashion, it becomes obvious that an adequate kinetic description of irreversible phenomena in magmatic systems can only be based upon a comprehensive model describing heterogeneous phase equilibria. Unfortunately, the additional parameters necessary to quantitatively model the "paths" of irreversible processes in magmatic systems are lacking. There is in particular a dearth of data regarding the diffusive transport of melt species in multicomponent silicate melts. Despite these limitations, the purpose of this paper is to develop the theoretical framework and extend the available numerical techniques toward a first step in quantitatively predicting crystal growth rates and thermal diffusion in igneous systems. In the subsequent discussion some of the kinetic and thermodynamic theory involved in modeling the growth of complex non-ideal solid solutions from multicomponent silicate liquids will be outlined and results will be provided on the calculation of the growth of plagioclase from basaltic melts. These growth rates and liquid/solid disequilibrium compositional trends are compared to those predicted by Lasaga (1982), who modeled the growth of plagioclase in the A b - A n binary. It will be seen that the multicomponent phase-loop geometry has a profound affect on the resulting disequilibrium compositional trends. In addition to multicomponent crystal growth calculations, this paper will present equations suitable for analyzing experimental results on chemical and thermal diffusion in melts and will demonstrate how data collected on coupled diffusive processes measured under stationary-state conditions (e.g., steadystate Soret arrays) can be used in the calibration of multicomponent solution models for silicate liquids.

Crystallization kinetics

Approaches and problems Kirkpatrick (1976) has suggested that the rate of increase in the volume of crystals, dV/dt, in a solidifying magma body, evaluated at time t, is given by ~

=4nU(t)Y(t)iI(t')[

! Y(t")dt"]Zdt ',

(1)

where U(t) is the volume fraction of melt remaining uncrystallized at t, Y(t) is the linear interface growth rate (the rate of advance of the crystal surface) and I(t) is the rate of crystallite nucleation at time t. In this equation, it should be noted that the temperature and pressure dependence of U, Y and I is implicitly stated, as in a cooling magma body both T and P are functions of time. Although Eq. (1) is valid only for isotropic, spherical crystals growing from a melt, where the rate of transport of melt components to

the growing crystallite surface far exceeds the rate of consumption of these components (i.e., interface controlled rather than diffusion controlled growth; an assumption which will be shown below to be invalid), it brings to light the mutual importance of nucleation and crystal growth rates in modeling quantitatively the solidification of a magma body. Of the two contributing effects, the rate of crystallite nucleation is perhaps the most difficult to quantify, in that it is dependent upon both the surface energy of a crystallite and on the activation energy for transport of components to the site of nucleation. The interface growth rate, on the other hand, is often written in the form

+ AG* Y(t)=c exp ( i f ) J 1

-A -- exp (R-~)],

(2)

where c is a constant, A G* is the activation free energy for molecular attachment and A is the chemical affinity associated with the crystallization of one mole of solid from the melt (Turnbull and Cohen 1960). The factor c exp (-A G*/R T) can be regarded as proportional to a diffusion coefficient for transport across the melt-crystal interface (Turnbull and Cohen 1960). If this diffusive factor is modeled with a Stokes-Einstein temperature and viscosity (~) dependence, then Eq. (2) can be written: Y(t) = - ~ [1 -- exp (~-TA)],

(3)

where Y~ is defined as the reduced growth rate and is a function of temperature (Kirkpatrick 1975; Jackson 1967). Our subsequent discussion of magmatic crystallization will be concerned with modeling interface growth kinetics in naturally occurring silicate liquids, and will focus attention on the compositional dependence of the term involving the chemical affinity in Eq. (3). This is not to diminish the importance of nucleation effects or the central issue of the rate of magmatic crystallization [e.g., Eq. (1)], but rather to focus in on the contribution the equilibrium description of melt behavior (Ghiorso et al. 1983) can play in describing the time-dependent compositions of minerals crystallizing from a supersaturated silicate melt. Before launching into the full multicomponent coupled reaction-diffusion growth equation, a simple example will illustrate these points. Consider the crystallization of binary albite (Ab)-anorthite (An) plagioclase solid solutions from a melt. For simplicity we will write the two end-member reactions involving hypothetical but convenient melt components: NaA1Si308 = NaA1Si308, (melt) (solid)

(4a)

CaA12Si20s = CaA12Si208. (melt) (solid)

(4b)

If for the purposes of this simple example we make the approximations that (1) the growth rate of the solid solution is independent of the relative rates of the diffusive fluxes and that (2) the composition of the crystallizing solid is given by the relative, uncoupled, interface growth rates of An and Ab, then the ratio of the mole fraction of An and

293 Ab precipitating at the crystal-melt interface is given from Eq. (2) by: -- AAn r=Xab

Yah

[

--Aab '

(5)

where the melt viscosity term has canceled and the affinity terms correspond to Eqs. (4a) and (4b). The reduced growth rates, Y~.An and Y~,Ab, can be considered functions of temperature, or more conveniently the degree of undercooling (A T) below the liquidus. Clearly, far from equilibrium (A ~>RT) r may be written

r = Y~,*n/g,*~, where the composition of the precipitating crystal is only a function of the reduced growth rates. A more likely condition in silicate melts, however, is the situation where A ~ RT (Part II). For sufficiently small A -A

A can be approximated by

A~--(Ah-TA~, where A K is the partial molar enthalpy of the end-member in the solid solution minus the partial molar enthalpy of the melt component. For a pure substance crystallizing from a liquid of identical composition A g and A h would be just the molar entropy and enthalpy of crystallization at T. Substituting into Eq. (7) we have at last: NAn (A gAb XAn~ zJ hAn -- TA SA. XAb = \A SAn XAb]eq A~AbAb-- T A gA~"

Equation (8) reveals that for conditions sufficiently close to equilibrium, such that A/RT,~ 1, the amount of undercooling, T~q-T, directly determines the composition of a crystallizing solid solution. If we substitute the identities NAb -]-XAn = 1, (Xn0o. + (XAn)oq = t,

A

Asi=Ag~ A

[OTln(aisolid/ailiquid)\ ~ )e'

and

RT' and XAn

(8)

9 /aTln(aisolid /a liquid i \ ZJh/= A ~o +gTln(a solid i /a liquia i ) -- RT~ OT )v'

AAn Yr,An

~=Y~.~- AE.~ ~,A~"

(6)

As crystal-melt equilibrium is approached the ratio r tends to the value (XAnfXAb)e q ; the composition of the solid which would form at equilibrium9 In the situation of a fixed amount of undercooling, the ratio of reduced growth rates may be considered constant. This value will be approximated by evaluating the limit of Eq. (6) as A T tends toward zero. In this limit Eq. (6) becomes

where Ag o and A ~ are the molar entropy and enthalpy of crystallization of the pure substance i, into Eq. (8), we obtain: NAn (Ag~ NAb -- ~ e q

/OAAn\

solid liquid _[dTln(aa. /aAn )~ ~-y -)/

[CqAAb\

Yr An (zJ SAn~ = rr,~b \A~bAb]eq,

-- T A gAl')b)'

(9)

(10)

k2=-K~

(11)

liquid), k a -__R T l n ( a a .solid/aan

(12)

and

where ASA. and dgAb are the partial molar entropies of crystallization of An and Ab, respectively, and the subscript equation denotes values at the liquidus temperature. Close to equilibrium, therefore Yr'An_ ( ~ A b / { g A n ~ Yr,Ab -- \A SAn/r kXAb/eq

which when substituted into Eq. (6) yields XAn (ASAb~ (XAn~ AAn ~AhAh= ~A~An]eq ~bAb]e q AAb"

~

where a TIn tas~ v,.t= --R~/ ~ ff~Ab / Ab :/] ,

XAn [XAn\ Yr An AA. lim ~ - ~ = ( ~ ) ~ | = ~ lim Ab \ Abfeq Ir, Ab AT~O NAb

AT~O

rr, Ab

(Ah~176 (XAb)eq (

(7)

Finally, if growth arises from undercooling the melt to a temperature T below the equilibrium temperature, T~q, then

k4 = R T In (a~id/a~quid).

(13)

Calculated values of A g~n, A h-~ A gob and A fiOb are provided in Table 1 at a variety of temperatures and i bar pressure. To pursue this example, let us evaluate Eq. (9) by considering the crystallization of plagioclase in the system CaA12Si2Os-NaA1Si3Os, where to an excellent approximation both liquid and solid-solutions can be treated as ideal. Then Eqs. (113)through (13) can be written: k 1= - - R ln(XAb solid liquid /XAb ), __ solid liquid k2 - - R ln(XA. /XA, ), liquid), k3=RTln(XA.solid/XA. and

(10') (11') (12')

294 Table 1. Plagioclase Compositions as a function of the degree of undercooling (Eq. 9)

Table 2. Plagioclase compositions as a function of the degree of undercooling (Lasaga 1982)

Equilibrium T (~ = 1400 An A H = 30752.2 cals/mol Ab A H = 17893.6 cals/mol

X]~ =0.344152 An AS= 16.76 cals/K-mol Ab AS= 12.84 eals/K-mol

Solid mole fractions An = 0.77713 Ab = 0.22287

Liquid mole fractions An = 0.34415 Ab = 0.65585

T (~

AHOn

AS ~

AU~

AS~

X~~ id

1395 1390 1385 1380 1375 1370 1365 1360 1355 1350 1345 1340 1335 1330 1325 1320 1315 1310 1305 1300

30691.2 30629.6 30567.7 30505.3 30442.3 30378.9 30315.0 30250.6 30185.8 30120.4 30054.6 29988.4 29921.6 29854.4 29786.7 29718.6 29650.0 29580.9 29511.4 29441.4

16.72 16.69 16.65 16.61 16.57 16.54 16.50 16.46 16.42 16.38 16.34 16.30 16.25 16.21 16.17 16.13 16.08 16.04 16.00 15.95

17851.2 17808.6 17765.8 17722.8 17679.6 17636.2 17592.6 17548.9 17504.9 17460.8 17416.4 17371.9 17327.2 17282.3 17237.3 17192.0 17146.5 17100.9 17055.1 17009.1

12.81 12.79 12.76 12.74 12.71 12.68 12.66 12.63 12.60 12.58 12.55 12.52 12.49 12.47 12.44 12.41 12.38 12.35 12.32 12.29

0.77713 0.77713 0.77714 0.77714 0.77715 0.77715 0.77716 0.77717 0.77717 0.77718 0,77719 0.77720 0.77721 0.77723 0.77724 0.77725 0.77726 0.77728 0.77729 0.77731

Data Sources: entropies of fusion for anorthite and albite (Stebbins et al. 1984); heat capacities of molten anorthite and albite (Stebbins et al. 1984); heat capacity of anorthite solid (Krupka et al. 1979); heat capacity of analbite solid (Robie et al. 1978)

k4 = R

T In(XA~/XA~ solid liquid).

(13')

If we assume that the bulk composition of the melt is unaffected by crystallization of the first minute quantities of plagioclase, then the liquid mole fractions can always be equated to their equilibrium values. Thus Eq. (9) becomes a function of one variable, XAn :

[A Ob- R In{XAb]] XA~ (1 -- XA,)

XA,

1407.51 1397.51 1387.51 1357.51

2.07375 2.18878 2.31106 2.72667

0 5.31479E-5 1.91991E-4 1.02268E-3

0.713684 0.753273 0.795355 0.938376

mately 0.78 2. If we examine the predicted plagioclase composition at successive five degree intervals of undercooling we observe a very minor, but increasing enrichment in anorthite content. This anorthite enrichment is certainly not as p r o n o u n c e d as that suggested by L a s a g a (1982), whose calculations are illustrated in Table 2 a, but does run contrary to the expectation that the supersaturated crystalline composition should follow a " s o l i d u s " p a t h which parallels the k n o w n solidus. It is obvious that the growing crystallite should "zone o u t w a r d s " to the equilibrium plagioclase composition at the lower T (Hopper and U h l m a n n 1974). The point of this rather lengthy example is that when Eqs. (9-13) are evaluated for m u l t i c o m p o n e n t melts, there is no valid reason to treat the liquid activity of A n or A b as ideal. Consequently, Eqs. (9-13) demonstrate that just as liquid activity/composition relations have a significant affect on the geometry of the plagioclase loop (Morse 1980, p. 93), they will also significantly affect the disequilibrium compositions of the first feldspar crystallites grown from supersaturated melts. The anorthite enrichment which accompanies undercooling in the binary m a y very well turn to albite enrichment in an undercooled basalt. In order to quantify this effect we must develop a m o r e rigorous m e t h o d of calculating the instantaneous solid compositions and chemical affinity associated with the precipitation of any mineral phase from a m u l t i c o m p o n e n t silicate liquid. The m e t h o d will be independent of the a w k w a r d assumptions of uncoupled growth rates which underlie our simple example and will consider the effect of chemical diffusion on overall growth rates. 2 Bowen gives XAn~0.72, (XSAn/XIAn)eq~--2.058. T h e s e values were adjusted in light of Table 1 in accord with our previous assumption of plagioclase/melt ideal mixing. Thus the values given in Table 1 are consistent with

0 = A h 0 b - T~qAg~ + RT, q ln(X~Ab/XZAb).

TAg~ (9')

L (x~,)oq J

Y(t)

and

"\'Ab/eq -

K

0 = A hOn - Teq A go + gT~q ln(X~n/XL),

kxA V /J~ {xA~

[AgOn--R

T(~

J

Numerical values of XA, can be extracted from Eq. (9') using a root-finder such as Newton's method. If we choose an equilibrium temperature of 1400 ~ C, then the d a t a on the plagioclase phase l o o p provided by Bowen (1913) together with the enthalpies and entropies tabulated in Table 1 allow the evaluation of Eq. (9') as a function of the degree of undercooling. In particular, at 1400 ~ C, (XAn)eq is approxi-

3 Lasaga's anorthite enrichment is so pronounced because of his simplifying assumption relating melt and solid anorthite content to a single distribution coefficient, K, chosen to be the liquid-solid equilibrium constant at the lower temperature. If a similar growth equation is written for the Ab end-member, then the albite mole fraction of the precipitating grain can be calculated analogous to Lasaga's treatment of the An end-member. This generates an inconsistency, in that the sum of the predicted albite and anorthite mole fractions in the solid will never be unity except at the liquidus temperature. For any finite degree of undercooling the calculated concentration of An in the precipitated grain is consequently ambiguous. This unrealistic assumption of a single distribution coefficient does not, however, affect the main conclusions of Lasaga's growth model (e.g., Lasaga 1982, Fig. 5)

295

General growth equations for plagioclase precipitation from multicomponent melts In this section the growth of plagioclase from a multicomponent melt is considered. The equations to be derived can be readily generalized to the estimation of the growth rates of other minerals from magmatic liquids. Let us choose the usual set of thermodynamic components (Ghiorso et al. 1983) to describe the melt composition, and in particular define a vector et which expresses concentrations of Si408, Al16/308, Ca4Si208, Na16/3Si8/308 and K16/3Si8/308 and possibly other melt components, in units of mol/cm 3. Similarly, a vector e~ can be defined whose elements express concentrations of CaA12Si/O8, NaA1Si308 and KA1Si308 in the plagioclase. If Vt and V~ are taken to be the molar volumes of liquid and solid respectively, then a constant matrix of stoichiometric coefficients, T ~, can be defined such that, on an equivalent compositional basis, TVI el= V~c ~.

(14)

It proves advantageous to also consider the plagioclase composition in terms of a vector of mole fractions, x. Doing so a diagonal matrix of reaction coefficients, v, may be defined as: v---E(Trx) e~,

(15)

i

and constant D, to Lasaga's (1982) Eq. (2). As Eq. (16) represents a system of second order partial differential equations, their solution requires the stipulation of initial and boundary conditions. Generalizing Lasaga's (1982) Eqs. (3) and (4) we stipulate an initial concentration profile corresponding to homogeneous melt of composition Co; el(x, 0) = %, and assume that the melt reservoir is infinite el(w, t) = Co.

~c l

X = Z x e r,

(19)

i

where el are unit vectors of length 3, then the solid flux vector at the interface is given by Jso.~ = - Y X c~.

(20)

This flux of solid components must be matched by the liquid flux given by Eq. (15) (21)

where the gradient has been evaluated at the crystal-liquid interface. From the definition of T the fluxes in Eqs. (20) and (21) are related by j=TTjsoUd, or

D0Ct ~x x=o + Y(vct-TrXcs) =0.

(22)

where it is often tacitly assumed that the diffusion matrix is not concentration (i.e., spatially or temporally) dependent in which case it can be brought outside the derivative. This assumption need not be made, and the consequences for arriving at a solution to Eq. (16) for concentration dependent D are discussed in the Appendix. Equation(16) is the basic multicomponent diffusion equation for crystal growth and reduces, for one component

Equations (22), (17) and (18) provide sufficient conditions for the solution of Eq. (16) for a given growth rate, Y The general method of solution for this system of partial differential equations is outlined in the Appendix and is provided in the case of constant D by Eq. (A-23). Evaluation of Eq. (A-23) requires expressions for the solid concentration vector c~ and the interface growth rate Y The vector es can be computed in terms of the vector el if we make the assumption that the concentration in the surface of the growing solid face is always that dictated by exchange equilibrium with the adjacent melt. This requirement is equivalent to maintaining homogeneous equilibrium in the solid and the melt while allowing disequilibrium between the two phases. As shown in the Appendix, the requirement of exchange equilibrium allows the calculation of a consistent set of solid mole fractions while simultaneously extracting a measure of the chemical affinity for the formation of the solid from the melt. By manipulating Eq. (A-11) the ith element of the vector es can be written in the functional form:

4 Here T might be defined such that:

% =f~ (T, P, el(0, t))

[3/8 3/8 1/4 0 ~ ] T~[5/8 3/16 0 3/16 1_5/8 3/16 0 0 3/16

where T and P are the temperature and pressure and el (0, t) is the interface melt composition. The activity/composition relations necessary for he evaluation of Eq. (23) can be cal-

j = --D ~-~-- Yvet,

(15)

(18)

The boundary condition demands that there be equal fluxes at the crystal-melt interface. In particular, if we define the diagonal matrix X as

_ Oct j = - D ~ ~ : o - Yvcl'

where e~ is a unit vector of order 5 or larger (there are at least five liquid components) with a I in the ith row. The elements of v describe the relative amounts of all the liquid components that must diffuse to the crystal liquid interface in order to precipitate the plagioclase of composition x. Following Lasaga (1982) we will assume that this interface is set at the space coordinate x = 0 for all times t. Then the vector of fluxes of liquid components at time t and at distance x from the crystal-liquid interface is given by

(17)

where Y is the interface growth rate (Eq. 3) and D is a matrix of multicomponent diffusion coefficients. It should be noted that because of our definition fixing this interface in space, Y and v are functions only of time. Applying a continuity argument to Eq. (15) we arrive at 0 (D -0cA ~ ) + _Yv Oct 0ct_~xx ax ~t

-

(16)

(23)

296 culated from the silicate liquid solution model of Ghiorso et al. (1983) and the ternary feldspar solution model of Ghiorso (1984). It follows from Eq. (23) that a solution of Eq. (A-23) for the liquid concentration vector at the crystalliquid interface (et (0, t)) implicitly yields the time dependence of the composition of the solid phase forming from the melt. Before proceeding to an evaluation of specific growth rate functions it is useful to discuss the general form of solutions to Eq. (16).

General comments concerning the form of solutions of Eq. (16) for constant D. Equation (A-23), the solution to Eq. (16) for constant D, involves product integrals of exponentials of the sum of the reaction coefficient matrix times the interface growth rate, IT.v, and the diffusion matrix, D. This equation is analyzed in some detail in the Appendix, where it is demonstrated that the evolution of the system away from the initial state (%) will proceed either monotonically or in an oscillatory fashion depending on the eigenvalues of Y.v+D. The evolution will be monotonic, and consequently give rise to normally zoned solid crystals, if the eigenvalues of this matrix sum are real numbers. On the other hand, complex eigenvalues of ~ + D will produce time oscillatory liquid concentrations at the solid interface and consequently oscillatory zoned solid crystals. The wavelength of the oscillations (width of the oscillatory zone) will be related to the size of the imaginary part of the complex eigenvalue and the amplitude (maximal variation in mole fraction across a band) will be a function of the size of the real part. It is therefore of some interest to discover under what condition(s) the matrix sum Yv+D will yield complex eigenvalues and whether such an occurrence can arise from pure diffusion, pure reaction, or is a phenomena characteristic only of coupled reaction-diffusion processes. We first investigate constraints on the contributing matrices to the sum. At constant temperature, Y is a monotonic function of time (e.g., Eq. 3), and the matrix v is diagonal and consequently possesses real eigenvalues. Thus the reaction matrix (Y-v) alone contributes toward making the liquid concentration a monotonic function of time at any distance x from the crystal-liquid interface. Any possibility of an oscillatory solution to Eq. (16) depends either directly on the structure of the matrix D or more subtlety on that of the matrix sum Yv + D. The structure of D can be appreciated by considering the phenomenological diffusion matrix, L, which relates the component flux, j, to the spatial gradient of a vector of liquid chemical potentials, p:

Onsager's (1931) reciprocity theorem requires L to be symmetric. However, H is defined by the activity composition relations of the phase under consideration, and even in the case of an ideal solution H is asymmetric 5. The asymmetry of H is generally a function of the bulk composition of the system and the molar volumes of system components. In the binary ideal case (see footnote), the further from the compositional midpoint of the system the more asymmetric H becomes. The general conclusion to be made of this discussion is that if L is guaranteed symmetric and H is in general asymmetric, then through Eq. (24) the matrix D will generally be asymmetric in near-equilibrium states of the system. Asymmetric matrices are not expected to have real eigenvalues. In near-equilibrium states, however, a further constraint can be placed on the structure of D. Prigogine's (1967) theorem of minimum entropy production demands that L be positive semi-definite (have real eigenvalues greater or equal to zero) in states close to thermodynamic equilibrium. In a remarkable proof published by Kirkaldy et al. (1963) it is demonstrated that the condition for thermodynamic stability implies that the matrix H is always positive semi-definite despite its obvious asymmetry ! Since the product of two positive semi-definite matrices is itself a positive semi-definite matrix, we conclude that the matrix D possesses only real eigenvalues in near equilibrium states. Liquid concentration vectors predicted by Eq. (A-23), under conditions where the interface growth rates are small enough such that diffusion dominates the exponential terms, will be monotonic functions of time and therefore cannot give rise to oscillatory zoned crystals. In those circumstances where the matrices Yv and D contribute in a roughly equal fashion to the matrix exponentials O.e., when the problem is said to be strongly coupled), then concentration profiles predicted by Eq. (A-23) admit the possibility of time oscillatory solutions. The reasons for this are discussed in the Appendix and in Othmer and Scriven (1969). Briefly, the eigenvalues of a matrix sum 5 Consider the simplest case of a two-component ideal solution. H is given by

I I_~Cl ~c2J #t and #2 are given by P l = # ~ + R T l n X 1, and

Dp

j= - L ~ .

#2=#~

As we can also write

The mole fractions may be written: 1 --C 2 /)2

De

1 --C 1 131

j= - D ~ ,

Xl=-l+c2(vl_v2 ) and X2=l+q(v2_v~)

it follows that if H is defined such that

where vl and v2 are molar volumes. The off-diagonal elements of H can be evaluated to yield

Dp Dc ~=H~

O#t

RTX1

0#2

-

then LH=D.

RTX2

and 0c2 c~ vl ~c~ e~ v~ which are clearly not equal (H is not symmetric) unless fortuitously (24)

In diffusive systems close to thermodynamic equilibrium,

V2 X 2 ~/31 X 1

at a particular bulk composition

297 cannot be predicted from the eigenvalues of the matrices that contribute to the sum, unless those matrices have rather restrictive (and in this case physically unrealistic) structures. Thus, the eigenvalues of Yv + D are unconstrained despite the constraints placed above on those of Yv and D. In practical terms, this means that when the rates of chemical diffusion and interface growth are competitive, oscillatory zoned crystals may develop without the necessity of invoking perturbations in the intensive variables of the system. Whether or not this mode of growth predominates depends most critically upon the off diagonal entries of the multicomponent diffusion matrix for magmatic liquids. In this regard, it is quite disappointing to realize the paucity of chemical diffusion coefficient data of any kind for magmatic systems. In preparation for evaluating Eq. (16) to calculate plagioclase growth in basaltic melts, diffusion data are summarized in the remaining paragraphs of this section. Available data on chemical diffusion in multicomponent silicate melts, obtained through the analysis of experimental diffusion couples, has been recently reviewed by Jambon (1983). Data on the diffusion of major elements in anhydrous natural melt compositions can be found in Smith (1974; the diffusion of Na, K), Yoder (1973; the diffusion of SiO2, FeO, CaO, KzO, MgO and TiO2), Alibert and Carron (1980; Si Ti, A1, Na, K, Ca, Mg, Fe) and Medford (1973; Ca) and in hydrous natural melts from Alibert and Delbove (1980; Si, A1) and Watson (1981; Na, Ca). Watson (1987) has recently reviewed diffusivities obtained from mineral dissolution experiments in melts. Such measurements are difficult to evaluate, however, as "apparent" diffusivities obtained from dissolution experiments in melts of fixed bulk composition, appear to be functions of the dissolving mineral (Brearley and Scarfe 1986). The diffusion concentration profiles reported in all of these experimental results indicate a high degree of coupling between diffusive species, giving rise to strongly asymmetric concentration profiles and in some cases "uphill" diffusion (Yoder 1973). Despite these observations, the experimental data have only been analyzed in the context of simple, uncoupled diffusion of a particular melt species referenced to an otherwise inert matrix. This method of analysis is often justified by Cooper's (1968) concept of an "effective binary diffusion coefficient" (EBDC) and the apparent lack of the quantity or quality of data necessary to define the off-diagonal terms of the diffusion matrix. It is therefore only possible to estimate the diagonal elements of D. Even with this conservative restriction there is considerable lack of consistency among experimental results (see discussion of Lasaga 1982). In future studies of diffusion in silicate melts, it would seem preferable to parameterize the experimental data including terms involving diffusive coupling, rejecting the EBDC approach. A representative solution to the multicomponent diffusion problem involving typical experimental boundary conditions is provided for this purpose in the Appendix. An algorithm is suggested for fitting constant and compositionally dependent D matrices to experimental concentration profiles. The argument that the D matrix is underdetermined as a consequence of limited analytical resolution can be answered if the non-linear least squares fits to the data are performed using algorithms that incorporate Singular Value Decomposition (Gill et al. 1981). These least squares techniques and their application to the calibration of D matrices, will be the subject of a future paper. For the purpose of evaluating Eq. (A-23), in lieu of ex-

perimental data, D will be approximatd by a diagonal matrix (see below). Unfortunately, this approximation prevents us from exploring numerically solutions to Eq. (16) involving oscillatory concentration profiles.

Evaluation of Eq. (A-23) for basaltic melts. In order to calculate plagioclase growth in basaltic melts, it is necessary to describe the interface growth function, Y. Lasaga (1982) modeled this time dependent quantity by assuming that the reduced growth rate function, Yr, for anorthite was applicable over the entire plagioclase binary. By calculating the affinity in Eq. (3) in terms of the precipitation of an ideal binary solid solution from an ideal binary melt, Y could be approximated from end-member thermodynamic data and the phase loop geometry (Muan 1979). Estimated growth rates obtained by Lasaga (1982) compare well with those measured by Kirkpatrick et al. (1979). We will assume that Lasaga's (1982) reduced growth rate function applies to the formation of ternary plagioclase crystallizing from basaltic melts. To evaluate Lasaga's Y~, the liquidus temperature, Tz, which corresponds to a particular liquid composition of interest, must be calculated in order that the degree of undercooling may be approximated. For multicomponent plagioclase crystallizing from basaltic melts the calculation of Tt is a formidable task. It involves the determination of the position of the plagioclase saturation surface in temperature/composition space. An algorithm for computing Tt is provided in the Appendix. It utilizes the activity/composition relations of Ghiorso et al. (1983) and Ghiorso (1984) and additionally provides an estimate of the liquidus composition of the plagioclase. The accuracy of T~ calculated in this fashion is estimated to be • ~ C. Evaluation of the chemical affinity for plagioclase precipitation necessitates determining the closest "distance" to the plagioelase saturation surface in temperature/plagioclase-eomposition space at a specific liquid bulk composition and a given degree of undercooling. An algorithm for performing this task is provided in the Appendix. The algorithm is closely related to that for computing the liquidus temperature, and as additional output, it specifies the instantaneous plagioclase composition that would form from the specified liquid at this degree of undereooling. In general, this solid composition will differ from the liquidus composition (see above). Figures 1 through 5 demonstrate solutions of Eqs. (16), (17), (18) and (22) for the case of plagioclase growth from an undercooled melt of olivine tholeiitic composition (Table 3; fractional crystallization example discussed in Part II). D was taken to be of the particularly simple form DI, where I is the identity matrix and D assumes the representative value of 10-6 cmZ/s.Magmatic viscosities were computed from the model of Shaw (1972). Figure 1 illustrates the effect of undercooling on plagioclase composition for crystals growing from olivine tholeiitic magma. The liquidus of the bulk melt is at about 1196~ C and the temperatures plotted on the abscissa correspond to liquidus temperatures of the melt directly adjacent to the crystal-melt interface. The two sets of curves correspond to plagioclase growth from melts supercooled to 1190~ C and 1185~ C respectively. The curves labeled "equilibrium" correspond to hypothetical solid compositions that would precipitate reversibly, at the indicated liquidus temperature from the interface melt composition. The other curves indi-

298 Table 3. Initial composition for plagioclase growth experiments Thingmuli Olivine Tholeiite (Carmichael 1964)" SiO2 TiO 2 A1203 Fe203 FeO

48.47 1.71 15.31 1.66 9.90

MgO CaO Na20 K20 HzO

8.77 11.57 2.30 0.20

0 .9_-g /

0.10

Minor oxides have been neglected and the ferric/ferrous ratio has been adjusted to correspond to the QFM buffer at the temperature of interest (Kilinc et al. 1983). The water content has been approximated (e.g., Part II) a

, 6 ~ , ~ ~ / ~

5~ t

i

,

i

,

=

,

i

i

/

''f

i

t

t

i

i

i

undercooling =

i

equilibrium mole traction in solid

/

Fig. 2. Calculated mole fractions of plagioclase components in crystals precipitating from olivine tholeiitic magma (Carmichael 1964),

at 5~ C undercooling below the plagioclase liquidus, plotted against hypothetical liquidus solid compositions calculated for melts directly adjacent to the crystal liquid interface. The arrows indicate the direction of time evolution of the system. Both the absissca and the ordinate span the following scales: Albite (0.180~.189), Anorthite (0.8054).814), Sanidine (0.005504).00575)

g

0

~

'

,,'88

,,~8

,,~,

Ciquidus temperature

1,~,

'

(*el

Fig. 1. Plagioclase growth from a magmatic liquid of olivine tholeiitic composition (after Carmichael 1964). Liquidus temperature and corresponding "equilibrium" plagioclase anorthite mole fractions are calculated for melt immediately adjacent to the crystalmelt interface. Calculated anorthite mole fractions of plagioclase precipitated from the undercooled melts are parameterized as a function of liquidus temperature (which proxies for time). The arrows indicate the direction of time evolution of the system cate calculated solid compositions for plagioclase precipitated at the crystal-liquid interface, parameterized in terms of "interface liquidus temperatures". Note that as the system evolves toward equilibrium, the two sets of curves converge, and that the 5 and 10 ~ C undercooling "equilibrium" curves are displaced due to a difference in specified Fe203/FeO ratios in the initial supercooled melts. It is clear from Fig. 1 that the composition of the precipitated plagioclase is always more albitic than that which would be produced at equilibrium. This is in direct contrast to the results we obtained above in our analysis of the simple two-component phase loop, and the growth calculations on the binary system reported by Lasaga (1982). The albite enrichment produced in these calculations emphasizes the importance of considering the correct phase loop geometry when modeling plagioclase growth from magmas. Albite enrichment in plagioclase grown from supercooled natural melts has been observed experimentally by Schiffman and Lofgren (1982) and suggested by studies of the Makeopuhi lava lake plagioclase by Evans and Moore (1968). Schiffman and Lofgren (1982) report albite enrichments of roughly 0.22 mol%/~ of undercooling (the data admit numbers between 0.11 and 0.31) for plagioclase grown in isothermal-drop crystallization experiments performed on Columbia River basalt. These numbers can be compared to the calculations reported in Fig. 1 by extrapolating the computed curves to the "steady state" condition where the

composition of the growing grain remains essentially constant. This condition occurs when the interface growth rates are slowed to the point that "interface equilibrium" is approached; the overall rate of crystal growth being essentially diffusion controlled in this state (Lasaga 1982). Intersection of the "equilibrium" and "undercooling" curves in Fig. 1 approximates a limit for this condition. For the 5~ C undercooling calculation reported in Fig. 1, the steady state plagioclase composition is Anso.s2 while that for the 10~ run is An8o.48. The calculated liquidus plagioclase composition is Ansi.33 which results in estimates of albite enrichment on the order of 0.10 mol% and 0.09 mol%/~ of undercooling. These estimates are a factor of two lower than the average of Schiffman and Lofgren (1982). Figure 2 illustrates this feature of albite enrichment in more detail for the 5~ undercooled calculation. Here it can be seen that e n r i c h m e n t in the potassium feldspar component a c c o m p a n i e s that of albite. In Fig. 3 a and b characteristics of the melt directly adjacent to the crystal-liquid interface are displayed. Fig. 3 a illustrates that although the crystallizing plagioclase is more albitic than that which would be produced on the hypothetical melt liquidus, it is still more anorthitic than the eventual equilibrium composition. Hence, initial subtraction of "anorthite" from the liquid is rapid enough to enrich the melt adjacent to the crystal in soda (Fig. 3 a). The rate of crystal growth at this stage is dominated by the interface growth rate. As the interface melt composition is lowered in "anorthite" content, the chemical affinity driving crystal growth (Fig. 3 b) and the interface melt liquidus temperature drop. This lowers the reduced growth rate and the overall interface growth rate, Y (see Fig. 5), so that the rate of growth at the crystal liquid interface begins to be controlled by the rate of c h e m i c a l diffusion, rather than the interface kinetics. Diffusive transport of melt species, in the absence of diffusive coupling, smoothes out the concentration profiles at the crystal-melt interface (Fig. 3 a) as this mechanism begins to dominate the overall crystal growth rate. The shape of the curves displayed in Fig. 3a and b illustrates the rapidity with which diffusion rates dominate the isothermal growth of plagioclase from basaltic melts.

299 ~

q

T Liquid

Concentration

10oC

underr

Table 4. Results of plagioclase growth calculations

r Profiles

o

-6

Liquid temperature (~ Plagioclase liquidus characteristics Temperature (~ Composition X An X Ab Final t (s) Diffusion coefficients (cm-cm/s) Growth rate (cm/s) Initial Final

=

l

o

o,2

0.6

04

a

08

Time (ms)

Affinity (cals/mol) Initial Final

240

11~4

~

|

=

5~ undercooUng

Iso

1191

g .~

Ii~

10*C undercoolJng

1185 ~

10~ undercooling

'

60

5~ undercooling 0.2

o~

~,6

Time

Viscosity Initial Final Precipitated plagioclase composition Initial X An X Ab Final X An X Ab

5~ C of undercooling

10~ C of undercooling

1190

1185

1195.86 0.8133 0.1812 0.1301 1.0.10- 6

1t95.67 0.8140 0.1806 0.001

6.45888.10- 3 2,67690-10- ~

2.10311.10- 2 5,09965.10 - 4

140.891 28.465

256.599 39.400

0.389982 0.384305

0.399446 0.388195

0.8117 0.1826 0.8090 0.1826

0.81ll 0.1832 0.8058 0,1883

(ms}

Fig. 3a, b. Changing liquid concentrations in undercooled melt directly adjacent to a growing plagioclase crystal. The initial bulk composition is that of olivine tholeiitic magma (e.g., Carmichael 1964). a Representative concentration profiles for the liquid components Si4Os (0.3316-0.3334-10-z), Al16/3Os (0.1513~0.1543 910-2), CazSi40 s (0.1393~0.1414 .10-z), Nat6/3Sis/308 (0.38060.3814 910-2) and K16/3Ni8/308 (0.2182-0.2182.10-4). Units are mol/cm 3 and ms. Profiles are coincident except for Na16/3Si8/308. Ii Liquidus temperature (upper two curves, lefthand scale) and chemical affinity driving crystal growth (lower two curves, right-hand scale) for melt directly adjacent to the crystal-liquid interface9 Time is in ms

v-

, " "

t

r

0.O2O

11S4

9

i

J:: 0.010

0"005

1.4t414

--

1188 ~{~

~

~

o ~ G

1185

, o

6o

12o

iso

Alfinily

~4o

~oo

(cals)

Fig. 5. Interface growth rate (lower curves) and melt liquidus temperature (upper curves) plotted as a function of the chemical affinity driving crystal growth at the plagioclase-melt interface. The calculations are plotted for two instances of undercooling of an olivine tholeiitic magma (composition after Carmichael 1964)

1.37900

1o-~

- o

0.2

04

o.8

I oa

Distance (u)

Fig. 4, lsochrons of the concentration of Ca4Si2Os in the liquid as a function of distance (measured in gin) away from the plagioclase-melt interface. Calculations correspond to an undercooling of l0 ~ C for a melt of olivine tholeiitic composition (Carmichael 1964) These features are not critically dependent a p o n the relative magnitude of the diagonal terms of the diffusion coefficient matrix utilized in this example. However, off-diagonal terms in the D matrix will influence the compositional

trends in the precipitated crystals as soon as diffusion dominates the overall growth rate. During this "diffusion controlled" phase of growth, plagioclase compositions are essentially those which would precipitate at the eventual interface equilibrium state of the undercooled system. Solid composition profiles should display a spike of anorthite enrichment followed by a record of "near-equilibrium" diffusion controlled growth. F o r the calculation performed at 10 ~ C of undercooling, Fig. 4 demonstrates the spatial dependence of the concentration of the calcium c o m p o n e n t in the liquid, as a function of time. The 1 ms concentration profile is clearly that expected from a diffusion d o m i n a t e d process. Table 4 provides a s u m m a r y of the results of the plagioclase growth simulations reported in Figs. 1-5. In all p r o b -

300 ability, these calculations should be representative, qualitatively, of the compositional trends accompanying the precipitation of plagioclase from undercooled mafic magmas. The greatest uncertainty in these results centers on the approximation of the diffusion matrix, D. Further quantitative modeling of crystal growth in magmatic systems awaits the experimental determination of the coupled coefficients describing multicomponent chemical diffusion in natural silicate melts. Thermal diffusion

In recent years considerable interest and experimental effort has been expended in exploring the possible significance of the "Soret effect" as a mechanism for magmatic differentiation (Walker et al. 1981; Walker and Delong 1982, 1984). In addition, the utility of using steady state concentration profiles, developed in response to a temperature gradient, to extract information concerning melt solution properties has recently been explored by Lesher and Walker (1986) and Lesher (1986). Measurements of steady-state Soret arrays provide a unique opportunity to investigate melt activity/composition relations over a broad temperature interval and should thus provide excellent constraints on the construction of configurational entropy models for melts. This should come as no great surprise if one considers that the evaluation of any transport phenomena is more critically dependent upon melt speciation than any description of the equilibrium state (i.e., phase equilibrium data). It is the purpose of this section to analyze available steady-state Soret data in the context of the solution model of Ghiorso et al. (1983) and to outline an algorithm for the estimation of Soret separations in magmatic systems. Theoretical description of the steady state

The formalism of irreversible thermodynamics is ideally suited for the analysis of "steady state" Sorer arrays in multicomponent melts. Utilizing this formalism, the steady state can be viewed as one of minimum entropy production for the system. We will write the rate of entropy production (a), corresponding to the flow of matter in response to an imposed temperature gradient, following De Groot (1951): a = (j~-x, + ~ Jk" Xk)/T,

(25)

k=l

where jq represents the flux of heat in the direction of the temperature gradient, x,, across a system consisting of a liquid described in terms of n components; the flux of each of the components, Jk, generates cor/centration gradients, Xk, in the time evolution of the system away from the homogeneous state. In the absence of external forces the gradients in Eq. (25) may be written (De Groot 1951, p. 98) x, = - (grad T)/T,

Ji = ~, Lik Xk + Liu Xu,

(28)

k=l

for all i liquid components, and j~ = ~. L~k xk + L . . x . .

(29)

k=l

Prigogine (1967) has shown that so called stationary or time-invariant states can develop in systems evolving toward equilibrium. These states are not characterized by equilibrium criteria (e.g., homogeneous temperatures throughout the system) but rather by a minimum rate of entropy production [Eq. (25)]. This rate of entropy production is minimized if all fluxes are set to zero in Eq. (25) except that which is associated, or conjugate, to the imposed temperature gradient, x,. Thus, steady state Soret concentration profiles are described by equations of the form: 0 = ~ L~k Xk+ Liu X,,

(30)

k=l

for all i liquid components. It is convenient to define an n by n square matrix L whose elements consist of the plienomenological coefficients modifying the gradients xk (k = 1, n), and a vector 1 whose n elements correspond to the Li,(i = 1, n). Equation (30) can then be rewritten: L [ - Tgrad ( ~ ) ] = - I ( g r a d / ) ,

(31)

where we have used Eqs. (26) and (27) and defined an nvector/t, whose elements consist of the chemical potentials of the liquid components. The matrix L has n - 1 independent rows and columns due to the conservation of mass requirement that the sum of the mass fluxes in the system be zero. Consequently, if Eq. (31) is viewed as a system of linear equations for evaluating gradients in liquid chemical potentials, then only n - 1 of these can be uniquely defined. If one of the liquid component fluxes is neglected, say the nth component, then the system [Eq. (31)] can be uniquely solved. As L is symmetric (Onsager 1931) elimination of the nth component generates a system of equations identical to Eq. (31) where all vectors and matrices are reduced in dimension by 1. These reduced vectors and matrices will be labeled with primes. As the stationary state is one of minimum entropy production the matrix L' must be positive definite 6; L is positive semi-definite (it possesses n - 1 real 6 The rate of entropy production may be written as a quadratic form: x

x

x~

I r]

xl

(26) Xn- 1

and Xk = -- T grad(Irk~T),

form of these relations is

Xu

(27)

where grad denotes the spatial gradient and gk is the chemical potential of the kth component in the liquid. In systems that are "close" to equilibrium, linear phenomenological laws (i.e., Fick's law, Fourier's law) can be written to relate the fluxes in Eq. (25) to the gradients x, and Xk. The general

If G has a minimum at the stationary state then the central matrix in the expression for a must be positive definite (a must be >0 for any choice of gradients; in the case where all gradients vanish, = 0, and we have equilibrium). As any minor of a positive definite matrix is itself positive definite, it follows that L' is a positive definite matrix. A minor of a matrix is any matrix partition formed by deletion of an arbitrary row and column

301 positive eigenvalues and one 0 eigenvalue). Since L' is positive definite its inverse is guaranteed to exist and Eq. (31) may be written

- T grad(K/T) = L ' - t l'(grad T)/T.

(32)

If we now confine our attention to a single spatial dimension, parameterized by x, then Eq. (32) may be written:

1 \dx T

dr _L,_ T2 dx]-

dr1 ll, d ~ T ,

or with rearrangement

d~x'

1 dT =(L'-II'--B') T dx"

(33)

Let us define a vector s equal to L' - 11' and rewrite Eq. (33) for the kth liquid component (dropping the primes):

d$lk 1 dT dx =(Sk--Pk) 7" dx"

(34)

Recognizing that both Pk and T are functions of x, Eq. (34) reduces to the ordinary differential equation

--dpk Sk-- Pk

dT T '

which has a solution of the form ln(Sk-- Pk)= in T+ Ok,

(35)

where ck is a constant. Equation (35) can be rewritten as /1 k =

a k + h k T,

(36)

where ak and b k a r e independent of T (and hence x). We have thus obtained the rather remarkable result that for steady state Soret arrays, each liquid component's chemical potential is a linear function of the temperature profile across the system.

Analysis of steady state Soret arrays We will apply the theory of the previous section to analyze the experimental results of Lesher and Walker (1986). They provide steady-state Soret concentration arrays for two natural and two synthetic compositions. In Lesher and Walker's experiments temperature gradients of 50~ ~ C/ mm were imposed along the axis of cylindrical capsules of approximately 5 mm length. These temperature gradients are approximately linear, and runs were carried out of sufficient duration to demonstrate the achievement of steady state conditions. We have applied the theory of the previous section to Lesher and Walker's experimental results by calculating liquid chemical potentials, at the appropriate T and P, using analyzed melt compositions tabulated for various points along the axis of the charge. The liquid activity/ composition relations are those of Ghiorso et al. (1983). These chemical potentials were then fitted, as a function of temperature, to linear equations of the form of Eq. (36). Coefficients for the equations so obtained are provided in Table 5. Residuals between "measured" and "calculated" chemical potentials for each liquid component can then be

calculated to generate a measure of (1) the quality of the experimental data, (2) the validity of the constant temperature gradient approximation, (3) the appropriateness of the steady state approximation, and (4) the validity of the activity composition relations of Ghiorso et al. (1983). The model residuals are displayed, for each liquid component analyzed in the four Lesher and Walker experiments, in Fig. 6 a-