Development of Compositional Patterns during the Growth of Solid

Article pubs.acs.org/crystal

Development of Compositional Patterns during the Growth of Solid Solutions from Aqueous Solutions: A Cellular Automaton Simulation Mário A. Gonçalves*,† and Manuel Prieto‡ †

Departamento de Geologia and CeGUL, Universidade de Lisboa, 1749-016 Lisboa, Portugal Departamento de Geología, Universidad de Oviedo, 33005 Oviedo, Spain

‡

ABSTRACT: Crystal growth of solid solutions (SS) in aqueous media frequently leads to the development of oscillatory compositional patterns. Oscillatory zoning (OZ) has recurrently attracted the attention of researchers hypothesizing different generation mechanisms, some of which have been implemented in computer models. Unfortunately, these models typically disregard thermodynamic factors such as the solubility of the SS end-members, which have been shown to be fundamental in the development of compositional patterns. Using a broader approach, this work aims to reproduce the compositional evolution of SS crystals, regardless of whether OZ occurs. We use a cellular automaton model that simulates diffusion and attachment of growth units to a crystal surface, taking (Cd,Ca)CO3 and (Sr,Ba)CO3 as example systems. The model demonstrates the importance of the solubility difference among the SS end-members. When that difference is small, the crystals grow fairly homogeneous in composition; however, they tend to develop a sharp compositional zoning when the end-member solubility products differ by orders of magnitude. A compositiondependent kinetic parameter, the threshold supersaturation, is also shown to strongly influence the crystallization behavior. The simulations are qualitative but reproduce the main features of the experimental patterns and provide a rough estimation of ionpartitioning during crystal growth under non-steady-state conditions.

1. INTRODUCTION The growth of minerals in ambient temperature aqueous systems often leads to the formation of solid solutions, which are mixed crystals where ions of different types substitute for each other in equivalent structural positions. Solid solution− aqueous solution (SS−AS) effects have gained increasing attention over the last few decades, not only because of their relevance in reconstructing past growth environments1 but also because of the role of solid solutions as sequestering phases for harmful ions.2 During growth, the substituting ions typically do not incorporate into the solid phase in the same stoichiometric proportion as in the aqueous phase. Therefore, in finite closed systems, both aqueous and solid compositions tend to change during the growth process, and the crystals develop concentric compositional zoning. The term “concentric” is taken in a broad sense to mean the compositional change measured along the normal to the growing solid front. The system may reach a “partial” equilibrium condition3 in which only the outer layers of the crystals maintain equilibrium with the aqueous solution. True equilibrium would require ongoing dissolution−recrystallization of the zoned crystal to form a solid solution of homogeneous composition at equilibrium with the aqueous solution. However, because the outer layers protect the inner zones from contact with the aqueous solution, the crystals can preserve their zoning for a long time. Concentric zoning can also occur during crystal growth in open transport-reaction © 2014 American Chemical Society

systems, particularly in constricted porous media where masstransfer occurs by diffusion.4,5 Concentric compositional zoning is a record of the reaction pathway undergone by the system during crystal growth; decoding such a record is not an easy task. Predicting the reaction path to be followed from a given initial set of conditions is also difficult. Although the distribution of the substituting ions between the solid and the fluid depends primarily on thermodynamic factors, the effective compositions are also determined by the growth mechanisms and kinetics. A number of models have been developed to account for the effective composition of solid solutions during nucleation6 and growth,7 but all of these models have limitations when applied to actual scenarios. For example, the surface-entrapment model by Watson7 is widely used to account for nonequilibrium partitioning during layer-by-layer growth; unfortunately, as discussed by the author, the number of unknown parameters is a barrier to making quantitative predictions. Moreover, these models deal with either “instantaneous” precipitation events or stationary growth states. However, crystal growth actually occurs under changing conditions which can lead to the Received: January 3, 2014 Revised: April 21, 2014 Published: April 24, 2014 2782

dx.doi.org/10.1021/cg500010p | Cryst. Growth Des. 2014, 14, 2782−2793

Crystal Growth & Design

Article

repetitive spatial pattern because of its own dynamics, even under steady growth conditions. In all of the previously cited literature, the hypothesized mechanisms are implemented in a computer model; if the model obtains oscillatory solutions, it is considered functional without further investigation. Nevertheless, none of the models is definitive or unique in explaining a given oscillatory pattern. In that respect, the sharp OZs observed in (Ba,Sr)SO4 and (Ca,Cd)CO3 crystals grown in an aqueous diffusion medium were found to be strongly influenced by the solubility difference between the end-members in both solid solutions.4,5,13 Both experimental systems attracted some attention from researchers which modeled the oscillatory patterns using different approaches.22−25,32 However, the solubility effect described by Prieto and co-workers5 was mostly disregarded, and some models used strongly supersaturated conditions instead,22 showing that other mechanisms were sought. Thus, despite their success in developing OZ patterns, most models do not describe the workings of these particular systems. The weaknesses of the models become clearer by considering that solid solutions with close end-member solubility products, such as Ba(CrO4,SO4) or (Ba,Sr)CO3, do not develop OZ when grown in an analogous diffusion medium.5 The only attempt to connect the SS−AS thermodynamics was performed by Lubashevsky and co-workers,33 who extended their own initial model25 to accommodate the solubility difference between the end-members. 2.2. The Cellular Automaton Approach. Although computer modeling of OZ has been mostly performed by numerically solving a set of nonlinear partial differential equations (PDEs), there are several advantages in using cellular automata for modeling such nonlinear problems.34 Cellular automata are flexible for implementing models and simulating the intricate dependencies observed in SS−AS reactiondiffusion systems that have been overlooked so far (e.g., coupling between thermodynamics, kinetics, external and intrinsic processes). This work addresses some of these issues by implementing a cellular automaton that builds upon the one previously proposed by L’Heureux and Katsev.24 The cellular automaton is developed to better understand the mechanisms that generate OZ and to study how the different model parameters and boundary conditions affect the growth history in SS−AS systems. This paper critically analyzes the conditions in which the results are obtained and outlines some caveats in cellular automaton modeling. 2.3. Thermodynamic Description of SS−AS Systems. A complete description of the thermodynamic solubility of a binary solid solution (B,C)A requires two mass-action equations, one for each end-member component, BA and CA, and can be described as follows:

development of complex compositional patterns, such as oscillatory zoning. In this work, we use a cellular automaton microscopic model that simulates the diffusion and attachment of growth units to a crystal surface. Both thermodynamic and molecular-scale mechanistic factors are considered when designing the cellular automaton rules, taking the nearly ideal otavite-calcite (CdCO 3 −CaCO 3 ) and strontianite-witherite (SrCO 3 − BaCO3) solid solutions as example systems. The model prediction is qualitative because it is able to reproduce the main features of the compositional patterns observed in crystallization experiments but not a perfectly matching composition. However, we can extract quantitative assessments from the model as well, such as a rough estimation of the effective ion-partitioning during crystal growth under nonsteady-state conditions.

