A Theoretical Treatment of the Precipitation of Doubly Substituted

Article pubs.acs.org/crystal

A Theoretical Treatment of the Precipitation of Doubly Substituted Solid Solutions in Aqueous Solutions C. Noguera,*,†,‡ B. Fritz,§ and A. Clément§ †

CNRS, Institut des Nanosciences de Paris, UMR 7588, 4 place Jussieu, 75005 Paris, France UPMC Univ Paris 06, INSP, UMR 7588, 4 place Jussieu, 75005 Paris, France § Université de Strasbourg/EOST, CNRS, Laboratoire d’Hydrologie et de Géochimie de Strasbourg, 1 rue Blessig, F-67084 Strasbourg Cedex, France ‡

ABSTRACT: We present a theoretical formalism, which, for the first time, accounts for nucleation and growth of mineral particles of variable composition involving two independent substitutions, in aqueous solutions. It is based on the classical nucleation theory, on a size-dependent algebraic growth law allowing growth, resorption, and ripening of particles simultaneously, and on conservation laws akin to a thermodynamically closed system. Truly kinetic conditions are provided to determine the composition of the new nuclei which form and the composition of the deposited layers on the particles which were previously formed. Devised for describing the precipitation of a large range of solid solutions with double substitution, it yields the time evolution of all ion activities in the aqueous solution, together with the particle population characteristics: number, size, and composition profile of particles as a function of time and of their time of nucleation. We show that four different scenarios may take place, depending on the number of end-members involved in the chemical mixing (three or four) and depending on the particle aspect ratio, which may imply quasi two-dimensional growth, like in clay minerals, or three-dimensional growth for more compact solids like salts. This theoretical approach, in the line of our previous works made to treat precipitation of minerals of fixed composition or binary solid solutions, provides a major advance with respect to presently available solid solution precipitation models.

I. INTRODUCTION In the natural environment, crystallization of solid solutions (SS) from aqueous solutions (AS) is an ubiquitous process in which ion partitioning between the AS and the solid phases yields complex chemical compositions which are very dependent on the conditions of formation (temperature, pressure, and fluid composition).1 Indeed, very rarely are natural minerals pure phases. If, for example, one considers calcite, which is the most abundant carbonate mineral in nature, precipitation of solid solutions involving substitution of Ca by one or several divalent metal ions often occurs in oversaturated groundwaters, lakes, or oceans.2,3 The process is of particular importance when the divalent ions are radiotoxic, such as 90Sr, Pb, or Cd, since it allows a retardation of contaminant migration in the natural environment.4,5 Other examples stem for simple salts, such as mixed compounds found in brines which belong to the NaCl−KCl−MgCl2−CaCl2 system,6 or zeolites, which present a large structural and compositional diversity, involving substitutions in the tetrahedral sites, between extra-framework cations and/or between water and vacancies.7 Finally, clay minerals present similar features, with many possible coupled or uncoupled cation substitutions in the tetrahedral, octahedral, and interlayer sites.8 However, while precipitation of pure endmembers has been the subject of much work in the last few decades, understanding of the thermodynamic and kinetic properties of complex solid solution−aqueous solution systems is much less advanced. © 2012 American Chemical Society

Understanding the thermodynamic properties of mineral solid solutions remains a subject of intense activity. From a general viewpoint, relationships between activity coefficients of the end-members and excess Gibbs energies are well established, even in the case of multiple substitutions or in the so-called reciprocal solid solutions.9 Much effort is made to assess the compositional range of stability of the minerals, taking into account thermodynamic factors (temperature, pressure), characteristics of the exchanging ions (their relative radius and charge, their hydration energy), and the possible structural difference of the end-members. Atomistic simulation techniques have been applied to this problem, including ab initio or force field methods,10,11 phonon spectrum calculations for free energies,12 cluster expansion methods or double cluster methods for treating disorder,13 etc. They provide excess mixing properties, information on miscibility gaps and possible ordered structures at intermediate compositions. However, although some exceptions may be found,14 the atomistic works have mainly been limited to binary solid solutions. Understanding SS−AS interactions presents additional difficulties when minerals of variable composition are involved. Many spectroscopic methods nowadays allow the determination of structure, local environments and composition of Received: December 7, 2011 Revised: April 17, 2012 Published: June 1, 2012 3444

dx.doi.org/10.1021/cg201621s | Cryst. Growth Des. 2012, 12, 3444−3457

Crystal Growth & Design

Article

particles produced naturally or by synthesis. A recent review15 emphasizes the difficulties to synthesize clay minerals under low temperature and pressure conditions, due to slow kinetics of precipitation, even if thermodynamic conditions are favorable in the reactors. For that reason, most of the experiments are conducted at higher pressures and temperatures, similar to hydrothermal conditions, but the window for suitable experimental conditions under reasonable aging times may be quite narrow.1 A number of empirical rules derived from experimental or field observations have tried to correlate metal partition coefficients to physical parameters such as the Gibbs free energies of formation of the corresponding metal cations,16,17 the solubility products of the pure end-members,4 or the chemical potential difference between substituent and host cations in solutions.18 As a general statement, very few data exist on solid solution partitioning as a function of precipitation rate,18,19 and the same is true a fortiori when several substitutions are involved. The thermodynamic equilibrium behavior of an SS−AS system is nevertheless a well-established concept.20 The conditions for thermodynamic equilibrium IM = aM = XMλM relate the true saturation state IM of the solution with respect to the M end-member (ratio between the ion activity product and the solubility product of M considered as a pure mineral), the activity aM of M in the SS, its molar fraction XM, and its activity coefficient λM at equilibrium. The fulfilment of as many such equations as end-members is needed to characterize the true equilibrium of the system. In some situations, mainly found when dissolution of a preformed SS takes place, a true equilibrium cannot be reached and the system remains trapped in a metastable state. The concepts of stoichiometric saturation and stoichiometric solubility product have been developed for such situations.20−22 A full kinetic description of the formation of solid solutions in aqueous solutions is still missing. Fritz23 and Tardy and Fritz24 were among the first to apply a solid solution approach for simulating clay-aqueous solution interactions. In their work, at each time step, the amount precipitated just balances the increase of supersaturation, so that the SS always remains at thermodynamic equilibrium with the AS and its composition varies with time. Their model was progressively developed for different types of clay minerals.25−28 For precipitation processes, geochemical codes, like KINDIS 2 9 and PHREEQC2,30 use a single rate equation to describe the time evolution of the total amount precipitated and make an assumption of stoichiometric saturation to account for the composition of solid solutions. As explained in ref 36, these models do not provide an actual kinetic treatment of nucleation and growth. Combining nucleation and growth processes, we have recently developed a formalism which accounts for the precipitation of minerals of fixed composition31−34 and the precipitation of minerals with a single substitution in their framework34,35 or in the interlayer region in the case of clay minerals.36 It goes beyond the previous models, since it combines the classical theory of nucleation of particles with size-dependent kinetic rate laws for growth and/or resorption of particles. The instantaneous value of the saturation state of the aqueous solution is the key parameter which drives the size evolution of stable/unstable particles, and the model keeps track of all particles in the system. In the case of SS precipitation,34,35 at each time step, it makes use of the stoichiometric saturation f unction to know the saturation state

of the aqueous solution with respect to all possible solid phases and kinetic conditions are provided to determine the composition of the new nuclei which form and the composition of the deposited layers on the particles which were previously formed. Although new in the field of geochemistry, a similar approach (limited to the nucleation step) had been previously worked out and validated to describe the formation of binary droplets by condensation in the gas phase.37−40 The present study provides a similar formalism, now relevant for the precipitation of solid solutions involving double exchange. Its numerical implementation, as an extension of our Nanokin code,33 will be the subject of a forthcoming paper. The paper is organized as follows: In section II, we present a classification of doubly substituted solid solutions according to whether substitutions occur on the same site, or on different sites with or without coupling between them. This leads to four scenarios of precipitation which are considered in the four subsequent sections: two-dimensional (2D) nucleation and growth of a four end-member reciprocal solid solution (section III), three-dimensional (3D) nucleation and growth of a four end-member reciprocal solid solution (section IV), twodimensional (2D) nucleation and growth of a three endmember solid solution (section V), and finally three-dimensional (3D) nucleation and growth of a three end-member solid solution (section VI). A discussion follows in section VII and a conclusion in section VIII.

II. CLASSIFICATION We start our analysis of doubly substituted solid solutions with a tentative classification. We will restrict ourselves to minerals which can be formed by precipitation in aqueous solutions. A. The Two Substitutions Occur on Different Sites and Are Uncoupled. When the two substitutions occur on different sites and are uncoupled, the solid solution has four end-members (A, B, C, D), associated with molar fractions XA = (1 − x)(1 − y), XB = x(1 − y), XC = y(1 − x), and XD = xy, and its generic chemical formula may be written as A(1−x)(1−y)Bx(1−y)Cy(1−x)Dxy. The XM (M = A, B, C, D) are linked not only by the usual relationship XA + XB + XC + XD = 1, but also by the cross-product equality XAXD = XBXC, which results from the fact that there are only two independent compositional parameters (x and y). Such solid solutions are referred to as reciprocal solid solutions.9 As will become obvious in the following, even when the mixing in each site is ideal, they may behave nonideally with end-member activity coefficients different from 1. Some examples of reciprocal solid solutions are found in the salt families with generic chemical formula (M,M′)[Y,Y′], in which M and M′ are two isovalent cations (monovalent as Na or K, or divalent as Ca, Sr, etc.) and Y and Y′ two isovalent anions (monovalent as Cl or Br, or divalent as CO3, SO4, CrO4, etc.), apatites, such as (Sr,Ca)5(PO4)3(F,OH), some minerals of the nesosilicate family, etc. Because of their compactness, a 3D nucleation and growth formalism is required. Also belonging to the reciprocal SS family are clay minerals with substitution in the T(tetrahedral) and O(octahedral) sites, and compensation of charge by cations X in the interlayer space. Montmorillonites [Si4−xAlx]T[Al2−2yM3y]OO10(OH)2Xx/m with substitution of Al 3+ by a divalent cation, or [Si4−xAlx]T[Al2−2yM2y]OO10(OH)2Xx/m with substitution of Al3+ by a three-valent cation M3+ (M = FeIII, for example), or [Si4−xAlx]T[Mg3−3yFe3y]OO10(OH)2Xx/m with Mg2+-Fe2+ substitution, etc., belong to this family. Because of their large 3445