2. BACKGROUND 2.1. Oscillatory Zoning: An Intriguing Case of Concentric Zoning. While a variety of concentric patterns can develop during crystal growth, computer modeling has been mostly focused on oscillatory zoning. That is not surprising; OZ in crystals is among the most striking patterns found in geological structures and appears in a vast number of mineral systems.8 The occurrence of patterns in nature has sparked interest in scientists for a long time, and the variety of patterns displayed is matched by an approximately equal amount of explanatory hypotheses for their formation. The rediscovery of chaos theory in the late 1960s and the subsequent development of the study of nonlinear dynamical systems opened up a wealth of opportunities to explore these systems through numerical modeling.9 In chemical systems, oscillatory patterns often arise as the result of autocatalytic or feedback processes for which significant insight has been gained. Computationally, a landmark was reached when the properties of a class of elementary cellular automata were extensively investigated, showing unexpected complexity behavior and pattern formation using a few simple rules.10 These were fundamental discoveries that have since guided much of the research on complex patterns. Experimentally produced OZ patterns have been obtained in calcite through the incorporation of trace elements11 and in sulfate and carbonate solid solutions.4,5,12−14 These experiments and the observations of the phenomenon in natural samples from different genetic environments15−19 have provided a wealth of hypotheses (often backed up by numerical models) to explain their causes. The hypotheses include (i) constitutional undercooling feedback mechanisms involving diffusion and growth kinetics,20 (ii) surface-compositioncontrolled autocatalytic growth combined with diffusive transport,21−25 (iii) the effect of fluctuations in the external environment during growth,26−28 (iv) positive feedback generated by the interplay of growth-inhibiting cations (e.g., Mn2+, Zn2+, Fe2+) with growth-induced H+ build-up at carbonate-water interfaces,29 (v) variable composition surface relaxation,30 and (vi) micromorphological instability of growing faces.14 OZ can occur by different mechanisms that are primarily classified as either external or intrinsic. In the former case, compositional oscillations result from periodic changes in the bulk growth environment, such as a variable mass flux through an open system or fluctuations in temperature and pressure. In contrast, intrinsic oscillations are self-organization phenomena31 in which a far-from-equilibrium system adopts a

a Bn+a An− = KBAaBA = KBAγBAXBA

(1)

and aCn+a An− = K CAaCA = K CAγCAXCA

(2)

where aBn+, aCn+, and aAn− are the activities of the corresponding ions in aqueous solution, KBA and KCA are the thermodynamic solubility products of the pure end-members, and aBA and aCA are the activities of both components in the solid solution. These activities are related to the mole fractions, XBA and XCA = 1 − XBA, by the solid-phase activity coefficients, γBA and γCA, which for ideal solid solutions become equal to 1. 2783



Article

Because eqs 1 and 2 describe thermodynamic equilibrium, so does any combination of these equations. Using the aqueous activity fraction, XB,aq = aBn+/(aBn+ + aCn+), and eqs 1 and 2 one can obtain the following expression:35 KBAγBAXBA XB,aq = (KBAγBA − K CAγCA )XBA + K CAγCA (3)

which tends to be either Cd- or Ca-rich by modifying XCd,aq by a small amount (some thousandths). Ion partitioning can also be described quantitatively using the distribution coefficient. For the binary (B,C)A solid solution, the equilibrium (eq) distribution coefficient of Bn+ ions is given by

Equation 3 relates the equilibrium aqueous activity fraction (0 ≤ XB,aq ≤ 1) to a given solid solution composition (0 ≤ XBA ≤ 1). If the difference between the solubility products of the endmembers is large (KBA ≪ KCA), a wide range of aqueous solution compositions will be in equilibrium with solids that are rich in the less soluble component (XBA > 0.9). Only when XB,aq is a small fraction (at least the same order of magnitude as the difference between the solubility products, KCA−KBA), the aqueous solution will be in equilibrium with more CA-rich solids. That is the case for the barite-celestite (BaSO4−SrSO4) solid solution in which the end-member solubility products differ by several orders of magnitude (KBrt = 10−9.98 and KClt = 10−6.63). The effect is more pronounced in the otavite-calcite (CdCO3−CaCO3) solid solution (KOtv = 10−12.1 and KCal = 10−8.48). The difference between end-member solubilities induces a strong preferential partitioning of one of the cations into the crystallizing solid solution. Equation 3 can be used to construct an XB,aq−XBA plot, which is known as a “Roozeboom diagram” and is useful for visualizing the equilibrium distribution of the substituting ions in the system. Figure 1 shows the equilibrium XCd,aq−XOtv

=

⎛X DB(eq) = ⎜ BA ⎝ XCA

a Bn+ ⎞ ⎟ aCn+ ⎠

eq

K CA γCA KBA γBA

(4)

This new equilibrium condition can be deduced by combining eqs 1 and 2. Values of DB > 1 indicate preferential partitioning of Bn+ toward the solid phase. On the contrary, for DB < 1, Bn+ tends to incorporate preferentially into the fluid. Note that DC = 1/DB. The quotient γCA/γBA is a correction factor that stands for the nonideality of the solid solution. Distribution coefficients are frequently expressed using total concentrations of the aqueous ions instead of activities. In such a case, eq 4 should include two additional correction factors, γBn+/γCn+ and αBn+/αCn+, where γi is the activity coefficient of the aqueous ion i and αi the uncomplexed fraction of the corresponding ion.40 In the case of the otavite-calcite solid solution DCd(eq) ≈ KCal/ KOtv ≈ 4.17 × 103, and comparable values occur for the partition of barium to the barite-celestite solid solution. During growth, the preferential incorporation of Cd ions (Ba in the barite-celestite example) into the solid solution can promote their rapid depletion in the vicinity of the crystal surface. As a result, the system suddenly switches toward the incorporation of Ca ions (Sr in the second example), which accumulates in the aqueous solution. This sharp transition can give rise to strong compositional gradients and OZ patterns in the solid.5 When the end-member solubility products differ by less than 1 order of magnitude, the crystals tend to grow homogeneously or with slight compositional gradients. This is observed in the case of the strontianite-witherite series, with DSr(eq) ≈ KWth/KStr = 5.13 (KWth = 10−8.56 and KSrt = 10−9.27). Thus, the driving mechanism for the development of OZ in SS−AS systems appears to be strongly anchored in the uneven incorporation of the cations into the precipitating solid solution. 2.4. Supersaturation and Nonequilibrium Distribution Coefficients. While the equilibrium relationships in eqs 1−4 provide a foundation for interpreting SS−AS processes, kinetic and mechanistic factors significantly influence the “effective” crystallization behavior. For all crystal growth scenarios, the fundamental link between thermodynamics and crystallization behavior is supersaturation. Supersaturation is a measure of the departure from equilibrium of a given system, which in the case of a pure solid (BA) is given by the quotient (saturation state Ω) between the ionic activity product (IAP = aBn+ aAn−) in the parent solution and the solubility product (KBA). In binary SS− AS systems, there are two disequilibrium conditions, one for each component, and the corresponding saturation states (ΩBA or ΩCA) represent a species of “partial” supersaturation:

Figure 1. Equilibrium Roozeboom diagram for the (Cd,Ca)CO3 solid solution, which is assumed to be ideal. The inset shows the “vertical” branch of the diagram at a larger scale.

diagram for the (Cd,Ca)CO3−H2O system. Otavite and calcite form a complete solid solution36−38 with nearly ideal (γCA ≈ γBA ≈ 1) mixing properties.37 Recent calculations show that the otavite-calcite SS has a positive free energy of mixing for temperatures below 100 K but that the SS is stable for the full range of Cd/Ca compositions at 200 K.39 Therefore, using an ideal model for this solid solution is a reasonable approach. Due to the low solubility of otavite compared with calcite, the Roozeboom curve approximates two perpendicular lines.5 Therefore, aqueous solutions in equilibrium with calcium-rich solid solutions are extremely poor in Cd2+. Only a narrow range of aqueous compositions (0.000027 ≤ XCd,aq ≤ 0.0022) can coexist in equilibrium with solid compositions in the range of 0.1 ≤ XOtv ≤ 0.9. Small changes in the fluid composition within this range result in dramatic changes in the solid composition,

Ω BA (XBA ) =

a Bn+a An− KBAγBAXBA

(5)

aCn+a An− K CAγCAXCA

(6)

and ΩCA (XCA ) = 2784



Article

composition compared to the previous substrate.46−48 This difference determines the elastic stress generated at the interface so that, in practice, the threshold supersaturation for developing a layer with a given composition will depend on the composition of the previous layer.

In contrast to stoichiometric solids, inspection of eqs 5 and 6 shows that in SS−AS systems supersaturation depends on the solid solution composition. A given aqueous solution can be supersaturated with respect to a range of solid compositions and undersaturated with respect to another range.35 Crystal growth occurs under nonequilibrium conditions, and this can lead to a significant deviation of the “effective” (eff) distribution coefficient from the equilibrium (eq) value, particularly at high growth rates and supersaturation.41,42 In general, there is an attenuation of the degree of preferential partitioning, and the solid phase becomes richer in the most soluble component.6 Thus, when DB(eq) > 1, the attenuation results in DB(eff) < DB(eq) and vice versa. A suitable simulation model should be able to reproduce qualitative growth patterns and the effective ion partitioning observed during growth. 2.5. Nucleation Threshold. According to the classical nucleation theory (CNT), the nucleation rate depends sharply on the supersaturation; there is a critical value (Ω*) that marks the breakdown of the solution metastability. The critical supersaturation is related to the solubility by the general rule that when the solubility of a substance is smaller, the interfacial tension and the critical supersaturation are higher.43 This factor must be considered when interpreting nucleation in SS−AS systems, where the solubility of the end-members can differ by orders of magnitude. For example, the critical supersaturation for the nucleation of barite has been estimated to be ∼25 times higher than that of celestite, and intermediate values can be expected for different members of the (Ba,Sr)SO4 solid solution.4 As a consequence, Sr-rich members can precipitate from aqueous solutions that are more supersaturated with respect to Ba-rich compositions because the critical supersaturation with respect to the more soluble Sr-rich members is lower. The previous ideas hold not only for homogeneous nucleation but also for heterogeneous nucleation and for the transitional supersaturation between spiral and two-dimensional (2D) nucleation mechanisms of crystal growth. In the case of SS−AS systems, the critical supersaturation for 2D nucleation depends on the solid solution composition44 and is lower for the most soluble end-member. When there are changes in composition during growth (with the subsequent defective lattice matching), the process is better described as an epitaxial growth phenomenon which involves a critical supersaturation barrier. However, in the following discussion we are going to use the term “threshold” instead of “critical” supersaturation. In CNT, the concept of critical supersaturation implicitly means that supersaturation is uniform through the system and is reached instantaneously. This concept does not apply to diffusing-reacting systems in which the supersaturation changes rapidly and metastable states are preserved.45 Instead, the threshold supersaturation (Ωτ) is defined as the maximum supersaturation or metastability limit that can be reached under a given set of changing conditions. Similar to most kinetic parameters, Ωτ does not have a unique value, but it instead depends on the system history. Changing the starting conditions modifies the system evolution and thus modifies Ωτ. In the case of SS−AS systems, the threshold supersaturation typically increases from the most to the least soluble end-member,4 which allows us to define a relative scale of metastability. In the simplest approach, a linear variation of Ωτ with the solid solution composition (XBA) is assumed.44 When the composition changes during growth, the development of a new epitaxial layer depends on the difference of the

3. METHODS 3.1. The Cellular Automaton Model. Cellular automata are a class of dynamical computational models that are discrete in space and time and whose behavior is completely specified by a set of rules that govern the relationships between constituents. These models were first proposed by the mathematician John von Neumann in the 1940s, but it was not until the seminal paper of Wolfram10 that cellular automata were definitely legitimized as a useful research tool in physics. One of the most common applications of cellular automata is qualitative discrete modeling, including their operational use as an alternative to PDEs. Models of reaction-diffusion systems fall into this category.49 For x > 0 representing the aqueous solution domain and ci(x,t) representing the molar concentration of species i (i = 1, ..., n) in aqueous solution at position x and time t, the one-dimensional (1D) growth of a crystal surface (with the frame of reference x = 0 at the crystal growing front) can be described by the diffusion-advection equation: ∂ci ∂ 2c ∂c = Di 2i + v i ∂t ∂x ∂x

(7)

where Di is the effective diffusion coefficient of species i in aqueous solution and v is the linear crystal growth rate. Equation 7 can be integrated across a finite boundary layer (adopting constant bulk concentrations beyond this layer); using the “crystal growing front” boundary condition results in a set of n+1 approximate ordinary differential equations (ODEs): one for each ci(x,t) and one for X, which is the mole fraction of the corresponding component in the solid solution. Therefore, if we impose the boundary layer approximation, the system of PDEs can be replaced by a simpler system of ODEs.22 Such a macroscopic model disregards the molecular processes that take place at specific locations on the interface, which can be important in modeling three-dimensional (3D) compositional patterns. Advantageously, cellular automata can be used to simulate microscopic interactions, and their macroscopic behavior can be described by a similar set of ODEs.24 The cellular automaton developed for this work models the experimental arrangement studied by Prieto and co-workers,5 taking the “ideal” otavite-calcite solid solution, (Cd,Ca)CO3 as an example. This solid solution was studied in laboratory experiments conducted in a U-shaped tube (Figure 2), in which each of the vertical branches was filled with either CaCl2 + CdCl2 or Na2CO3 solutions. The dissolved ions counter-diffuse from these solution reservoirs through a porous silica hydrogel (∼95.6 wt % water in interconnecting pores) that fills the horizontal branch of the tube. Because gel properties inhibit

Figure 2. Scheme of the crystallization experiments. A U-shaped tube is filled with silica hydrogel in its horizontal branch and with CdCl2 + CaCl2 and Na2CO3 solutions in the left and right vertical branches, respectively. The solutes are left to diffuse along the hydrogel, as shown by the arrows, and crystals start to form approximately midway between both reservoirs. 2785



Article

advection and convection, the crystallization medium is purely diffusive. Porous media are effective at suppressing nucleation,45 and this experimental design therefore maintains high levels of supersaturation during crystal growth. Crystals tend to nucleate about midway in the gel column, which is the place where supersaturation builds up faster because the counter-diffusing ions meet at this point. Several crystals were obtained showing concentric zoning; the patterning differs depending on the Ca:Cd ratio and the concentration used in the parent solution. The cellular automaton does not consider the nucleation stage but instead assumes that a crystal nucleus greater than the critical size is already present. Thus, the model focuses on reproducing the experimental patterns and defining the controlling variables. The simulation is also applied to the witherite-strontianite (Ba,Sr)CO3 system to check the effect of the solubility difference on the compositional patterning. 3.2. Definitions, Topology, Boundary and Initial Conditions. Cellular automata are models composed of m lattice sites with a set of i (i = 1, ..., m) states, si, each representing the state of a lattice site. The sites are arranged in several possible configurations. A lattice site may be either empty (si = 0) or occupied (si ≠ 0), in which case the state value, si, can also represent the identity of the particle occupying the lattice site or any other property of that particle that is defined by a unique integer value. For the sake of clarity, the term “cell” is used to refer to any lattice site in the cellular automaton, either occupied or empty, and the term “particle” refers to any state, sj ≠ 0, of the j cell belonging to the subset of the p occupied cells (j = 1, ..., p). Associated with each cellular automaton model there is also a set of k (k = 1, ..., l) rules, Rk, that determine transitions of state in each lattice site as a function of its neighbor sites in the grid. Models can represent 1D to 3D grids with different symmetries. In the current work, the cellular automaton model has a 2D square symmetry (square cells), although not necessarily with a square domain. Each step in the iteration process simulates time’s passage, in which a uniformly distributed random sequence of p numbers, Mr ∈{1,2,3,4}, r = 1, ..., p, is generated. Those numbers determine the direction of movement of each of the p particles toward each of the four possible directions. Particles move into a new position according to the set of rules, Rk, determining the configuration that the model will have at that particular time. The details of these rules depend on whether the particle is freely moving in the solution or whether it is in the neighborhood of the crystal surface, and will be explained later. Because of the geometry of the experimental setup used in the crystallization experiments, the cellular automaton considers a rectangular 2D geometrical domain folded into a cylinder sheet (Figure 3). According to this geometry, the following boundary conditions apply: (a) Input of dissolved ions occurs through boundary A. In each iteration step, whenever any particle either moves outside the solution domain through A or becomes attached to the crystal surface (at C, see below), new particles are incorporated through A. Newly input particles are randomly generated according to the atom fraction in the parent (bulk) solution, X[Cd] = Cd/(Ca + Cd), such that the total number of particles inside the domain at any time frame, t, is conserved. In contrast to the aqueous activity fraction defined in section 2.3, X[Cd] involves total concentrations (free ions + complexes) of the substituting atoms. (b) Boundaries B and B′ are periodic, and particles are free to move in and out such that whenever a particle moves through boundary B it reappears in the same relative position on boundary B′, and vice versa. This boundary condition defines the cylindrical topology of the model. (c) In boundary C, only crystal precipitation is possible according to the precipitation rules. A variant of this model was also developed in which this boundary condition was slightly modified, acting as an open boundary where particles could move out with a very small probability; this would represent a scenario in which a fraction of the ions do not stay near the crystal surface indefinitely and eventually pass beyond the area where the crystal is growing. The cellular automaton has one great single domain, whose area (volume) is maintained constant, representing the “solution”. The “crystal” grows at boundary C, and whenever a full layer of the crystal