Article

Just as an example, let us consider a montmorillonite of chemical formula [Si4−xAlx]T[Al2−2yMg3y]OO10(OH)2Kx. Its four end-members are A = pyrophillite (x = y = 0), B = muscovite (x = 1, y = 0), C = talc (x = 0, y = 1) and D = phlogopite (x = 1, y = 1). According to values quoted in ref 23, at T = 25 °C: KPyrophillite = 101.06, KMuscovite = 1014.56, KTalc = 1021.56, and KPhlogopite = 1038.22. As a consequence, U/RT = −7.28 and α = 10−3.16, a value largely different from 1. In the following, we will make the assumption that the cation interactions within each site are ideal; that is, there is no excess Gibbs energy of mixing ΔGexc. However, we will see in the following that, even in that case, reciprocal solid solutions behave nonideally,9 when U ≠ 0, since at equilibrium the activity coefficients of the end-members in the SS, function of the composition (x,y), are equal to

aspect ratios, these minerals are to be treated with a 2D nucleation and growth formalism (see below and ref 33). Special cases of 2D reciprocal solid solutions are clay minerals with one substitution in the T or in the O site and another in the interlayer space. An example is [Si4−x0Alx0]T[Al2−2yMg3y]OO10(OH)2Kx0(1−x)Cax0x/2, with x0 a constant. The treatment of the interlayer substitution is different from that of the framework substitution, due to different time scales involved, as explained in the following. B. The Two Substitutions Occur on Different Sites and Are Coupled. When the two substitutions occur on different sites and are coupled, the SS have three end-members. Their general formula can be written A(1−x)Bx(1−y)Cxy, with A, B, and C the three end-members. The three molar fractions XA = (1− x), XB = x(1−y), and XC = xy obey the relationship XA + XB + XC = 1. Illites or montmorillonites with coupled exchange in their framework (either tetrahedral or octahedral sites) and in the interlayer region, such as, for example, (Si 4−x Al x )Al2O10(OH)2Kx(1−y)Naxy compounds, belong to this family. Coupling occurs due to charge conservation and a 2D model of nucleation and growth applies. C. The Two Substitutions Occur on the Same Site. When the two substitutions occur on the same site, the SS have three end-members and the generic compound formula is A(1−x)Bx(1−y)Cxy. The XM molar fractions (M = A, B, C) obey XA + XB + XC = 1. Clay minerals with double substitution in either their tetrahedral or octahedral sites belong to this family and their nucleation and growth have to be described with a 2D model. An example among others is [Si4]T[Al,Fe,Mg]OO10(OH)2. Also belonging to this family are more isotropic compounds such as salts (for example, carbonates (Ba,Sr,Cd)CO3 or (Ca,Mg,Mn)CO 3) or some complex oxides. For those compounds, a 3D model of nucleation and growth applies. In the following we will derive the formalism required to describe the precipitation of these various types of SS. We will take into account the number of end-members (four or three) and the dimensionality (2D or 3D) of the nucleation and growth processes which depend on the structure of the minerals under consideration. This leads to four different scenarios, whose formalism is successively established.

ln λA = xy ln α ln λB = −y(1 − x) ln α ln λC = −x(1 − y) ln α ln λD = (1 − x)(1 − y) ln α

The λM are not independent of each other. Since ΔGexc = 0, ΠMλMXM = expΔGexc/RT = 1 or equivalently ∑MXM ln λM = 0 which is easily checked from eq III.2. In the following (sections IIIA, IIIB, and IIIC) we will describe nucleation, growth, and feedback effect on the aqueous solution, respectively, both in the case of double substitution in the framework of the clay mineral or when the y exchange takes place in the interlayer region. A. Nucleation. Since the present section deals with precipitation of clay minerals, nucleation and growth will be considered as 2D processes (see ref 33 and Appendix 2). We will account for nucleation in the framework of the classical nucleation theory.41 The change in Gibbs free energy in the nucleus formation involves the difference between the chemical potentials of the n growth units in the AS and in the SS, proportional to the particle volume, and a surface term. In addition, since the surface energy of the nucleus is expected to vary with composition, the latter cannot be uniform throughout the nuclei but presents some (algebraic) surface enrichment,38−40 so that the Gibbs adsorption isotherm equation (equivalent of the Gibbs−Duheim equation for adsorption processes) is fulfilled.44 We thus introduce algebraic excess quantities nMs of the end-members at the surface of the nucleus and write the total change in Gibbs free energy ΔG(N)(n,x,y) in the formation of a 2D nucleus:

III. FOUR END-MEMBER RECIPROCAL SOLID SOLUTION WITH 2D NUCLEATION AND GROWTH This section mainly concerns clay minerals, described as reciprocal SS of general formula A(1−x)(1−y)Bx(1−y)Cy(1−x)Dxy. A peculiarity of reciprocal SS is the relationship that exists between the true saturation states of the M end-members (M = A, B, C, D): KK IAID ⎛U ⎞ ⎟ = B C = α = exp⎜ ⎝ RT ⎠ IBIC KAKD

(III.2)

ΔG(N )(n , x , y) = Zn1/2v1/2(x , y)σlat(x , y) +

∑ [XMn + nMs]ΔμM M

(III.3)

The first term is the surface Gibbs energy term. Zn v (x, y) represents the lateral area ( (x, y) of the 2D nucleus and σlat(x, y) the surface energy of the lateral faces. Z is a constant factor, function of the number Nlat of lateral faces, the thickness e, and a form factor s of the nucleus (Appendix 2). The number n of growth units in the nucleus is such that nv(x, y) = sel2, with l = (nν(x, y)/se)1/2 the lateral size of the particle and v(x, y) the volume associated to one growth unit of composition (x, y). In the following, n will be referred to as the particle “size”, and we will assume that v(x, y) depends linearly on the end-member growth unit volumes ν(x, y) = ∑MXMνM. 1/2 1/2

(III.1)

The ratio α, independent of the aqueous solution composition, is a function of the solubility products KM and thus of the standard Gibbs free energies of dissolution ΔG0M of the endmembers (KM = exp(−ΔG0M/RT), R the gas constant, and T the temperature) and U is equal to U = ΔG0D + ΔG0A − ΔG0B − ΔG0C. For the sake of readability, we proceed assuming that the stoichiometry coefficients of the sites where substitution takes place are both equal to 1. In Appendix 1, we generalize the approach when this is not the case. 3446



Article

The first term in eq III.9 is equal to ( m(x, y)σlat(x, y)/2; that is, half the surface energy of the nucleus when ΔG(N)(n, x, y) is maximum with respect to n. The determination of the critical nucleus composition requires considering separately two cases, according to whether the y substitution reflects cation exchange in the clay framework or in the interlayer space. 1. y Exchange in the Clay Framework. When y exchange takes place in the clay framework, x and y play equivalent roles and the nucleation barrier can be minimized with respect to both of them. The derivatives of ΔGm(x, y) contain two parts: the first one involves the excess quantities nMs and the derivatives of σlat(x, y); it has to vanish because it represents the Gibbs adsorption isotherm equation.44 The remaining terms, when put equal to zero, lead to the following system of coupled equations, which determine the composition (x*, y*) of the critical nucleus:

The second term in the right-hand side of eq III.3 represents the condensation energy of the growth units which constitute the nucleus. It depends on the change of chemical potentials of the end-members ΔμM = −RT ln(IM/λMXM). The term which does not involve the excess quantities may be recast under the form:

∑ XMΔμM = −RT ln I(x , y) M

(III.4)

with I (x , y ) =

⎛ IM ⎞ XM IAP ⎟ = K (x , y ) ⎝ λMXM ⎠

∏⎜ M

(III.5)

I(x, y) is the instantaneous stoichiometric saturation function, equal to the ratio between the ion activity product IAP and the stoichiometric solubility product K(x, y). Its arguments are the compositional parameters (x, y), but it is also an implicit function of time, via the time evolution of the IM. It allows to know the saturation state of the aqueous solution with respect to all possible solid phases. The deviation of I(x, y) from the value 1 is the driving kinetic factor which will determine all characteristics of nucleation and growth. When I(x, y) > 1 for a given range of (x, y) compositions, formation of nuclei with these compositions may occur. The size of the critical nuclei is determined by the condition that ΔG(N)(n, x, y) is maximum with respect to n, as usual for the nucleation of minerals with fixed composition.41 The maximum values ΔGm(x, y) represent the barrier to nucleation of particles with composition (x, y). It has been recognized37 that, as a good approximation, mostly critical nuclei with the composition (x*, y*) which maximizes the nucleation frequency are formed. In a steady state approximation, the latter reads F = F0(x, y) exp −[ΔGm(x, y)/RT]. In most cases, the dependence of F0(x, y) on composition is unknown, although some functional forms have been proposed in the literature.42,43 The nucleation equations written below assume that F0(x, y) does not vary significantly with composition, so that it keeps a constant value F0. Under this condition, (x*, y*) is determined by the minimization of the nucleation barrier. We give in Appendix 3 the generalized equations when F0 varies as a function of x and y. With respect to n, ΔG(N)(n, x, y) is maximum for n = nm(x, y) such that n m (x , y ) =

(1 − y*)ΔμA + y*ΔμC (1 − y*)vA + y*vC (1 − x*)ΔμA + x*ΔμB (1 − x*)vA + x*vB

( m(x , y) ( m(x , y)

∑ nMs

∂σlat(x , y) + ∂y

∑ nMs

M

M

∂ΔμM ∂x ∂ΔμM ∂y

=0 =0 (III.11)

+ (1 − y*)(nCs + nDs) − y*(nCs ln λC + nDs ln λD) 0=

∑ nMsvM M

(III.8)

0=

and the maximum value of ΔG(N)(n, x, y) is equal to u(x , y)RT ΔGm(x , y) = + ln I(x , y)

∂σlat(x , y) + ∂x

0 = )y − y*(nAs + nBs) + (1 − y*)(nAs ln λA + nBs ln λB)

(III.7)

The associated lateral size of the nucleus is equal to lm(x , y) =

(1 − x*)vC + x*vD

0 = )x − x*(nAs + nCs) + (1 − x*)(nAs ln λA + nCs ln λC) + (1 − x*)(nBs + nDs) − x*(nBs ln λB + nDs ln λD)

2 Z2v(x , y)σlat (x , y )

n m (x , y )v (x , y ) se

(1 − x*)ΔμC + x*ΔμD

allow the excess surface quantities to be determined, when coupled to the definition of the Gibbs dividing surface38 ∑MnMsvM = 0 and to the condition ∑MnMsΔμM = 0, which indicates that the nucleation barrier is independent of the nMs. In the Gibbs dividing surface model, this assumption is consistent with the meaning of excess quantities which do not modify extensive thermodynamic quantities such as free energy and volume. The system of equations which leads to the nMs is thus:

(III.6)

4(kBT )2

=

(1 − y*)vB + y*vD

Note that, in this set of equations, the ΔμM = −RT ln(IM/ XMλM) depend upon x* and y* via the activity coefficients λM and the molar fractions XM. The Gibbs adsorption isotherm equations:

with u(x , y) =

(1 − y*)ΔμB + y*ΔμD

(III.10)

u(x , y) ln 2 I(x , y)

=

∑ nMsΔμM M

∑ nMsΔμM M

(III.12)

with the short-hand notation ) x = Bx(x*,y*) and ) y = By(x*, y*)