Figure 3. Geometry of the cellular automaton domain, showing the boundaries (A, B, B′, and C) and the nature of each boundary (A: open; source; B and B′: periodic; conservative; and C: closed; sink). The nature of the B (B′) boundary gives the domain a cylindrical topology. is formed, an additional layer of “solution” is appended at boundary A in order to maintain the solution area (i.e., the solution domain thickness) constant during all the simulation. The solution domain used in most simulations has 300 × 30 cells, but some simulations also used 300 × 60 or 300 × 90 cells to study the influence of the domain size on the results. The cellular automaton has the option to start with different initial conditions with respect to the composition inside the column. The particles inside the solution domain can either (a) be uniformly distributed along the solution domain or (b) be unevenly distributed along the same domain following a linear gradient with decreasing concentration from boundary A to boundary C. In all cases, the maximum number of particles that are allowed to be inside the solution domain is limited by the constant bulk concentration (i.e., the fixed concentration defined outside the boundary A). Therefore, in either of the previous options, the simulation can start with the maximum number of particles or with a fraction of that number. In the latter case, particles are introduced at a rate greater than they are consumed until the maximum number is reached during the simulation. 3.3. Rules. The model has two fundamental sets of rules: a set of movement rules and a set of precipitation rules. The symmetry of the domain defines four possible directions in which each particle can move, each one assigned an integer value Mr. A particle moves randomly to one of the four neighboring empty cells with the same probability, p, using a sequence of pseudorandom numbers generated from a set of n (n = 4) integer values, ai representing each one of the possible directions of motion:

p(Mr = ai) = 1/n , for i = 1, ..., n The movement of the particles is asynchronous, which means that if a particle is about to move to a recently occupied cell by a previous particle in the same iteration, then a new integer value is randomly assigned to that particle among the remaining possible integer values defining the allowed directions in which the particle can move. The rules also change slightly for particles moving in the neighborhood of the growing crystal face (typically one cell away from a surface site) and are described below. First, it is necessary to define how the crystal surface grows and the conditions for a particle to bind to the crystal surface. As it will be discussed later in section 3.4 in more detail, for the sake of simplicity the model ignores the physical presence of the anionic species (although they are considered for the calculation of the ion activity product). In practice, the procedure is equivalent to use CA or BA growth units according to the end-member stoichiometry. A particle 2786



Article

Figure 4. Configurations of a particle (light-gray cell) near the crystal surface. In a, the particle is only a layer away from a surface site and can move to any of the neighboring cells with equal probability. In b−g, the arrows indicate the highest probable movement of the particle among the possible options (in b, the most probable movements toward the left or the right have equal probabilities). In h, the particle rests over a precipitating site. (or growth unit) can bind (sorb) to a surface site whenever the following conditions are met: (a) A surface site (an empty cell sharing at least one boundary with the crystal surface) can bind a particle from solution if it shares two boundaries (kink) or three boundaries (pit) with the crystal surface. When it shares a single boundary (terrace), it can only bind a particle if its other two neighboring empty cells in the same layer are also surface sites. This implies the existence of a minimum of three consecutive structural units in the surface of the crystal to bind a particle from the solution onto the central surface particle; it also implies that particles never bind on top of a kink site. (b) Only a maximum of two growing surface layers (incomplete layers) of the crystal are admissible at a time. Therefore, surface sites satisfying condition “a” in the front layer become active only when the inner layer is filled. These conditions prevent the crystal surface from growing rough, which avoids the growth of a high surface energy crystal and makes implicit the idea that the crystals grow mostly by step advancement or 2D nucleation. Therefore, particles can stand on a surface site that is not a binding site, or they can share boundaries with kink sites and pits. It is in these cases that the movement rules change, by favoring particle movement toward favorable binding sites. Figure 4 shows seven configurations (b−h) in which particles are near binding sites. The arrows indicate the most probable movement the particles will take in the next step. In b, c, and d, one of the directions is forbidden to the particle, but since they all have the same probability to be chosen, the probability assigned to the forbidden one is equally distributed to the direction(s) indicated by the arrows in the figure, keeping the probability for the other directions unchanged. In e, f, and g, all directions are possible, but because the particle is over a kink or a pit the probability to choose the direction toward those binding sites increase. The last configuration shows a particle at a binding site. In such a case, an incorporation probability (precipitation; see below) is assigned instead. If the particle is not incorporated in the crystal, then it moves from that site (desorbs) as the usual movement rules determine. Once a particle from the solution is at a binding site, it can be incorporated into the crystal, and a small set of rules define a simple surface diffusional process before the effective precipitation occurs (Figure 5). The probability of incorporating a particle into the crystal structure is a function of the neighboring solution composition driven by the partition thermodynamics, which for an ideal solid solution (γCA = γBA = 1) can be described by (eq 3):

XBA =

Figure 5. Surface diffusion. Once the particle is in a precipitating site and attaches to the surface, it moves toward a kink site if it is one cell away from that site. In b the particle can move toward the left or the right with equal probability. In c the particle does not move and is fixed there. This rule favors the growth of the surface by step advancement. XBA. However, if B and C ions exist in the neighboring solution but the concentrations are below the corresponding threshold supersaturation (see below), the incorporation probability of BA or CA units into the solid is zero. 3.4. Concentration and Threshold Supersaturation. Because the precipitation rules of the cellular automaton use eq 8, we need to know XCd,aq and hence the activities and concentrations of Cd2+ and Ca2+ in the aqueous solution. The first step is to define a relationship between the particle number and the chemical concentration. Kier and co-workers50,51 studied several physical properties of water using cellular automata and were able to simulate these properties with an automaton in which water molecules occupied 69% of the cells. In other words, the relative occupancy of the total “water volume” by water molecules would be 0.69. The rest (0.31) would be empty space, which could be partially occupied by solutes. With this occupancy and adding solutes to the empty cells, Kier and co-workers simulated various phenomena in aqueous solutions (including diffusion of solutes) that are consistent with experimental results.52,53 By using an occupancy factor of 0.69 for water molecules in the automaton and considering the molar volume of water, we can provide a relationship between the number of particles of each ion in the automaton domain and its concentration. The addition of the solute normally increases the solution volume, but because the concentrations used here are small, the volume change is assumed to be negligible. Particles in the cellular automaton represent the Cd2+ and Ca2+ ions in solution only. For the sake of simplicity, the model does not consider the physical presence of carbonate ions, but it instead assumes that the concentration of CO32− is equal to the sum of the concentrations of Ca2+ and Cd2+. This assumption is supported by the experimental evidence that shows that in counter-diffusion systems, crystallization occurs in a region of the diffusion column where the concentrations of the reacting cations and anions are similar. This “equality range” condition obeys kinetics54 and has been experimentally checked for a number of sparingly soluble carbonates and sulfates.55 From the concentration data, the activities of the ions in solution and the supersaturation are calculated using the Debye−Hückel model for electrolyte solutions. However, the calculations are not simple because the diffusion column contains a high concentration of

K CAXB,aq (K CA − KBA )XB,aq + KBA

(8)

For the configuration of 300 × 30 cells, the composition of the “neighboring” solution is computed for the 50 columns closest to the crystal surface front. Then, the value of XBA calculated from eq 8 is used as the incorporation probability because 0 ≤ XBA ≤ 1. Similarly, the incorporation probability of a CA unit, pCA, is pCA = 1 − pBA = 1 − 2787



Article

Figure 6. Simulations carried out with Ca:Cd ratios of 1:3 (a), 1:1 (b), and 2:1 (c). In each case, the upper image is a snapshot of the automaton with the solution domain (in blue) on the left and the crystal on the right. The yellow and red stripes correspond to Ca- and Cd-rich zones, respectively. The lower graphic shows the corresponding compositional profile (XOtv) of the crystal from core (at 0) to rim (d = distance in cells). The BSE images correspond to central sections of solid solution crystals grown in gels under conditions similar to the simulation parameters, where white and dark areas show Cd-rich and Ca-rich domains, respectively. The data points show some experimental microanalysis (XOtv). The solution bulk for the experiment shown in each subplot is (a) 0.5 mM Ca and 1.5 mM Cd; (b) 1 mM Ca and 1 mM Cd; (c) 2 mM Ca and 1 mM Cd (BSE images by A. Fernández-González, used by permission). dissolved NaCl. In the experiments by Prieto and co-workers,5 the silica hydrogel was obtained by the acidification of a Na2SiO3 solution with 1 N HCl. Under these conditions, NaCl forms as a soluble byproduct (0.2−0.4 mol/kg water, depending on the gel density). Sodium and chlorine ions also get into the diffusion column from the parent solutions, but that is a minor contribution when compared to the “environmental” NaCl. The concentration of Na+ is only used in calculating the activity of the components in solution through its effect on the ionic strength. However, the high concentration of Cl− produces an additional effect that cannot be disregarded. In the presence of 0.2 m NaCl, most dissolved Cd occurs as cadmiumchlorine complexes (CdCl+, CdCl20, CdCl3−, and CdOHCl). The complexation (see section 2.3) factors, α Ca2+ and αCd 2+, are approximately 0.97 and 0.1, respectively. Such an effect implies a much higher value of the {Ca2+}:{Cd2+} activity ratio compared to the

analytical Ca:Cd ratio, which has also been considered to calculate the supersaturation. As discussed in section 2.5, nucleation, 2D growth, and epitaxial growth in nonhomogeneous, evolutionary systems require a threshold supersaturation, Ωτ, to be overpassed. This threshold supersaturation can vary substantially depending on the supersaturation rate45 and is usually higher for less soluble compounds. For example, from 3D nucleation experiments Putnis and co-workers45 obtained values ranging from 1.6 to 5.5 for gypsum, 90 to 665 for witherite, 446 to 3640 for strontianite, and 1800 to 10930 for Barite. In SS−AS systems, the concept of threshold supersaturation allows defining a relative scale of metastability as a function of the crystal surface composition. The present model, developed for the ideal (Cd,Ca)CO3 solid solution, assumes a simple linear variation of Ωτ with composition. Let nsT be the total number of particles (growth units) in the exposed crystal 2788



Article

surface. The number of those particles that are a BA or CA unit is CA τ τ given by nBA s or ns , respectively. Let ΩBA and ΩCA be the threshold supersaturation for the precipitation of pure BA (e.g., otavite) and CA (e.g., calcite), respectively. Then, the threshold supersaturation as a function of surface composition is given by

Ωτ (nsBA ) =

nsBA nsT

⎛ n BA ⎞ + Ω τBA ⎜1 − s T ⎟ ns ⎠ ⎝

(9)

⎛ n CA ⎞ τ ⎜1 − s T ⎟ + ΩCA ns ⎠ ⎝

(10)

4.1. Single-Open Boundary Model. In this model, boundary C is closed and only allows for the precipitation of the solid solution phase. Several simulations were run with different initial Ca:Cd ratios in aqueous solution. The results obtained depend not only on the concentration of the cations in solution but also on the saturation state with respect to both end-members. Figure 6 shows a set of results that formed oscillatory patterns starting from Ca:Cd ratios of 1:3, 1:1, and 2:1. In each of these simulations, the core of the crystal domain is always Cd-rich. As the system evolves, the precipitation of a Ca-rich layer follows, and the system progresses in an alternating but nondeterministic (in contrast with analytical deterministic solutions, the oscillations differ from each other) succession of Cd-rich and Ca-rich layers, whose width depends primarily on the Ca:Cd ratio in the parent solution. The result arises from the interplay between the accumulation of particles in the vicinity of the interface and the rate of precipitation, aided by the threshold supersaturation. Because Cd particles partition strongly toward the solid (which is assumed to start as a pure Cd nucleus seed), they are removed from the neighboring solution at a greater rate than their replenishment by diffusion. Calcium particles only precipitate when their concentration is above the threshold supersaturation for calcite (over an otavite-rich surface). Normally, when that happens, the neighboring solution is extremely depleted in Cd (or even with no Cd at all) so that the activity fraction of Ca can be 1 or close to 1. Under those conditions, there is a massive precipitation of Ca particles over the surface. Indeed, as the surface gets Ca-rich, the threshold supersaturation for calcite decreases (eq 10), which favors continual Ca precipitation. Thus, the threshold supersaturation rule results in an autocatalytic effect. Eventually, Ca also gets depleted in the vicinity of the interface, and the Ca-rich growth period stops. The process is then repeated all over again. Figure 6 shows that the widths of the Cd-rich layers depend on the Ca:Cd ratio in the parent solution. At high Ca:Cd ratios, the Ca-rich layers are never as pure as the Cd-rich ones, and they tend to be wider. As to the Cd-rich ones, they tend to become narrower. These patterns are very similar to those obtained experimentally, as shown in the accompanying backscattered electron (BSE) images. The data points displayed in the images correspond to electron probe microanalyses (EPM) and are also in good agreement with the calculation results. The simulations generally show that the development of the oscillatory pattern depended critically on how close to the threshold supersaturation of calcite the system is maintained. Ideally, the concentration of Ca should move above and below that threshold, so that when the concentration remains systematically above the threshold supersaturation, the oscillatory pattern is gradually lost. Figure 7 shows the result of a simulation with a Ca:Cd ratio of 1:1 using the same conditions as in Figure 6b but doubling the concentrations. Again, this result is consistent with the experimental observations. However, the model fails to produce oscillatory patterns for the parent solutions with low concentrations and Ca:Cd ratios greater than 2. This failure results from a limitation of the single-open boundary model in reproducing the real experiments. In the single-open version, the particles never move out of the solution domain unless they incorporate into the solid. When the bulk concentration of Ca in the parent solution is below the threshold supersaturation for calcite, the Ca particles that accumulate in the boundary layer continue in solution indefinitely. Because the particle that substitutes for

and

Ωτ (nsCA ) =

nsCA nsT