(III.9) 3447



Article

RTBx (x , y) = ( m(x , y)x(1 − x)

∂σlat(x , y) ∂x

RTBy (x , y) = ( m(x , y)y(1 − y)

∂σlat(x , y) ∂y

aqueous solution changes with time. For a 2D particle of size n and composition (x, y) in its outer layer, ΔG(G)(n, x, y) is formally similar to the expression written for nucleation, but interfacial ionic activities have to be used in the expressions of the end-member saturation states, the stoichiometric saturation function, and the change in chemical potential of the M endmembers, which will be be denoted Ii(x, y), IMi and ΔμMi = −RT ln(IMi/XMλM), respectively:

(III.13)

2. y Exchange in the Interlayer Space. Cation exchange is considered to be fast when compared to other geochemical processes involved in water−rock interaction processes. According to literature data relevant for clay minerals,45−47 it is usually less than a few hours, while typical time scales of nucleation and growth are in the [20−4500] day range. Treating simultaneously kinetic processes with largely different time scales usually leads to simulations requiring excessive amounts of computer time. A solution to this problem consists of assuming that the fastest process is quasi instantaneous, which, in the present case, amounts to considering that cation exchange reactions in the interlayer region are instantaneous and may be considered at stoichiometric equilibrium with the aqueous solution. Under this assumption, the stoichiometric saturation condition with respect to y may be used, which, expressed as a function of the ΔμM yields: ΔμB = ΔμD ; ΔμA = ΔμC

ΔG(G)(n , x , y) = Zn1/2v1/2(x , y)σlat(x , y) +

M

(1 − y*)vA + y*vC

=

(III.14)

ln Ii(x , y) =

(III.15)

These two equations (eqs III.14 and III.15) allow the determination of (x*, y*). The excess surface quantities nMs are determined from the resolution of the Gibbs adsorption isotherm equation with respect to x, as in eq III.12 (first line). However, regarding y, no excess is expected since stoichiometric equilibrium with the aqueous solution is assumed, which amounts to writing:

n n y* = Ds = Cs 1 − y* nBs nAs

(III.16)

In both cases (exchange in the interlayer region or in the clay framework for y), the size of the critical nucleus n* = nm(x*, y*) reads: n* =

u(x*, y*) ln 2 I(x*, y*)

(III.17)

and the nucleation barrier ΔG* is equal to ΔGm(x*, y*): u(x*, y*) ΔG* = RT ln I(x*, y*)

u(x , y) n

(III.20)

The second conditions yield a set of equations similar to those written for nucleation (eqs III.10, III.12−III.16), but in which interfacial chemical potential changes ΔμMi and the (x, y) composition of the growing layer have to be inserted, instead of ΔμM and (x*, y*). An important point to note is that the size of the particle enters the expression of the composition of the deposited layers. As a consequence, at a given time t of the precipitation process, when the particle population is composed of many particles of different sizes, although the aqueous solution is the same for all of them, the layers deposited on the different particles have different compositions. Interfacial ionic activities are required to solve the growth equations. They are obtained from the conditions of mass balance at the surface of the particle,49 which have to reproduce the (x, y) composition of the growing layer. These conditions involve the diffusion constants of the participating ions in the aqueous solution, but their precise form is highly dependent on the solid solution under consideration. Just as an example, we write them for a montmorillonite of chemical composition [Si4−xAlx][Al2−2yMg3y]O10(OH)2Kx. We denote the anionic part O10(OH)2 of the mineral by the symbol Â and its effective diffusion constant by DÂ . With obvious notations, mass balance at the fluid−particle interface requires the following relationships between the bulk and interfacial molalities:

ΔμB (1 − y*)vB + y*vD

(III.19)

In the following, we will make the assumption of local (stoichiometric) equilibrium at the particle−solution interface. It is valid when short-range transport across the interface is rapid enough to equilibrate the ions in the liquid and solid layers in contact, that is, when the interface motion is slow enough.48 Local equilibrium is reached when the derivatives of ΔG(G)(n, x, y) with respect to n, x, and y vanish. The first condition relates the interfacial stoichiometric saturation function of the aqueous solution Ii(x, y) to the particle size n

while the condition that ΔGm(x, y) be minimum with respect to x reads: ΔμA

∑ [XMn + nMs]ΔμMi

(III.18)

Equations III.10, III.12, III.17, and III.18 provide all the characteristics of the nucleation step for double exchange in the clay framework, while eqs III.14, III.15, III.16, III.17, and III.18 have to be used if y exchange takes place in the interlayer region. B. Growth. The description of growth relies on the expression of the variation of Gibbs free energy ΔG(G)(n, x, y) for the n growth units which are to be incorporated in the particles. Specific equations have to be written to express the change of size of the particle and the evolution of the composition of the deposited layers as the composition of the

DSi [(Si)i − (Si)] = (4 − x)D Â [(Â )i − (Â )] DAl [(Al)i − (Al)] = (2 − 2y + x)D Â [(Â )i − (Â )] DMg [(Mg)i − (Mg)] = 3yD Â [(Â )i − (Â )] DK [(K)i − (K)] = xD Â [(Â )i − (Â )] (III.21)

which are complemented by the expression of the ratio between interfacial and bulk saturation functions: 3448


Crystal Growth & Design ln

Ii(x , y) = I (x , y )

⎛ 1 u(x , y) ⎜⎜ − ⎝ n

Article

⎞ ⎟ ⎟ n m (x , y ) ⎠

frequencies, size, length and compositions are equal to F(t1), n(t1, t), l(t1, t), x(t1, t), y(t1, t), and nMs(t1, t), respectively. In a closed system, the nucleation and growth processes of all particles, which are dependent upon I(x, y), exert a feedback effect on the AS, which modifies I(x, y). The amounts of endmember formula units which have been withdrawn from the AS at time t are equal to

1

(III.22)

If diffusion constants are known, these equations, together with speciation equations in the aqueous medium fully determine all interfacial activities. They take a simplified form if cation diffusion in the AS is much faster than anion diffusion: DĈ ≫ DÂ for all Ĉ cations (Ĉ = Si, Al, Mg, K), due to their small size. This implies that interfacial and bulk molalities are equal (Ĉ i) ≈ (Ĉ ), and the last equation then reduces to ln[Ii(x, y)/I(x, y)] ≈ ln[Â i]/[Â ]). In the time scale of the precipitation process, and especially when it takes place at low temperature, we can assume that solid state diffusion is negligible. The composition of a given layer thus remains fixed once formed. Because of the evolution of the aqueous composition in the course of the precipitation process, the particles thus display inhomogeneities of composition and composition profiles, which can be revealed, post mortem, by electron microscopy characterization.50 All quantities which characterize a particle are functions of two time indices: t1 the time at which the particle was formed, and t the time of observation. The size, length, composition, and surface excesses thus read n(t1, t), l(t1, t), x(t1, t), y(t1, t), and nMs(t1, t), respectively. The same is true for the local saturation indices, which depend on the particle upon consideration: IMi(t1, t). The stoichiometric saturation function, on the other hand, is only dependent on t. In the regime of increasing particle size dl(t1, t)/dt > 0, the rate equation used for the growth of particles of fixed composition31 can be straightforwardly generalized. If growth is limited by interfacial reactions, it reads (τ the short-hand notation for (t1, t), and κ a linear growth constant):

δQ M(t ) =

F(t1)nMs(t1 , t ) dt1

+

∫0

+

∫0

t

F(t1)(n*(t1) − 1)XM(t1 , t1) dt1 t

∫t

F(t1) dt1

t

dt 2

1

dn(t1 , t 2) XM(t1 , t 2) dt 2 (III.26)

The three terms represent the contributions of surface excess, critical nucleus formation, and their (algebraic) growth, respectively. From these quantities, one can calculate all activities in the AS, given a speciation model, and the values of IM(t), IMi(t1, t), I(x, y, t), and Ii(x(t1, t), y(t1, t)). Treating the feedback effect self-consistently allows the full determination of the precipitate and aqueous solution characteristics at all times.

IV. FOUR END-MEMBER RECIPROCAL SOLID SOLUTION WITH 3D NUCLEATION AND GROWTH This section extends the above formalism to four end-member reciprocal solid solutions in which nucleation and growth are 3D processes, as occurs for doubly substituted solid solutions of salts for example (cf. section II). The approach is very similar to that describing double substitution in the clay framework. The main change is in the expression of the surface Gibbs free energy for nucleation and growth (Appendix 2):

⎛ ⎡ u(x(τ ), y(τ )) ⎤1/2 ⎞ dl(τ ) = κ ⎜⎜I(t , x(τ ), y(τ )) − exp⎢ ⎥ ⎟⎟ ⎦ ⎠ ⎣ dt n ( τ ) ⎝

ΔGsurf = Zn2/3v 2/3(x , y)σ ̅(x , y)

(IV.1)

Zn2/3(x, y)v2/3(x, y) represents the effective surface area ( (x, y) of the particle, with Z a geometric factor, function of its shape and σ̅(x, y) its mean surface energy. The change in Gibbs free energy during nucleation is thus equal to

(III.23)

while, if it is limited by diffusion in the AS, κ should be replaced by κ/l(τ). In eq III.23, it is the saturation index relative to the composition x(t1, t), y(t1, t) of the deposited layer, which has to be used. The time variation of the particle size n(t1, t) can be deduced from eq III.23, using the relationship n(τ)v(x(τ), y(τ)) = sel2(τ). During resorption (dl(t1,t)/dt < 0), layers formed at anterior times are progressively dissolved. A layer of lateral size l(t1, t) which reaches the particle/solution interface at a given time t, was deposited at time t2 such that

l(t1 , t 2) = l(t1 , t )

t

∫0

ΔG(N )(n , x , y) = Zn2/3v 2/3(x , y)σ ̅(x , y) +

∑ [XMn + nMs]ΔμM M

(IV.2)

Its maximum value with respect to n occurs for n = nm(x, y) such that nm(x , y) =

(III.24)

t2 is specific to the particle and to the time t and should thus to be written t2(τ). At time t2(τ), the layer composition was equal to x(t1, t2(τ)), y(t1, t2(τ)). Consequently, in the resorption regime, using the short-hand notation τ′ = {t1, t2(τ)}, the growth rate reads:

2u(x , y) ln 3 I(x , y)

(IV.3)

with: u(x , y) =

⎛ ⎡ u(x(τ′), y(τ′)) ⎤1/2 ⎞ dl(τ ) = κ ⎜⎜I(t , x(τ′), y(τ′)) − exp⎢ ⎥ ⎟⎟ ⎣ ⎦ ⎠ dt n(τ ) ⎝

4Z3v 2(x , y)σ ̅ 3(x , y) 27(kBT )3

(IV.4)

and it is equal to ΔGm(x , y) =

(III.25)

u(x , y)RT ln 2 I(x , y)

+

∑ nMsΔμM M

(IV.5)

The first term is equal to one-third the surface Gibbs free energy of the particle: ( m(x, y)σ̅(x, y)/3, ( m(x,y) denoting the value of the nucleus area associated to ΔGm(x, y).

C. Feedback Effect on the Solution. At a given time t, the particle population consists in all the particles which have nucleated at anterior times t1 < t, and whose nucleation 3449



Article

single site of a clay mineral. As previously, we will discriminate the case of y exchange in the interlayer region, which will be assumed quasi-instantaneous, from exchange processes in a site of the clay framework, whose dynamics will be fully accounted for. A. Nucleation. Since 2D nucleation of clay particles is considered here, the Gibbs free energy expression for nucleation is formally identical to what was written for four end-member 2D solid solutions, but with only three terms in the sums over M. The stoichiometric saturation function I(x, y) which enters the expression of ΔG(N)(n, x, y) reads:

The system of coupled equations, which determine the composition (x*, y*) of the critical nucleus is unchanged (eq III.10) as well as the system of equations which yields the nMs (eq III.12), because their formal expressions are only dependent on the number of end-members (four in the present case). The nucleation barrier reads: u(x*, y*) ΔG* = 2 RT ln I(x*, y*)

(IV.6)

which allows estimation of the nucleation rate F = F0 exp −(ΔG*)/(RT) and the size of the critical nucleus: n* =

2u(x*, y*) ln 3 I(x*, y*)

I (x , y ) = (IV.7)

M

Regarding growth, the condition of maximum of ΔG (n, x, y) with respect to n yields the relationship between the interfacial saturation index Ii(x, y) and the size of the particle: (G)

⎛ 2u(x , y) ⎞1/3 ln Ii(x , y) = ⎜ ⎟ ⎝ ⎠ n

(IV.8)

ΔμA vA

=

(1 − y*)ΔμB + y*ΔμC (1 − y*)vB + y*vC

ΔμC − ΔμB

(1 − y)D Na [(Na)i − (Na)] = (1 − x)DCl [(Cl)i − (Cl)]

vC − vB

(1 − y)DK [(K)i − (K)] = xDCl [(Cl)i − (Cl)]

=

(1 − x*)ΔμA + x*ΔμB (1 − x*)vA + x*vB

(V.2)

The excess surface quantities nMs are related to the derivatives of the surface energy σ(x, y), via the Gibbs adsorption isotherm equation. This leads to

(1 − y)DBr [(Br)i − (Br)] = yDCl [(Cl)i − (Cl)] (IV.9)

Added to speciation equations in the aqueous solution and to the ratio of interfacial and bulk saturation states: ⎞ ⎛ 1 Ii(x , y) 1 ⎟⎟ = (2u(x , y))1/3 ⎜⎜ 1/3 − 1/3 I (x , y ) n m (x , y ) ⎠ ⎝n

(V.1)

As a consequence of the assumption of ideality, the activity coefficients λM are equal to 1. The characteristics of the nuclei at the maximum of ΔG(N)(n, x, y) with respect to n, nm(x, y), u(x, y), lm(x, y) and ΔGm(x, y) are thus given by eqs III.6−III.9. On the other hand, as in section III, the determination of the composition of the critical nuclei requires discriminating whether the y substitution takes place in the interlayer space or in the clay framework. 1. y Exchange in the Clay Framework. When y exchange takes place in the clay framework, the composition x* and y* of the critical nuclei is obtained by minimizing ΔGm(x, y) with respect to x and y, neglecting the variations of σ(x, y) and the nMs. This leads to the system of coupled equations:

All equations which determine the composition of the deposited layer during growth are identical to those during nucleation (eq III.12), but in which interfacial chemical potential changes ΔμMi and composition (x,y) of the growing layer have to be used, instead of ΔμM and (x*, y*). The conditions of mass balance at the solution−particle interface, which allow interfacial activities to be determined, are specific to each individual case. Just to exemplify them for one 3D reciprocal solid solution, we write them for (Na,K)(Cl,Br) with obvious notations for the diffusion constants Di:

ln

⎛ IM ⎞ XM ⎟ ⎝ XM ⎠

∏⎜

0 = )x − x*nAs + (1 − x*)(nBs + nCs) 0 = )y − y*nBs + (1 − y*)nCs 0=

(IV.10)

∑ nMsvM (V.3)

M

they provide a full determination of all interfacial activities. The equations for the size variation of the particles with time, whether in the growth regime or in the resorption regime, are unchanged (eqs III.23 and III.25), as well as equations of feedback on the aqueous solution (eq III.26).

with ) x and ) y defined as in section III (eq III.13). The last equality corresponds to the condition of Gibbs dividing surface. The nucleation barrier does not depend on the nMs since the latter fulfill the equation ∑MnMsΔμM = 0. 2. y Exchange in the Interlayer Space. When y exchange occurs in the interlayer space, the stoichiometric saturation condition with respect to y has to be used:

V. THREE END-MEMBER SOLID SOLUTIONS WITH 2D NUCLEATION AND GROWTH This section deals with the nucleation and growth of clay minerals that may be described by an ideal solid solution with three end-members A1−xBx(1−y)Cxy. This may be the result of a substitution in the interlayer region correlated to a substitution in either T or O site, or to a double substitution in a single site. A typical example of the first case is the montmorillonite [Si4−xAlx][Al2]O10(OH)2Tx(1−y)/nVxy/m with a substitution x in the tetrahedra and compensation of charge in the interlayer region due to a mixture of Tn+ and Vm+ cations, possibly non3+ isovalent. [Si 4] T [Al 2−2y Fe 2xy Mg 3y(1−x) ] OO 10 (OH) 2 is one (among others) representative of double substitution in a

ΔμB = ΔμC

(V.4)

and together with the minimization of ΔGm(x, y) with respect to x (first equation of V.2), which, due to the equality ΔμB = ΔμC, simplifies into ΔμA vA

=

ΔμB (1 − y*)vB + y*vC

(V.5)

The excess surface quantities nMs are determined from the Gibbs adsorption isotherm equation with respect to x, as in the 3450



Article

Figure 1. Two representations of the stoichiometric saturation function I(x, y) of an aqueous solution with respect to the (Sr,Ba,Pb)CO3 solid solution, as a function of x and y. Left panel: thermodynamic equilibrium I(x0, y0) = 1. Right panel: supersaturation I(x0, y0) = 100. The composition (x0, y0) at thermodynamic equilibrium is represented by a white triangle, and the composition of the critical nucleus (x*, y*) by a black triangle. In both graphs, the activities of the cations (in mol/kg H2O) are equal to aSr = 1.31 × 10−4, aBa = 1.43 × 10−4, and aPb = 6.0 × 10−9. The solubility products of the end-members are equal to KSrCO3 = 10−9.271, KBaCO3 = 10−8.562, KPbCO3 = 10−13.13. The activity of CO3 anions is adjusted to obtain the desired value of I(x0, y0).