Therefore, the threshold supersaturation for component i (BA or CA) increases linearly as the surface composition gets richer in the other solid solution component, reaching a maximum Ωτi when nis = 0. The threshold supersaturation is 1 (for each end-member) when the exposed surface is a pure phase of the same composition. The difference in the threshold supersaturations used in our simulations ranged from 1 to 2 orders of magnitude, similar to the difference measured by Prieto and co-workers4 in the case of the (Ba,Sr)SO4 system. Values of 900 for ΩτBA (otavite) and 30 for ΩτCA (calcite) seem reasonable for crystal growth and epitaxy, taking into account the difference of solubility between both end-members. Anyway, we notice that the relevant parameter is not the absolute value but the relative difference between these two extreme values. It is worth noting that our cellular automaton does not deal with actual crystal growth rates. That is also the case in other cellular automata24 developed to account for oscillatory compositional patterns. However, the resulting compositional zoning is a kinetic effect. The incorporation of growth units on a layer of different composition requires a certain supersaturation threshold to be exceeded, which implicitly implies a kinetic control of the process by 2D or 3D surface nucleation. Such a threshold is set to 1 when BA grows on a BA pure surface, and the same occurs for the CA endmember. A pure (stoichiometric) crystal can grow by a spiral growth mechanism, which does not require a supersaturation threshold to be exceeded. Under spiral growth conditions, the kinetics depends, among other factors, on both the frequency of arrival of growth units to the crystal surface and the incorporation probability (eq 8) of the corresponding (BA or CA) component. Both factors are intrinsic to the simulation but are not decoded into growth rates. The cellular automaton rules does not work with the classical, macroscopic rate laws that in the case of stoichiometric solids relate growth rate to supersaturation. Dealing with solid solutions, kinetic modeling must relate growth rate and crystal composition to aqueous composition and supersaturation, which is by far more complex. The present model simply aims to describe the compositional evolution during crystal growth of binary solid solutions from aqueous solutions under nonsteady-state conditions.

4. RESULTS AND DISCUSSION Simulations were run using two different versions of the cellular automaton model, which are described as single-open and dualopen boundary versions. As noted in section 3.2, the boundary C of the cellular automaton (Figure 3) can act as either a closed boundary where only crystal precipitation is permitted or an open boundary where in addition to crystal precipitation, a particle can be assigned a nonzero probability that it will move out of the model domain. Because boundary A is always an open boundary, the terms “single” and “dual” designate each model version. The model of L’Heureux and Katsev24 is similar to the single-open boundary version with respect to the boundary conditions, but it is fundamentally different with respect to the precipitation rules, which are the major cornerstones of our model. 2789



Article

Figure 8. Simulation carried out with a Ca:Cd ratio of 1:1 (1 mM Ca and 1 mM Cd). Only a small area (in blue) of the solution domain is shown. There are significant differences in both the automaton image and the compositional profile with respect to the equivalent simulation in Figure 6b.

Figure 7. Simulation conducted using a Ca:Cd ratio of 1:1, as in Figure 6b, but using double Ca and Cd bulk concentrations (2 mM Ca and 2 mM Cd). Only a small area of the solution domain (in blue) is shown. Under these conditions, the solution was kept above the supersaturation threshold for calcite during the simulation. Although the crystal still exhibits an incipient oscillatory zoning, the composition profile is quite irregular with maxima and minima that are less pronounced than are those in Figure 6b.

around that boundary. Instead, a simple probabilistic rule applies to particles when they are in a precipitating position. As presented in section 3, a particle incorporates the solid phase with a given probability p ∈ [0,1] when it reaches a favorable surface site. If the particle is not incorporated into the solid, then it returns back into solution in the single-open boundary model. In the dual-open version, there is also a small constant probability that the particle is simply withdrawn from the system and replaced by a new particle at boundary A. At a molecular scale, the size of a flat growing surface is almost infinite in comparison with the size of the cations and anions in solution. Therefore, this probability must necessarily be very small. Simulations using probabilities >0.005 virtually inhibit the development of sizable Ca-rich layers, unless the Ca fraction in solution is greater than 0.75. The most realistic simulations are obtained using a probability of 0.001 for a particle to move out of the system if it is not incorporated into the solid phase. This value was chosen after a systematic investigation and is more or less halfway from the upper limit mentioned earlier (0.005) and a lower limit (0.0005) in which the results were nearly identical to the single-open boundary model. Therefore, C becomes an open boundary in this version of the cellular automaton. Changing the boundary condition at C produces some modifications in the obtained patterns. Figures 9 and 10 show two simulations that were performed using Ca:Cd ratios of 1:1 and 5:1, respectively. The Cd-rich layers increase their relative width for the same Ca:Cd ratio in solution (Figure 9). Moreover, in contrast to the single-open version, we can obtain oscillatory patterns using Ca:Cd ratios > 2:1 (Figure 10). This important improvement shows how sensitive these systems are to changing (or adding) critical variables. 4.3. Dimension of the System. The final question addressed with these simulations was the influence that the dimension of the system could have on the final results, particularly whether increasing the length of the crystal growing front (boundary C) would bear any influence on the patterns obtained. The results of analogous simulations with solution domains that are 300 × 30, 300 × 60, and 300 × 90 do not show significant variations, although the width and peak compositions of the Cd- and Ca-rich bands may vary,

one that is incorporated in the solid can be either a Ca or a Cd with a probability equal to their fraction in the solution bulk, with time, the system will get locked into a state in which only Ca particles remain in the whole domain of the solution and the system stops evolving. This is not a realistic result; while surface crystal fronts may represent a barrier at the molecular scale, they do not impede the particles from moving forward along the diffusion column. The simulations shown in Figures 6 and 7 started with the maximum number of particles; i.e., each initial concentration in the solution domain was homogeneous and equal to the bulk concentration (see the last paragraph in section 3.1). Starting with fewer particles (lower concentration) or with concentration gradients in the solution domain favors the deposition of a wider initial Cd-rich zone that grows at a lower rate, but it has no other major influences on the final outcome. However, the presence of these gradients is similar to the gel experiments, where the reacting ions diffused along the column and gradually increased their concentration in the crystallization zone, and the width of the Cd-rich crystal core was generally wider than the rest of the compositional bands (see the BSE images in Figure 6). Figure 8 illustrates another interesting effect that is a consequence of the nondeterministic character of the cellular automaton model. The simulation in Figure 8 was performed starting from the same conditions as the simulation whose results are shown in Figure 6b (Ca:Cd = 1:1). The general aspects of both zoning patterns are quite similar, but the details differ significantly. This phenomenon occurs in crystal growth experiments, in which there is typically a lack of correlation between the specific zoning patterns of nearby crystals of the same precipitate.13 4.2. Dual-Open Boundary Model. Allowing the particles to move out of the system at boundary C provides a more realistic model in which particles are not permanently trapped near the crystal surface but can instead move around the crystal border and proceed along the column. Because boundary C is defined by the growing crystal surface alone, the cellular automaton cannot use the movement rules to move a particle 2790



Article

although the range of compositions shifts toward the most concentrated component. 4.5. Compositional Patterns and Effective Ion Partitioning during Growth under Non-Steady-State Conditions. The simulations presented in the previous sections show that the development of oscillatory patterns in a solidsolution crystal results from the interplay between the difference of the solubilities of the end-members and the threshold supersaturation. As previously explained, the supersaturation threshold was originally defined as the limit to be exceeded for nucleation (3D, 2D, epitaxial) to occur in systems with a changing supersaturation. The supersaturation threshold depends primarily on solubility in the same way as the critical supersaturation43 in the classical nucleation theory: the lower the solubility the higher the interfacial tension and the higher the critical supersaturation. However, in the case of crystal growth of solid solutions, the threshold supersaturation depends also on the compositional mismatch between subsequent layers. The concept of threshold supersaturation allows the generation of realistic patterns, as shown for the case of the calcite-otavite system. The oscillatory patterns become progressively degraded as the supersaturation increases (e.g., Figure 7), particularly when the solution composition is dominated by the more soluble component (i.e., at high Ca:Cd ratios, Figures 6c and 10). The development of an oscillatory sequence is especially sensitive to oscillations around the threshold supersaturation of calcite, which is in the key mechanism to switch the system from a Ca-rich to a Cd-rich growth period. These observations agree with the results obtained in several types of experiments. For example, Katsikopoulos58 performed precipitation experiments in strongly agitated solutions at very high supersaturation with respect to the whole range of the calcite-otavite series. The obtained precipitates were homogeneous, with compositions determined by the initial concentration of the parent solutions. All of the experimental compositions indicate that the effective distribution coefficients of Cd are significantly lower at high supersaturation. Katsikopoulos58 measured DCd(eff) values ranging between 7.0 and 31.2, with an average value of 11.5, which is well below the equilibrium value DCd(eq) ≈ 4.17 × 103 (sections 2.3 and 2.4). This author also conducted gel experiments42,58 using diffusion columns of different lengths (8, 18, and 28 cm). A shorter diffusion column results in a higher supersaturation rate (i.e., a faster increase of the supersaturation with time45) in the crystallization zone and a lower DCd(eff). Congruently, the average DCd(eff) values at nucleation were determined to be 29.5, 41.7, and 72.4 for the short, intermediate, and long columns, respectively. The crystals grown using intermediate and large columns show a clear concentric zoning, but the compositional gradients are much less pronounced when the crystals were grown using short columns. The results obtained with the cellular automaton are consistent with these observations. At high supersaturation, the concentric zoning vanishes and the crystal composition reflects only slight variations around the average solution composition. The cellular automaton also allows us to estimate effective distribution coefficients at specified times (iterations) in each simulation. Solution compositions can be obtained from the particles present in a neighboring layer near the crystal surface, and the solid compositions must consider only the outermost layers. In this way, the evolving DCd(eff) can be tracked sequentially. The estimations range from 30 to 422,

Figure 9. Simulation carried out with a Ca:Cd ratio of 1:1 (1 mM Ca and 1 mM Cd) using a dual-open model. The pattern is similar to that in Figure 6b but with more prominent Cd-rich layers and narrower Ca-rich layers.

Figure 10. Simulation performed with a Ca:Cd ratio of 5:1 (2.5 mM Ca and 0.5 mM Cd) using a dual-open model. Despite the high Ca:Cd ratio, oscillatory zoning develops. The patterning is clearly dominated by Ca-rich domains with narrow, well-defined, Cd-rich layers approximately 50−100 units apart from each other (see compositional profile).

particularly when the Ca:Cd ratio deviates significantly from unity. 4.4. Solid Solutions with Close End-Member Solubility Products. As explained in section 3, similar simulations have also been performed using the witherite-strontianite (Ba,Sr)CO3 series to evaluate the effect of the solubility difference on the resulting compositional patterns. There appears to be complete miscibility between the end-members at all temperatures and pressures.56 Moreover, the variation of the unit cell parameters with composition is virtually linear,57 so using an ideal model for the (Ba,Sr)CO3 solid solution is a reasonable approach. Therefore, the probability of incorporating a particle into the crystal structure has been calculated using eq 8 as in the (Cd,Ca)CO3 example. However, in this case, the solubility products (KWth = 10−8.56 and KStr = 10−9.27) differ by less than 1 order of magnitude, which has a significant implication for the results. Figure 11 shows a simulation carried out using the dualopen version with a Ba:Sr ratio of 1:1, i.e., using conditions that give a sharp concentric pattern in the otavite-calcite system. As can be observed, there are some variations (within the 0.4−0.6 range) in the solid composition, but the results are in agreement with the experimental observations (see BSE image) in which no zoning pattern arises. The results are similar starting from different Ba:Sr ratios in the bulk solution, 2791



Article

Figure 11. Growth of a (Ba,Sr)CO3 solid solution crystal. Simulation performed using a Ba:Sr ratio of 1:1 (1 mM Ba and 1 mM Sr) and a solution domain of 300 × 60. There are some variations in the solid composition, but the results are consistent with the experimental observations: no zoning pattern arises. The black and white picture shows the BSE image of the central section of a (Ba,Sr)CO3 crystal grown in silica hydrogel (Ba:Sr = 1:1). The lower graphic shows the corresponding compositional profile (XStr) of the crystal from core (at 0) to rim (d = distance in cells) (BSE image by A. Fernández-González, used by permission).

given growth unit (BA or CA) depends not only on the equilibrium partition coefficient but also on the corresponding threshold supersaturation. Such a different incorporation probability results in effective distribution coefficients that differ from the equilibrium value. With both parameters implemented in the automaton rules, oscillatory patterns arise, particularly when the system is maintained close to the threshold supersaturation for the most soluble end-member. The most sophisticated (dual-open) version of the model, which allows particles to move out of the system, provides more realistic simulations for high Ca:Cd ratios. Finally, the calculated effective distribution coefficients are consistent with the experimental values; this indicates that the proposed cellular automaton can predict reaction paths in the studied systems. This type of simulation could be immediately extendable to other SS−AS crystallization systems, as long as we are able to define the boundary conditions and the movement and precipitation rules suitable for the involved crystal growth environment. To address our main goal, we carefully chose two solid solution systems which are ideal or nearly ideal. Thus, the interplay between the supersaturation threshold and endmember solubility difference was demonstrated to be fundamental to produce complex patterns similar to the experimental observations. However, nonideality effects can also explain the generation of instabilities during the growth process.33 Small deviations from ideality can be easily accommodated in the cellular automaton rules, particularly when the solid solution is complete, but nonideality effects usually imply the presence of miscibility gaps and other complications (ordering, ordering and unmixing, etc.). Modeling the growth behavior in these types of systems is a challenge that remains beyond the scope of the present work. In our opinion, cellular automata can be an excellent tool to progress along this way, but achieving meaningful models still requires a significant amount of work.

which is reasonable because all of them fall between the equilibrium value and the measurements obtained in the gel experiments (note that these measurements correspond to the nucleation event, which occurs at a higher supersaturation than the subsequent growth process). The cellular automaton cannot compete in precision with the macroscopic, kinetic models developed59,60 to determine effective distribution coefficients under steady state conditions, at a constant bulk supersaturation. However, the method can approximately monitor the effective partitioning for non-steadystate scenarios like those involved in the gel experiments (at least during the 3 or 4 months that typically a gel experiment lasts) or in open reaction-transport systems where the supersaturation is continuously changing. To our knowledge there are no models available in the literature that deal with these changing conditions, and the cellular automaton rules could be improved to be an acceptable tool for that purpose.

5. CONCLUSIONS Cellular automata have been shown to be a useful tool to simulate the crystallization behavior of solid solutions from aqueous solutions in reaction-diffusion systems. The proposed version evolved from the former model by L’Heureux and Katsev24 and demonstrates the importance of considering the difference of solubilities between the end-members of the solid solution. When that difference is small, the crystals grow homogeneously in composition; however, they tend to develop a sharp compositional zoning when the end-member solubility products differ by more than 1 order of magnitude. However, a large difference alone is insufficient to simulate the experimental observations, particularly the development of oscillatory patterns. An additional parameter, the threshold supersaturation,45 plays a fundamental role in the crystallization behavior. The threshold supersaturation is a kinetic parameter related to an energy barrier in the same way as the critical supersaturation in CNT. The probability of incorporation of a 2792



■

Article