first equality of eq V.3, while the excesses associated with y fulfill the stoichiometric saturation condition: nCs y* = nBs 1 − y*

The systems of coupled equations allowing the determination of the composition of the critical nucleus or the composition of the deposited layer during growth are only dependent on the number of end-members of the solid solution and are thus similar to those written in section V. Only the interfacial saturation state, which is dependent upon dimensionality has to be changed into its 3D form:

(V.6)

In both cases (exchange in the interlayer space or in the clay framework), the nucleation barrier and the size of the critical nucleus read: u(x*, y*) ΔG* = RT ln I(x*, y*)

n* =

⎛ 2u(x , y) ⎞1/3 ln Ii(x , y) = ⎜ ⎟ ⎝ ⎠ n

u(x*, y*) ln 2 I(x*, y*)

(VI.1)

The determination of interfacial activities follows the same lines as in previous sections, as well as the treatment of the feedback effect on the aqueous solution.

(V.7)

B. Growth and Feedback Effect on the Aqueous Solution. The derivation is formally similar to that performed in section III. The relationship between the interfacial saturation state Ii(x, y) and the particle size is only dependent upon dimensionality and is thus identical to that for 2D growth of a four end-member reciprocal solid solution (eq III.20) The composition of the deposited layer and surface excesses are given by equations formally identical to those for nucleation (eqs V.2 and V.3 for y exchange in the clay framework or eqs V.4−V.6 for y exchange in the interlayer region) with (x*, y*) being replaced by (x, y) and ΔμM by ΔμMi. The determination of interfacial activities follows the same lines as in section III, as well as the treatment of the feedback effect on the aqueous solution.

VII. DISCUSSION In the four previous sections, we have provided a formalism which describes nucleation and growth of solid solutions with double substitutions. Although it relies on some approximations, it represents a major advance on various aspects, compared to treatments which are used in some popular geochemical codes, such as KINDIS29 and PHREEQC2.30 In the nucleation step, it takes into account the correct relationship between the critical nucleus size and the stoichiometric saturation function I(x, y) of the solution, rather than assuming, as in previous kinetic models, the sudden appearance of an arbitrary amount of the new phase as soon as I(x, y) exceeds the critical value for nucleation. In the growth regime, a size-dependent law is used (rather than a law proportional to (I − 1) to some power), which allows supercritical particles to increase in size and subcritical ones to resorb, in the process of Ostwald ripening.51,52 Only in that way can thermodynamic equilibrium be reached (at infinite time) with a single particle gathering the whole available matter. Through the present formalism, it is possible to obtain the time evolution of the crystal size distribution function (CSD) or keep track of the number of particles nucleated at each time t1 and their subsequent evolution of size. The preceding versions of our Nanokin code already involved such characteristics. However, they were only valid for precipitation of minerals of fixed composition,33 or binary solid solutions (with exchange in the framework of a particle35 or in the clay interlayer region36). The present derivation allows an extension of this code for

VI. THREE END-MEMBER SOLID SOLUTIONS WITH 3D NUCLEATION AND GROWTH In this section, more isotropic solid solutions such as some carbonates, for example, (Sr, Ba, Pb)CO3 with aragonite structure, sulfates, for example, (Ba, Sr, Zn)SO4 with barite structure, or complex oxides of compact shape are considered, for which a 3D model of nucleation and growth applies and double substitution occurs in a single site. Regarding nucleation, the expression for ΔG(N)(n, x, y) is formally similar to that written for four end-member solid solution with 3D nucleation, with the mere difference that the sum over the end-members M involves only three terms. The critical nucleus characteristics nm(x, y) and ΔGm(x, y) at the maximum of ΔG(N)(n, x, y) are thus given by the same relationships (eqs IV.3−IV.5). 3451



Article

applications to a wide range of doubly substituted solid solutions. An important advance of the present approach is in the treatment of the compositions of the critical nuclei and deposited layers. The composition (x*, y*) of the critical nuclei is determined by a kinetic criterion, namely, that it maximizes the nucleation frequency. Similarly, the composition of the layers deposited at time t on the particles is determined by the minimization of the change of Gibbs energy ΔG(G) for the growth units which enter the layer. It depends on interfacial activities and on the size of the particles. As a consequence, the layers deposited on the different particles have different compositions. These compositions thus differ from the stoichiometric equilibrium values (x0, y0) deduced from the equalities IM = λMXM. We have previously exemplified such a situation in our simulation of (Ba,Sr)CO3 and (Ba,Sr)SO4 precipitation.35 In a Roozeboom plot, the time evolution of their x* composition started far from the stoichiometric equilibrium line and merged the latter only at “infinite” times when thermodynamic equilibrium was recovered. The same situation is found here. Although we are not yet able to perform a full simulation of the precipitation of a doubly substituted solid solution, including growth and feedback effect on the aqueous solution, we have solved the equations which give the critical nucleus composition (x*, y*) in the case of (Sr, Ba, Pb)CO3, using tabulated values for the end-member solubilities and molar volumes.30 Figures 1 and 2 highlight the behavior of the stoichiometric saturation function I(x, y) and the critical nucleus composition (x*, y*), when the solid solution is in contact with hypothetical aqueous solutions, in which the activities of the cations are fixed and the CO3 activity aCO3 is varied. This is a particularly simple case, since the

function I(x, y) is multiplied by a constant factor when aCO3 increases. From Figure 1, left panel, which represents its dependence at thermodynamic equilibrium, it clearly appears that I(x, y) has a convex shape, and that a single (x, y) couple fulfills I(x0, y0) = 1, while for all other couples I(x, y) < 1. The composition of the solid solution at thermodynamic equilibrium can thus be obtained from the simultaneous condition of maximum of the stoichiometric saturation function I(x, y) and I(x, y) = 1. Of course, writing these conditions, using eq III.5, leads to the known set of equalities IM = λMXM for all endmembers M. In Figure 1, right panel, the I(x, y) function is represented in the case of a high supersaturation state I(x0, y0) = 100, for which particles of all compositions (x, y) may precipitate. However, the kinetic condition that the nucleation frequency be maximum selects a single composition (x*, y*), very different from (x0, y0). To better visualize how systematically (x*, y*) deviates from (x0, y0), in Figure 2, we have represented the evolution of (x*, y*) as I(x0, y0) varies from 100 to 1, thus simulating a return to thermodynamic equilibrium. This has been done for eight different aqueous solution compositions, involving different cation activity ratios. In each case, the path of (x*, y*) in the (x, y) plane ends at the thermodynamic equilibrium values (x0, y0). Cation activities in path #1 were used to plot Figure 1. As may have appeared in the course of the theoretical sections, it was fruitful to present the four scenarios of solid solution precipitation together, because they share common features (Table 1). The dimensionality of the precipitation process (2D or 3D) is reflected in the functional dependence of the surface Gibbs free energy with respect to the particle size n, entering the changes in Gibbs free energy for nucleation and growth ΔG(N)(n, x, y) and ΔG(G)(n, x, y). The n1/2 and n1/3 dependencies of the surface contributions, akin to 2D and 3D processes, respectively, are responsible for the expressions of the critical nucleus size, the nucleation barrier as a function of the saturation state of the solution, and the interfacial saturation state as a function of particle radius: u(x*, y*) ΔG* = RT ln I(x*, y*)

n* =

u(x*, y*) ln 2 I(x*, y*)

⎛ u(x , y) ⎞1/2 ln Ii(x , y) = ⎜ ⎟ ⎝ n ⎠

(VII.1)