(29) Wang, Y.; Merino, E. Geochim. Cosmochim. Acta 1992, 56, 587− 596. (30) Lee, W. T.; Salje, E. K. H. Eur. Phys. J. B 2000, 13, 395−398. (31) Ortoleva, P. Geochemical Self-Organization; Oxford University Press: New York, 1994. (32) Mues, T.; Heuer, A.; Burger, M.; Lubashevsky, I. Phys. Rev. E 2010, 81, 051605. (33) Lubashevsky, I.; Mues, T.; Heuer, A. Phys. Rev. E 2008, 78, 041606. (34) Karapiperis, T. J. Stat. Phys. 1995, 81, 165−180. (35) Prieto, M. Rev. Mineral. Geochem. 2009, 70, 47−85. (36) Chang, L. L. Y.; Brice, W. R. Am. Mineral. 1971, 56, 338−341. (37) Kö ningsberg, E.; Hausner, R.; Gamsjäger, H. Geochim. Cosmochim. Acta 1991, 55, 3505−3514. (38) Böttcher, M. E.; Gehlken, P. L. Appl. Spectrosc. 1997, 51, 130− 131. (39) Wang, Q.; de Leeuw, N. H. Mineral. Mag. 2008, 72, 525−529. (40) Curti, E. Coprecipitation of Radionuclides: Basic Concepts, Literature Review and First Applications; PSI Report 97-10; Paul Scherrer Institut: Villigen, 1997. (41) Tesoriero, A. J.; Pankow, J. F. Geochim. Cosmochim. Acta 1996, 60, 1053−1063. (42) Katsikopoulos, D.; Fernández-González, A.; Prieto, M. Mineral. Mag. 2008, 72, 433−436. (43) Sangwal, K. J. Cryst. Growth 1989, 97, 393−405. (44) Pina, C. M.; Enders, M.; Putnis, A. Chem. Geol. 2000, 168, 195− 210. (45) Putnis, A.; Prieto, M.; Fernández-Díaz, L. Geol. Mag. 1995, 132, 1−13. (46) Astilleros, J. M.; Pina, C. M.; Fernández-Díaz, L.; Putnis, A. Surf. Sci. 2003, 545, L767−L773. (47) Shtukenberg, A. G.; Astilleros, J. M.; Putnis, A. Surf. Sci. 2005, 590, 212−223. (48) Pérez-Garrido, C.; Fernández-Díaz, L.; Pina, C. M.; Prieto, M. Surf. Sci. 2007, 601, 5499−5509. (49) Weimar, J. R.; Boon, J. P. Phys. Rev. E 1994, 49, 1749−1752. (50) Kier, L. B.; Cheng, C. K. J. Chem. Inf. Comput. Sci. 1994, 34, 647−652. (51) Kier, L. B.; Seybold, P. G.; Cheng, C. K. Modeling Chemical Systems Using Cellular Automata; Springer: Amsterdam, 2005. (52) Kier, L. B.; Cheng, C. K. J. Math. Chem. 1995, 21, 71−77. (53) Kier, L. B.; Cheng, C. K.; Testa, B.; Carrupt, P. A. J. Pharm. Sci. 1997, 86, 774−779. (54) Henisch, H. K.; García-Ruiz, J. M. J. Cryst. Growth 1986, 75, 203−211. (55) Prieto, M.; Fernández-Díaz, L.; López-Andrés, S. J. Cryst. Growth 1991, 108, 770−778. (56) Speer, A. Rev. Mineral. 1983, 11, 145−190. (57) Chang, L. L. Y. J. Geol. 1965, 73, 346−368. (58) Katsikopoulos, D. Crystallization of Metal-Bearing Carbonate Solid Solutions with the Structure of Calcite at Ambient Conditions: The Cases of Cd2+, Mn2+, and Co2+. PhD Thesis, University of Oviedo, Oviedo, 2009. (59) Sangwal, K. Additives and Crystallization Processes: From Fundamentals to Applications; Wiley: New York, 2007. (60) Shtukenberg, A. G.; Punin, Y. O.; Azimov, P. Am. J. Sci. 2006, 306, 553−574.

AUTHOR INFORMATION

Corresponding Author

*Address: Departamento de Geologia, Faculdade de Ciências ́ C6, Piso 4, Campo Grande, da Universidade de Lisboa, Edificio 1749-016 Lisboa, Portugal. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was partially supported by the Spanish Ministry of Innovation and Science of Spain (Grant: CGL2010-20134). We thank Angeles Fernández-González for kindly providing the BSE images. The authors wish to acknowledge the two anonymous reviewers and their insightful comments that helped us shape a better and clearer manuscript.

■

REFERENCES

(1) Böttcher, M. E.; Dietzel, M. In Ion Partitioning in AmbientTemperature Aqueous Systems; Prieto, M., Stoll, H., Eds.; European Mineralogical Union and Mineralogical Society of Great Britain and Ireland: Twickenham, 2010; EMU Notes in Mineralogy, Vol. 10, Chapter 4, pp 139−187. (2) Prieto, M.; Astilleros, J. M.; Fernández-Díaz, L. Elements 2013, 9, 195−201. (3) Glynn, P. D.; Reardon, E. J.; Plummer, L. N.; Busenberg, E. Geochim. Cosmochim. Acta 1990, 54, 267−282. (4) Prieto, M.; Putnis, A.; Fernández-Diaz, L. Geol. Mag. 1993, 130, 289−299. (5) Prieto, M.; Fernández-González, A.; Putnis, A.; Fernández-Diaz, L. Geochim. Cosmochim. Acta 1997, 61, 3383−3397. (6) Pina, C. M.; Putnis, A. Geochim. Cosmochim. Acta 2002, 66, 185− 192. (7) Watson, E. B. Geochim. Cosmochim. Acta 2004, 68, 1473−1488. (8) Shore, M.; Fowler, A. D. Can. Mineral. 1996, 34, 1111−1126. (9) May, R. M. Nature 1976, 261, 459−467. (10) Wolfram, S. Rev. Mod. Phys. 1983, 55, 601−644. (11) Reeder, R. J.; Fagioli, R. O.; Meyers, W. J. Earth-Sci. Rev. 1990, 29, 39−46. (12) Putnis, A.; Fernández-Diaz, L.; Prieto, M. Nature 1992, 358, 743−745. (13) Fernández-González, A.; Prieto, M.; Putnis, A.; López-Andrés, S. Mineral. Mag. 1999, 63, 331−343. (14) Shtukenberg, A. G.; Punin, Y. O. Mineral. Mag. 2011, 75, 169− 183. (15) Yardley, B. W. D.; Rochelle, C. A.; Barnicoat, C. A.; Lloyd, G. E. Mineral. Mag. 1991, 55, 357−465. (16) Jamtveit, B.; Wogelius, R.; Fraser, D. G. Geology 1993, 21, 113− 116. (17) Jamtveit, B.; Hervig, R. L. Science 1994, 263, 505−508. (18) Wogelius, R. A.; Fraser, D. G.; Wall, G. R. T.; Grime, G. W. Geochim. Cosmochim. Acta 1997, 61, 2037−2051. (19) Jamtveit, B. In Growth, Dissolution and Pattern Formation in Geosystems; Jamtveit, B., Meakin, P., Eds.; Kluwer Academic Publishers: Dordrecht, 1999; Chapter 3, pp 65−84. (20) L’Heureux, I.; Fowler, A. D. Am. Mineral. 1994, 79, 885−891. (21) Ortoleva, P. Earth-Sci. Rev. 1990, 29, 3−8. (22) L’Heureux, I.; Jamtveit, B. Geochim. Cosmochim. Acta 2002, 66, 417−422. (23) Katsev, S.; L’Heureux, I. Phys. Rev. E 2002, 66, 066206. (24) L’Heureux, I.; Katsev, S. Chem. Geol. 2006, 225, 230−243. (25) Kalischewski, F.; Lubashevsky, I.; Heuer, A. Phys. Rev. E 2007, 75, 021601. (26) Holten, T.; Jamtveit, B.; Meakin, P.; Cortini, M.; Blundy, J.; Austrheim, H. Am. Mineral. 1997, 82, 596−606. (27) Holten, T.; Jamtveit, B.; Meakin, P. Geochim. Cosmochim. Acta 2000, 64, 1893−1904. (28) Katsev, S.; L’Heureux, I. Phys. Rev. E 2000, 61, 4972−4979. 2793


Development of Compositional Patterns during the Growth of Solid

Recommend Documents