in 2D with u(x,y) given by eq III.7, and u(x*, y*) ΔG* = 2 RT ln I(x*, y*) ⎛ 2u(x , y) ⎞1/3 ln Ii(x , y) = ⎜ ⎟ ⎝ ⎠ n Figure 2. Evolution of the (Sr,Ba,Pb)CO3 critical nucleus composition (x*, y*) as I(x0, y0) varies from 100 to 1, for eight different aqueous solution compositions, characterized by different ratios of cation activities (in mol/kg H2O). Path #1: aSr = 1.31 × 10−4, aBa = 1.43 × 10−4, aPb = 6.0 × 10−9; Path #2: aSr = 1.36 × 10−5, aBa = 5.73 × 10−5, aPb = 4.80 × 10−10; Path #3: aSr = 1.36 × 10−4, aBa = 5.73 × 10−3, aPb = 2.40 × 10−9; Path #4: aSr = 1.36 × 105, aBa = 5.73 × 10−4, aPb = 2.40 × 10−9; Path #5: aSr = 1.36 × 10−5, aBa = 5.73 × 10−5, aPb = 1.01 × 10−9; Path #6: aSr = 1.36 × 10−5, aBa = 5.73 × 10−5, aPb = 2.40 × 10−9; Path #7: aSr = 1.36 × 10−5, aBa = 5.73 × 10−6, aPb = 2.40 × 10−9; Path #8: aSr = 9.49 × 10−5, aBa = 5.73 × 10−5, aPb = 6.96 × 10−9.

n* =

2u(x*, y*) ln 3 I(x*, y*)

(VII.2)

in 3D, with u(x,y) given by eq IV.4. On the other hand, the equations obeyed by the composition of the critical nuclei or growing layers (x, y and surface excesses nMs) are not dependent on dimensionality but rather on the number of end-members. They are thus identical for 2D and 3D SS with four end-members, and similarly for 2D or 3D SS with three end-members. And since local equilibrium is assumed during growth, their expressions for the critical nucleus and the growing layer are formally equivalent, provided local activities are used in the latter. 3452



Article

Table 1. Scenarios of Precipitation of Doubly Substituted Solid Solutions from an Aqueous Solution, with Some Examples type of SS

2D minerals

3D minerals

4 end-member SS

y exchange in clay framework [Si4−xAlx]T[Al2−yMy]OO10(OH)2X(x+y)/m y exchange in interlayer space [Si4−x0Alx0]T[Al2−2xMg3x]OO10(OH)2Kx0(1−y)Cax0y/2

(Na,K)(Cl,Br)

3 end-member SS

y exchange in clay framework [Si4]T[Al2−2yFe3+ 2xyMg3y(1−x)]OO10(OH)2 y exchange in interlayer space (Si4−xAlx)Al2O10(OH)2Kx(1−y)Naxy

(Sr,Ba,Pb)CO3

which will allow actual simulations to be performed. It may also be useful for describing precipitation of solid solutions with double substitutions produced by precipitation from a melt, as often occurs in metallurgy or in the geological environment.

We have focused here on precipitation processes in an aqueous solution. Some solid solutions can however be formed in other conditions, in particular, by crystallization from a melt. This is a very often encountered situation in metallurgy where alloys are produced from metallic melts. In the geological environment, complex minerals may also be formed by crystallization from magmas. The driving force for nucleation (difference of chemical potential between the crystalline and liquid phases) is then a function of the undercooling, difference between the actual temperature, and the melting temperature. In a previous work focusing on minerals of fixed composition,31 we have shown that under some simplifying assumptions, the mathematical description of nucleation and growth is similar whether precipitation takes place in an aqueous solution or results from crystallization from a melt. This suggests that the present work may also be useful in the treatment of first-order crystallization processes of complex solid solutions.

■

APPENDIX 1 In this appendix, we give the generalization of the nucleation and growth equations for the precipitation of a solid solution in which the stoichiometric coefficients of the sites on which substitution occurs are not equal to 1 but may take any p and q values. The main difference lies in the expression of the configurational mixing entropy of the solid solution: S = −Rp[x ln x + (1 − x) ln(1 − x)] − Rq[y ln y + (1 − y) ln(1 − y)]

(X.3)

which can no longer be expressed as a function of the molar fraction XM as in the case when p = q = 1:

VIII. CONCLUSION We have presented a theoretical formalism which, for the first time, accounts for nucleation and growth of particles of variable composition involving two independent substitutions, in aqueous solutions. As in our previous works on the precipitation of minerals of fixed composition or binary solid solutions, it relies on the classical nucleation theory, on a sizedependent algebraic growth law allowing growth, resorption, and ripening of particles simultaneously, and on conservation laws akin to a thermodynamically closed system. Truly kinetic conditions are provided to determine the composition of the new nuclei which form and the composition of the deposited layers on the particles which were previously formed. Devised for describing the precipitation of a large range of doubly substituted solid solutions, it yields the time evolution of all ion activities in the aqueous solution, together with the particle population characteristics: number, size, and composition profile of particles as a function of time and of their time of nucleation. We have shown that four different scenarios may take place, depending on the number of end-members involved in the chemical mixing (three or four) and depending on the shape of the particles which are produced. A typical application includes the formation of clay minerals, with cationic substitutions, either coupled or uncoupled, in their framework and/or in their interlayer region. For them, a 2D description of the processes of nucleation and growth applies. For more compact solids like salts or complex oxides such as ilmenites, etc., the 3D formalism has to be used. This theoretical approach provides a major advance with respect to presently available solid solution precipitation models. It is currently being introduced in our Nanokin code,

S = −R ∑ XM ln XM

(X.4)

M

but rather reads: S = −R ∑ XM ln YM

(X.5)

M p

q

p

q

p q

with YA = (1−x) (1−y) , YB = x (1−y) , YC = (1−x) y , and YD = xpyq for a four end-member reciprocal solid solution and by YA = (1−x)p, YB = xp(1−y)q, YC = xpyq for a three endmember solid solution. Using this notation, the stoichiometric saturation state of the aqueous solution with respect to a solid solution of composition (x, y) reads: I (x , y ) =

⎛ IM ⎞ XM ⎟ ∏⎜ λ Y M ⎝ M M⎠

(X.6)

With the hypotheses we made, the activity coefficients λM are given by eq III.2 in the case of reciprocal solid solution and are equal to 1 in the case of ideal three end-member solid solutions. The equilibrium conditions read: IM = λMYM

∀M

(X.7)

λMYM is the activity of the end-member M in the solid solution at equilibrium. The critical nucleus characteristics are derived as in the main text. The equations giving the composition of the critical nuclei, expressed as a function of the ΔM, are unchanged but ΔM=−RT ln[IM/λMYM]. They are thus given by eqs III.10, III.14−III.15, V.2, or V.4−V.6, depending upon the type of SS under consideration or by the equations given in Appendix 3 if one takes into account the dependence of F0 on composition. 3453



Article

the substrate (Wadh the adhesion energy). Their aspect ratio is then given by the Wulff-Kaishev theorem:

The excess quantities, on the other hand, depend upon the site stoichiometry, in the following way for four end-member solid solutions:

σ − 0.5Wadh e = bas l σlat

) 0 = x − x*(nAs + nCs) + (1 − x*)(nAs ln λA + nCs ln λC) p + (1 − x*)(nBs + nDs) − x*(nBs ln λB + nDs ln λD) 0=

0=

If, instead of the basal face, one of the lateral faces is in contact with the substrate, the equilibrium shape involves three inequivalent dimensions l, l′, and e. The ratio between these lengths is then equal to

)y

− y*(nAs + nBs) + (1 − y*)(nAs ln λA + nBs ln λB) q + (1 − y*)(nCs + nDs) − y*(nCs ln λC + nDs ln λD)

e σbas

∑ nMsvM ∑ nMsΔμM M

(X.8)

and as 0= 0= 0=

)x − x*nAs + (1 − x*)(nBs + nCs) p )y q

− y*nBs + (1 − y*)nCs

∑ nMsvM M

=

l l′ = σlat σlat − 0.5Wadh

(X.9)

for three end-member solid solutions. Similar relationships are derived for growth, as explained in the main text.

■

ΔGsurf = Aσ + Ai (σ − Wadh)

APPENDIX 2 This appendix justifies the generic expressions of the surface energy term entering the change in Gibbs free energy during nucleation or growth which was written ΔGsurf = Zn2/3v 2/3(x , y)σ ̅(x , y)

(XI.6)

The equilibrium shape is a spherical cap of curvature radius ρ and wetting angle θ with the substrate. θ is given by the YoungDupré equation −σ cos θ = σ − Wadh.54 For a strong adhesion to the substrate, the wetting angle is equal to 0°, while it is equal to 180° when no wetting occurs (which is also the case for homogeneous nucleation). Writing:

(XI.1)

in the case of 3D nucleation (eq IV.1) and ΔGsurf = Zn1/2v1/2(x , y)σlat(x , y)

(XI.5)

With these elements in mind, one can consider heterogeneous nucleation of particles of various shapes. Usually having very few quantitative information on the anisotropy of surface energies for minerals in contact with aqueous solutions, it is useful to make the tentative assumption that the particles display their equilibrium shape during the processes of nucleation and growth. We will see, in the following, that this assumption allows a single effective surface energy to enter the formalism of nucleation and growth. One character common to all the cases under consideration in the following is that the nucleation barrier is lowered when particles nucleate on a foreign solid (surface Gibbs free energy is reduced), so that the probability of heterogeneous nucleation is usually much higher than the probability of homogeneous nucleation.41 Spherical Particles on a Substrate. Considering a particle which has a contact area A with the aqueous solution and Ai with the substrate, the total surface energy involves a surface and an interface terms.

M

0=

(XI.4)

A = 2πρ2 (1 − cos θ ); Ai = πρ2 sin 2 θ n = 4πρ3 Φ(θ )/3v ; Φ(θ ) = (1 − cos θ )2 (2 + cos θ )/4

(XI.2)

in the case of 2D nucleation (eq III.3). We first consider the 3Dcase.

(XI.7)

ΔGsurf reads:

A. 3D Particles

ΔGsurf = σn2/3v 2/3(36π Φ(θ ))1/3

Most solids, except amorphous ones, are non-isotropic and their external shape reflects the relative energies of their low index faces, as recognized by Wulff.53 Indeed, Wulff theorem states that, at equilibrium, the distance from the center of a particle to its external facets is proportional to the surface energy of these surfaces. For example, according to Wulff theorem, the aspect ratio of tetragonal particles with square basal faces (basal dimensions l × l and thickness e) is given by the ratio between the surface energies of the basal and lateral faces: σ e = bas l σlat (XI.3)

(XI.8)

Tetragonal Particles. Tetragonal particles nucleating with their basal face in contact with a substrate have a volume V = nv = l2e and an aspect ratio given by the Wulff−Kaishev theorem: e/l = (σbas −0.5Wadh)/σlat. The surface contribution to the change in Gibbs free energy during nucleation is equal to ΔGsurf = l2(2σbas − Wadh) + 4leσlat, that is, 2 ΔGsurf = 6n2/3v 2/3(σlat (σbas − 0.5Wadh))1/3

(XI.9)

If tetragonal particles nucleate on one of their lateral faces, their volume is equal to V = nv = ll′e and ΔGsurf = 2ll′σbas + 2l′eσlat + le(2σlat − Wadh)

This result can be extended to the case of particles in equilibrium with a substrate on which they lie on their basal face. In this case, σbas is relevant for the face in contact with the aqueous solution, and σbas − Wadh for the one in contact with

(XI.10)

Considering the equilibrium ratio between l, l′, and e given by the Wulff−Kaishev theorem, ΔGsurf is equal to ΔGsurf = 6n2/3v 2/3(σlatσbas(σlat − 0.5Wadh))1/3 3454

(XI.11)



Article

Hexagonal Particles. A similar derivation can be made for hexagonal particles nucleating on their basal face (side length l, thickness e). The volume of the particles is equal to V = nν = 3el2√3/2. At equilibrium, the ratio between l and e is given by the Wulff−Kaishev equation: e/l =

3 (σbas − 0.5Wadh)σlat

2) for tetragonal or hexagonal particles, respectively. It is important to note that in such a situation, the surface energy of the basal face no longer enters the nucleation and growth equations.

■

(XI.12)

ΔGsurf is equal to ΔGsurf = 3√3l (σbas − 0.5Wadh) + 6leσ1at. As a consequence: 2

In this appendix, we generalize the expression for the composition of the critical nuclei, when the prefactor F0(x, y) of the nucleation frequency varies with (x, y). When this is the case, the conditions of maximum for the nucleation frequency do not reduce to the conditions of minimum of the nucleation barrier, but rather read:

2 ΔGsurf = 361/3 3 n2/3v 2/3(σlat (σbas − 0.5Wadh))1/3

(XI.13)

To summarize, when 3D nucleation takes place, the change in surface Gibbs free energy can be written under the general form: ΔGsurf = Zn2/3v 2/3σ ̅

(XI.14)

RT

with (Φ(θ) given by eq XI.7): Z = (36π Φ(θ ))1/3 ; σ ̅ = σ for spherical particles,

APPENDIX 3

∂ log F0(x , y) ∂ ln F0(x , y) ∂ΔGm ∂ΔGm = = ; RT ∂x ∂x ∂y ∂y (XII.1)

(XI.15)

Whether for 2D or 3D nucleation:

2 Z = 6; σ ̅ = (σlat (σbas − 0.5Wadh))1/3

∂ΔGm n RT ⎡ ∂v ∂ ln I ⎤ = m ⎢ln I −v ⎥ ∂x ∂x ⎦ v ⎣ ∂x

(XI.16)

for tetragonal particles in contact with a substrate on their basal face, Z = 6; σ ̅ = (σlatσbas(σlat − 0.5Wadh))1/3

with nm given by eqs III.6 and IV.3 for 2D and 3D nucleation, respectively. As a function of the ΔμM, and using v(x,y) = ∑MXMvM, the previous expression reads:

(XI.17)

for tetragonal particles in contact with a substrate on their lateral face, and finally: 2 Z = 361/3 3 ; σ ̅ = (σlat (σbas − 0.5Wadh))1/3

∂ΔGm n = m ∂x v

(XI.18)

∑ M,M ′

for hexagonal particles in contact with a substrate on their basal face.

Expressions similar to eqs XII.2 and XII.3 hold when derivation is made with respect to y. The equations that allow the determination of the composition (x*, y*) of the critical nuclei are thus the following: • Reciprocal solid solution with y exchange in the framework:

Lamellar minerals, like clay minerals, have a very small aspect ratio which varies during growth in such a way as to nearly conserve the layer thickness. For them, the assumption that all dimensions of the particle obey the same growth law, made in the 3D formalism, is non-valid. For example, in the Wulff− Kaishev equation, which gives the equilibrium shape of an hexagonal particle lying on its basal face (see above), the 2D limit is obtained by simultaneously assuming53 that σbas − 0.5Wadh = 0 and e → 0. The former condition means that the interfacial energy between the particle and the substrate is vanishingly small, which is fulfilled when the particle nucleates and grows on a substrate of same nature (homoepitaxy). It is also necessary that the layer thickness be negligible with respect to other dimensions. In that case, the change in surface Gibbs free energy during nucleation becomes:

vRT ∂ ln F0(x*, y*) ∂x nm = [(1 − y*)ΔμA + y*ΔμC ][(1 − y*)vB + y*vD] − [(1 − y*)ΔμB + y*ΔμD ][(1 − y*)vA + y*vC] vRT ∂ ln F0(x*, y*) nm ∂y = [(1 − x*)ΔμA + x*ΔμB ][(1 − x*)vC + x*vD] − [(1 − x*)ΔμC + x*ΔμD ][(1 − x*)vA + x*vB]

(XI.19) 2

for tetragonal particles (V = nv = el ) or ΔGsurf = 6leσlat

(XI.20)

(XII.4)

• Reciprocal solid solution with y exchange in the interlayer region: The first equation in eq XII.4 is unchanged, while the second one reads:

for hexagonal particles (V = nv = 3el √3/2). As a function of n, ΔGsurf reads: 2

ΔGsurf = Zn1/2v1/2(x , y)σlat(x , y)

⎡ ∂X ∂X ⎤ vM ′ΔμM ⎢XM ′ M − XM M ′ ⎥ ⎣ ∂x ∂x ⎦ (XII.3)

B. 2D Particles

ΔGsurf = 4leσlat

(XII.2)

(XI.21)

ΔμB = ΔμD

1/2

with Z = Nlat(e/s) , Nlat being the number of lateral faces of the particles (Nlat= 4 and 6 for tetragonal or hexagonal particles, respectively), and s being a form factor equal to 1 or (3(√3)/

(XII.5)

• Three end-member ideal solid solution with y exchange in the framework: 3455



Article

I(x, y)

Instantaneous stoichiometric saturation function with respect to a SS of composition x, y Ii(x, y) Interfacial instantaneous stoichiometric saturation function with respect to a SS of composition x, y KM Solubility product of the end-member M K(x, y) Stoichiometric solubility product of a SS with composition x, y kB Boltzmann constant (1.3806504 × 10−23 J/ K) l, l′ Basal lengths of hexagonal and tetragonal particles (m) lm(x, y) Basal size of particles of composition x, y when ΔGN(n, x, y) is maximum with respect to n (m) n Number of growth units in a particle, referred to as its “size” nm(x, y) Number of growth units in a particle of composition x, y when ΔGN(n, x, y) is maximum with respect to n (m) n* = nm(x*, y*) Size of the critical nucleus (growth units) nMs Excess surface quantities of end-member M at the surface of a particle (growth units) NA Avogadro’s number (6.0221418 × 10 23 mol−1) ΔQM(t) Total amount of end-member M (formula units) withdrawn from the AS at time t R Gaz constant (8.314472 J/K/mol) (R = kBNA) S Configurational mixing entropy in the SS (J) T Temperature (K) t Time (s) vM Growth unit volume of end-member M (m3) v(x, y) Volume of one formula unit of a SS with composition x, y (m3) Wadh Adhesion energy of a particle on a substrate (J/m2) XM Molar fraction of the end-member M in the SS x, y Composition of the SS x*, y* Composition of a critical nucleus κ Linear growth constant (m/s) λM Activity coefficient of the end-member M in the SS ΔμM Change of chemical potential of the endmember M (J) ΔμMi Interfacial change of chemical potential of the end-member M (J) ρ Radius of spherical particles (m) σlat(x, y) Surface energy of the lateral faces of SS particles with composition x, y (J/m2) σbas(x, y) Surface energy of the basal faces of SS particles with composition x, y (J/m2) σ̅(x, y) Mean surface energy of SS particles with composition x, y (J/m2)

vRT ∂ ln F0(x*, y*) ∂x nm = ΔμA [(1 − y*)vB + y*vC] − vA[(1 − y*)ΔμB + y*ΔμC ] vRT ∂ ln F0(x*, y*) ∂y nm = [(1 − x*)ΔμA + x*ΔμB ][vC − vB] − [ΔμC − ΔμB ][(1 − x*)vA + x*vB]

(XII.6)

• Three end-member ideal solid solution with y exchange in the interlayer region: the first equation in eq XII.6 is unchanged, while the second equation reads: ΔμB = ΔμC

(XII.7)

The equations that determine the surface excesses and the composition of the deposited layer during growth are those written in the main text.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS We acknowledge support from the CNRS program PEPS 2011 ”Physique Théorique et ses Interfaces”. LIST OF SYMBOLS aM Activity of the end-member M in the SS DÂ , DĈ Effective diffusion constants of the anionic part of a mineral (Â ) or of a cation (Ĉ ) (m2/ s) e Thickness of hexagonal or tetragonal particles (m) F Nucleation frequency (number of nuclei/s/ kg H2O) F0(x, y) Prefactor of the nucleation frequency of a particle with composition x, y (number of nuclei/s/kg H2O) ΔGN(n, x, y) Change of Gibbs free energy during the deposition of n formula units of an SS of composition x, y, on the outer layer of a particle (J) ΔGG(n, x, y) Change in Gibbs free energy during the growth of a particle of size n and composition x, y in its outer layer (J) ΔGm(x, y) Maximum value of ΔGN(n, x, y) with respect to n (J) ΔG* Nucleation barrier, ΔGm(x*, y*), for a critical nuclei of composition (x*, y*) (J) ΔG0M Standard Gibbs free energy of dissolution of the end-member M (J) ΔGexc Excess Gibbs free energy of mixing (J) ΔGsurf Surface contribution to Gibbs free energy (J) IM Saturation state of the AS with respect to the end-member M IMi Interfacial saturation state of the AS with respect to the end-member M

■

REFERENCES

(1) Andrieux, P.; Petit, S. Appl. Clay Sci. 2010, 48, 5−17. (2) Barnes, I.; O’Neil, J. R. Geochim. Cosmochim. Acta 1971, 35 (7), 699−718. (3) Plummer, L. N.; Mackenzie, F. T. Am. J. Sci. 1974, 27 (4), 61−83. (4) Rimstidt, J. D.; Balog, A.; Webb. J. Geochim. Cosmochim. Acta 1998, 62, 1851−1963. 3456



Article

(40) Gaman, A. I.; Napari, I.; Winkler, P. M.; Vehkamäki, H.; Wagner, P. E.; Strey, R.; Viisanen, Y.; Kulmala, M. J. Chem. Phys. 2005, 123, 244502−244512. (41) Markov, I. V. In Crystal Growth for Beginners: Fundamentals of Nucleation, Crystal Growth and Epitaxy; World Scientific, London, 1995. (42) Nielsen, A. F. Crystal Growth; Peiser, H. S., Ed.; Pergamon, Oxford, 1967. (43) Pina, C. M.; Putnis, A. Geochim. Cosmochim. Acta 2002, 66, 185. (44) Zangwill, A. In Physics at Surfaces; Cambridge University Press: Cambridge, U.K., 1990. (45) Sparks, D. L.; Carski, T. H. Appl. Clay Sci. 1985, 1, 89−101. (46) Tertre, E.; Beaucaire, C.; Juery, A.; Ly, J. J. Colloid Interface Sci. 2010, 347, 120−126. (47) Tertre, E.; Prêt, D.; Ferrage, E. J. Colloid Interface Sci. 2011, 353 (1), 248−256. (48) Aziz, M. J. Mater. Sci. Eng 1988, 98, 369−372. (49) Maugis, P.; Gouné, M. Acta Mater. 2005, 53, 3359−3367. (50) Prieto, M.; Fernandez−Gonzalez, A.; Putnis, A.; FernandezDiaz, L. Geochim. Cosmochim. Acta 1997, 61, 3383−3397. (51) Myhr, O. R.; Grong, O. Acta Mater. 2000, 48, 1605−1615. (52) Ratke, L.; Voorhees, P. W. W. In Growth and Coarsening: Ostwald Ripening in Material Processing; Springer: Germany, 2002. (53) Müller, P.; Kern, R. Surf. Sci. 2000, 457, 229−253. (54) Israelachvili, J. N. In Intermolecular and Surface Forces, with Applications to Colloidal and Biological Systems; Academic Press: New York, 1985.

(5) Shao, H.; Dmytrieva, S. V.; Kolditz, O.; Kulik, D.; Pfingsten, W.; Kosakowski, G. Appl. Geochem. 2009, 24, 1287−1300. (6) Wood, J. R. Geochim. Cosmochim. Acta 1975, 39, 1147−1163. (7) Neuhoff, P. S.; Ruhl, L. S. Chem. Geol. 2006, 225, 373−387. (8) Meunier, A. Argiles; Collection Géosciences; Editions Scientifiques Gb: Paris, 2002; ISBN 284703014X. (9) Ganguly, J. In Solid Solutions in Silicate and Oxide systems; Geiger, C. (EMU notes in Mineralogy, Eötvös University Press: Budapest, 2001; Chapter 3, pp 37−71). (10) Jung, D. Y.; Vinograd, N. L.; Fabrichnoya, O. B.; et al. Earth Planet. Sci. Lett. 2010, 295, 477−486. (11) Kulik, D. A.; Vinograd, V. L.; Paulsen, N.; Winkler, B. Phys. Chem. Earth 2010, 35, 217−232. (12) Gale, J. D. J. Chem. Soc., Faraday Trans. 1997, 93, 629−637. (13) Vinograd, V.; Winkler, B. Rev. Mineral. Geochem. 2010, 71, 413− 436. (14) Vinograd, V. L.; Paulsen, N.; Winkler, B.; van de Walle, A. CALPHAD 2010, 34, 113−119. (15) Zhang, D.; Zhou, C.-H.; Lin, C.-X.; Tong, D.-S.; Yu, W.-H. Appl. Clay Sci. 2010, 50, 1−11. (16) Sverjensky, D. A. Geochim. Cosmochim. Acta 1984, 48, 1127− 1134. (17) Sverjensky, D. A. Geochim. Cosmochim. Acta 1985, 49, 853−864. (18) Wang, Y.; Xu, H. Geochim. Cosmochim. Acta 2001, 65 (10), 1529−1543. (19) Tesoriero, A. J.; Pankow, J. F. Geochim. Cosmochim. Acta 1996, 60, 1053−1063. (20) Prieto, M. Rev. Mineral. Geochem. 2009, 70, 47−85. (21) Thorstenson, D. C.; Plummer, L. N. Am. J. Sci. 1977, 277, 1203−1223. (22) Glynn, P. D.; Reardon, E. J. Am. J. Sci. 1990, 290, 164−201. (23) Fritz, B. Sci. Geol. Mém. 1981, 65, 197p. (24) Tardy, Y.; Fritz, B. Clay Miner. 1981, 16, 361−373. (25) Blanc, Ph.; Bieber, A.; Fritz, B.; Duplay, J. Phys. Chem. Miner. 1997, 24, 574−581. (26) Nourtier-Mazauric, E.; Guy, B.; Fritz, B.; Garcia, D.; Clément, A. Oil Gas Sci. Technol. - Rev. IFP 2005, 60 (2), 401−415. (27) Fritz, B.; Richard, L.; McNutt, R. H. Water-Rock InteractionWRI7: Proceedings of the 7th International Symposium; July 13−18, 1992, Park City, Utah; Khakara, Y. K., Maest, A. S., Eds.; A.A. Balkema Publishers: Rotterdam, The Netherlands, 1992; pp 927−930. (28) Richard, L. Introduction des signatures isotopiques dans la modélisation des interactions fluides−roches. Ph.D. Thesis; University Louis Pasteur: Strasbourg, 1993. (29) Madé, B.; Clément, A.; Fritz, B. Comput. Geosci. 1994, 20 (9), 1347−1363. (30) Parkhurst, D. L.; Appello, C. A. J. User’s Guide to PHREEQC (Version 2). A Computer Program for Speciation, Batch-Reaction, OneDimensional Transport, and Inverse Geochemical Calculations; U.S. Department of the Interior: Washington, DC, 1999. (31) Noguera, C.; Fritz, B.; Clément, A.; Baronnet, A. J. J. Cryst. Growth 2006, 297, 180−186. (32) Noguera, C.; Fritz, B.; Clément, A.; Baronnet, A. J. J. Cryst. Growth 2006, 297, 187−198. (33) Fritz, B.; Clément, A.; Amal, Y.; Noguera, C. Geochim. Cosmochim. Acta 2009, 73, 1340−1358. (34) Fritz, B.; Noguera, C. Rev. Mineral. Geochem. 2009, 70, 371− 410. (35) Noguera, C.; Fritz, B.; Clément, A.; Amal, Y. Chem. Geol. 2010, 269, 89−99. (36) Noguera, C.; Fritz, B.; Clément, A. Geochim. Cosmochim. Acta 2011, 75 (12), 3402−3418. (37) Reiss, H.; Shugard, M. J. Chem. Phys. 1976, 65, 5280−5293. (38) Laaksonen, A.; McGraw, R.; Vehkamäki, H. J. Chem. Phys. 1999, 111, 2019−2027. (39) Noppel, M.; Vehkamäkia, H.; Kulmala, M. J. Chem. Phys. 2002, 116, 218−228. 3457


A Theoretical Treatment of the Precipitation of Doubly Substituted

Recommend Documents