Abstract

Supersaturation-Nucleation-Time (S-N-T) diagrams are shown to be a useful tool to predict nucleation during reactive-transport processes in porous media. Such diagrams can be determined experimentally or estimated from theoretical calculations based on classical nucleation theory. With this aim, a ‘pragmatic’ understanding of the nucleation rate equation is adopted here and the meaning and magnitude of the interfacial tension and induction time discussed. Theoretical diagrams and experimental data are shown to match fairly well as long as there is an appropriate choice of the ‘relevant’ volume for induction-time calculations.

Introduction

Nucleation is in the scientific spotlight, as shown by the number of reviews published on the subject in the last few years, particularly in relation to crystal nucleation in solution (Benning and Waychunas, 2008; Vekilov, 2010; Yi and Rutledge, 2012; Gebauer et al., 2014). These additions to the literature are no doubt related to the detection of stable solute species, referred to as pre-nucleation clusters (PNCs), which are supposed to mediate the development of the crystal phase from the aqueous solution (Gebauer et al., 2008; Gebauer and Cölfen, 2011). That breakthrough challenges our idea of the mechanisms involved in the nucleation process as defined in classical nucleation theory (CNT). The use of state-of-the-art techniques (Pouget et al., 2009) and molecular simulation methods (Demichelis et al., 2011) reveals the nature of PNCs and making the statement “now you see them” (Meldrum and Sear, 2008) convincing, at least in the case of calcium carbonate. Other inorganic systems such as certain iron oxy-hydroxides (e.g. Yuwono et al., 2010), calcium phosphates (e.g. Wang et al., 2012), and silica could be reasonably accommodated in the new nucleation paradigm, although the recognition as PNCs of certain oligomers, amorphous nano-clusters, and other “primary” particles is a source of controversy (Gebauer et al., 2014). Moreover, the stability of many of these species and their effect on nucleation barriers remains a point for discussion (De Yoreo, 2013). Finally, even if ‘primary’ particles exist in a given aqueous system, the formation of precursor phases by aggregation may not occur (Baumgartner et al., 2013), and the process can still be described within the CNT framework.

Much work needs to be done to extend the aggregation-based mechanism to inorganic systems other than CaCO3, and therefore nucleation will continue to be a leading topic in the years to come. In the meantime, kinetic modelling of nucleation is typically performed using the classical nucleation rate equation, with some exceptions in systems that are known to involve amorphous precursor phases (e.g. Rodríguez-Blanco et al., 2011). CNT is used not only in industrial crystallization protocols (Sangwal, 2007) but also in geochemical modelling (Fritz and Noguera, 2009) and in fundamental studies of precipitation and co-precipitation on mineral surfaces (Fernández-Martínez et al., 2013; Shtukenberg et al., 2005).

In CNT, the nucleation rate (the number of nuclei, N, formed per unit volume and time: J = N/Vt) depends exponentially on the free energy change (ΔGc) involved in the formation of a nucleus of critical size, i.e.  
J=Γexp(ΔGckT)
(1)
where k is the Boltzmann constant, T the absolute temperature, and Γ a pre-exponential factor. In turn, ΔGc depends significantly on the supersaturation (S), according to  
ΔGc=βω2σ3(kTlnS)2
(2)
where ω is the molecular volume in the solid phase, β is a factor that depends on the shape of the nucleus, and σ is the interfacial tension. Most readers are likely to be aware of the weak points implicit in CNT, such as the idea of applying macroscopic concepts (surface and volume) to molecular-scale objects or the assumption that the interfacial tension does not depend on the size of the nuclei. While there is a general agreement on these weak points, the method by which the controlling parameters in equations 1 and 2 are estimated is widely disregarded. For example, the interfacial tension is recognized to represent all the excess energy accumulated in the vicinity of the solution–cluster interface, but in practice, σ is an ‘artificial’ parameter that cannot be determined by direct and unambiguous experiments (Sönhel, 1982). In fact, the opposite approach must be adopted: the only interfacial tension values that are relevant in modelling nucleation kinetics are those obtained in nucleation experiments, i.e. from nucleation-rate or induction-time measurements. Although the term ‘interfacial tension’ is used interchangeably in the literature, the nucleation-derived σ has a different meaning and value from the “surface energies” obtained by contact-angle measurements (Wu and Nancollas, 1999) or water-adsorption calorimetry (Forbes et al., 2011). Moreover, the CNT-derived σ values depend on the aqueous speciation model, the expression chosen for the supersaturation, and the shape factor chosen for the nuclei (Prieto et al., 2012). This fact was well known by the pioneers in this field (Söhnel, 1982) but is frequently disregarded in the literature. Therefore, to model the nucleation kinetics in a given system, all the input data must be consistent with the value of σ chosen for the calculation. Thus, the CNT expression in equation 1 can be imagined as a species of fitting function in which the main fitting parameter is σ. The nucleation barrier depends heavily (ΔGc ∝ σ3) on this empirical magnitude, whereas other less influential parameters (particularly those included in the pre-exponential factor) can be tuned to optimize the results. From that point of view, ΔGc plays a similar role to the energy barriers to form activated complexes, which are typical in chemical kinetics. Moreover, phenomena such as the singular nucleation behaviour of sparingly soluble compounds in ionic solutions (Kowacz et al., 2010) could be modelled by CNT, provided that ad hoc nucleation-derived σ data were available.

According to this pragmatic understanding of CNT, the present study models the nucleation behaviour observed in previous experiments (Putnis et al., 1995), taking baryte as an example system. Such experiments were conducted in U-shaped tubes in which the reacting ions counter-diffuse through a column of porous silica hydrogel. Because the gel properties inhibit advection and convection, the crystallization medium is purely diffusive. Crystals tend to nucleate approximately midway through the diffusion column, which is the zone where the supersaturation builds up most rapidly because the counter-diffusing ions ‘meet’ there. Hydrogel media have been used widely to simulate precipitation and reactive transport in natural environments (Fernández-Díaz et al., 1996; Prieto et al., 2002) and are currently used in bio-mineralization studies (Sancho-Tomás et al., 2013; Nindiyasari et al., 2014). In nature, precipitation frequently occurs in pores of rocks, soils and sediments. Moreover, crystallization of salts in porous building materials (masonry, cement, mortar) has been recognized as crucial for their weathering and decay (Rodríguez-Navarro et al., 2002). Developing predictive models for the nucleation behaviour in these types of systems is worthwhile, therefore. With this aim, we use a type of supersaturation-nucleation-time diagram (S-N-T) that allows the interpretation of nucleation pathways in solution systems. S-N-T diagrams are reminiscent of the temperature-transformation-time (T-T-T) diagrams used to interpret sub-solidus transformations in mineralogy and metallurgy (Putnis, 1992) and can be used similarly. The choice of baryte as an example system is not arbitrary. Baryte-scale formation is a problem in many industrial processes where knowledge of the nucleation kinetics is essential for developing anti-scale strategies. However, while most precipitation studies have been performed using vigorously stirred solutions, in many scenarios the nucleation of baryte occurs in pores or small cavities where the confined solution can become more highly supersaturated than an analogous free solution (Putnis and Mauthe, 2001; Rodríguez-Ruiz et al., 2014). The present study confirms that hydrogels are excellent media to emulate nucleation of minerals in such conditions of confinement.

Methods

Precipitation experiments

The experiments addressed here were described broadly in the original papers by Putnis et al. (1995) and Prieto et al. (1990). The U-tube arrangement in which two solution reservoirs are separated by a diffusion column of silica hydrogel is shown in Fig. 1. The gel contains ~95.6 vol.% solution within interconnecting pores of diameters 0.1–0.5 μm. However, the pore size is not uniform (Henisch, 1988) and secondary pores with diameters >10 μm are common. The initial gel pH is 5.5, and the working temperature is 25°C. In such a device, the experimental inputs are the U-tube dimensions and the initial concentration of the reactants, while the primary outputs are the location of the crystallization zone (cz) and the time (tw) at which the first crystallites are observed by optical microscopy at ×500 magnification. These data, as provided by Putnis et al. (1995) and related papers, are shown in Table 1. The secondary parameters have been recalculated here to make them consistent with the overall calculation procedure. Solution speciation, mass-transfer, and supersaturation have been computed using the geochemical code PHREEQC (Parkhurst and Appelo, 2013) and its default database. PHREEQC allows modelling of one-dimensional diffusion with the same algorithm that was used by Henisch and García-Ruiz (1986) to model the diffusion and precipitation patterns in gel columns. The newly obtained parameters (Ωth and RΩ) differ from the original ones but correlate similarly with each other.

For sparingly soluble minerals, the thermodynamic supersaturation is given by the quotient between the ionic activity product in the current solution and the thermodynamic solubility product. In the case of baryte:  
Ω={Ba2+}{SO42}KBrt
(3)

In Table 1, Ωth represents the value of Ω at tw, and column RΩ lists the supersaturation rate (dΩ/dt) at tw. The waiting time tw can be considered as representative of nucleation: once a crystal nucleates, the growth to reach a visible size (×500) is very fast due to the high supersaturation level in the nucleation area (Table 1). Even so, the imprecision in determining the nucleation event can be expected to depend on the observer's care and thoroughness. Here, we have chosen deliberately an over-estimated value of ±6 h, considering a 3 h monitoring interval. This imprecision is not very significant (<2%) when compared with the tw values but needs to be incorporated in the calculations. The last column in Table 1 (tR) represents the time elapsed to reach Ωth at a constant rate equal to RΩ.

CNT calculations

There are numerous compilations of CNT-derived σ values in the literature, but most of them come from the same original papers by Nielsen (1967), Garten and Head (1973) and a small number of others. These original data were determined using concentrations instead of activities to express the supersaturation or simply introducing correction factors for the degree of dissociation. Here, to use a value consistent with the PHREEQC code, σ has been recalculated from some of these original data. Garten and Head (1973) determined a value of 107 mJ/m2 considering spherical nuclei (β = 16π/3) and a supersaturation scale (s±)) based on mean ionic activities (a±) in the aqueous solution, where a±=(a+ν+aν)1ν, and ν = ν+ + ν is the number of ions in the solute formula. In contrast, Nielsen (1967) determined a value of 135 mJ/m2 using an Ω scale and cubic nuclei (β = 32). We can convert the Garten and Head (G&H) value to Nielsen's scale using the expression:  
σ(Neilsen)=σ(G&H)ν23π(π6)13
(4)

On the basis of Nielsen, the interfacial tension obtained by G&H is converted to 137 mJ/m2 (ν = 2 for baryte), a value very close to Nielsen's result, which confirms the consistency of both data sets. Starting from the experimental data by G&H, a value of 134 mJ/m2 has been recalculated using an Ω-PHREEQC scale and spherical nuclei. Note that, despite the appearance, this value differs significantly from Nielsen's, which was given for cubic nuclei. The spherical approximation is physically reasonable for the size range considered, in which the particles occurring in many processes tend to be spherical (Liu, 1999).

In CNT, the pre-exponential factor is also a key term that represents the rate of attachment of monomers (growth units or molecules) to the critical nuclei. Different approximations to derive this factor can be seen in the literature (Walton, 1969; Lasaga, 1998; Kashchiev, 2000; Sangwal, 2007). When the attachment is controlled by the diffusion of growth units from the solution bulk to the nucleus surface, a comprehensive formulation (Kashchiev, 2000) results in:  
Γ=ZDN1rcAcN0
(5)
where N0 and N1 are concentrations that represent the number of nucleation sites and the number of monomers per unit volume of fluid, respectively. D is the diffusion coefficient of the monomers in the fluid phase, and Ac is the surface area of the critical nucleus, assumed to be a sphere of radius rc (with rc = 2ωσ/kT lnΩ). The term DN1/rc stands for the incoming diffusion flux of monomers to the nucleus surface. Finally, the Zeldovich factor (Z) arises from the steady-state treatment of the problem:  
Z=(ΔGc3πkTnc2)12
(6)
where nc is the number of monomers in the critical nucleus. For homogeneous (HON) nucleation, every molecular position can be considered a potential nucleation site (Kashchiev and Rosmalen, 2003). Therefore, N0 can be equated to 1/vw, where vw is the volume of a water molecule in the solution. The calculations were implemented in Mathcad (MathSoft Inc.) using the parameters compiled in Table 2. The concentration of monomers (N1) depends on the supersaturation and was calculated in each case using PHREEQC. Simulations of 3D-heterogeneous nucleation (HEN) on foreign micro-particles have been performed using a unified operational formula in which HON is treated as a limiting case of HEN with σ replaced by an effective interfacial tension defined as:  
σeff=Φσ
(7)
where the coefficient Φ is a number between 0 and 1 (Sangwal, 2007; Liu, 1999) that depends on both the substrate–nucleus interaction and the geometry of the assemblage.

Calculation results and discussion

Critical supersaturation and threshold supersaturation

The classical nucleation rate equation is based on the assumption that supersaturation is reached instantaneously. For a given supersaturation, the nucleation rate adopts a specific value, and we can represent a function that relates the two parameters. Figure 2 represents such a function for baryte. As shown, the nucleation rate increases so dramatically with supersaturation that the occurrence of a critical value is usually admitted. Above this value, the nucleation becomes catastrophic (J → ∞), and below this value, nucleation decreases quickly to zero. The critical supersaturation is considered to be representative of the metastability limit, which marks the width of the so-called metastable zone. However, the critical supersaturation needs to be defined by choosing an arbitrary value of J (typically 1 s−1 cm−3), and therefore, the metastability limit is ambiguous. Moreover, in most scenarios, supersaturation is not reached instantaneously but increases continuously until nucleation (whether HEN or HON) occurs. That phenomenon also occurs in counter-diffusion experiments, for which the concept of threshold supersaturation was first defined (Putnis et al., 1995; Prieto et al., 1990). The threshold supersaturation, Ωth, is the effective supersaturation at which nucleation occurs and has been shown to depend on the reaction path followed by the system.

Figure 3 shows the increase in supersaturation in cz as a function of time in the case of experiment BRT-1 (see Table 1 also). The highest value is Ωth. The supersaturation increases to reach an almost constant value of RΩ. The inset shows at a larger scale that there is a linear trend for long diffusion times. The supersaturation rate, RΩ, is given by the slope of this line. Figure 3 also shows tR, i.e. the waiting time when supersaturation increases from zero to Ωth at a constant rate equal to RΩ. According to the experiments of Putnis et al. (1995), inspection of Table 1 shows that the greater the supersaturation rate, the greater the threshold supersaturation. The variation of Ωth with RΩ fits an empirical relationship that used to be familiar in the industrial crystallization literature (Putnis et al., 1995). The theoretical understanding of that correlation still requires deeper insight, however.

Induction time and waiting time

The induction time, i.e. the time (ti) that elapses between the ‘instantaneous’ creation of supersaturation and the detection of nucleation, is another important concept in CNT. This time is related to the ability of the solution to remain supersaturated and, with some misgivings, can be considered an “experimental observable” (Kashchiev and van Rosmalen, 2003). In contrast, the waiting time (tw) represents the total time involved in the transport–reaction process, from the beginning of the experiment to the moment at which Ωth is reached and nucleation occurs. Following the comprehensive review by Kashchiev and van Rosmalen (2003), “when the appearance of the very first supernucleus in the solution volume is the event that brings the solution out of metastability”, the induction time is given by:  
ti=1JV
(8)
where V is the solution volume. As Kashchiev and van Rosmalen (2003) noted, equation 8 implies setting N = 1 at t = ti in the expression J = N/Vt. The approach in equation 8 is suitable only if we use experimental techniques that allow counting of N or detection of the appearance of a single supernucleus. Although this point is frequently disregarded in the literature however, equation 8 is no longer valid when the ti values are obtained by techniques that detect a mass of many nuclei.

‘We can see one’

As discussed below, porous media are highly effective at suppressing nucleation. As a consequence, crystallization occurs at high supersaturation, and the nucleation density (N/V) is typically very low. For this reason (and because of the hydrogel transparency), the gel technique provides certain research opportunities that nucleation in free solutions cannot provide. The low nucleation density allows us to detect and count the very first crystallites in the crystallization zone. In fact, certain typical papers on crystallization in gels (Prieto et al., 1990) report not only tw and cz but also the number of crystallites observed. That feature makes equation 8 suitable for exploring the nucleation behaviour in this study.

Figure 4 shows an S-N-T diagram calculated for HON (solid line) using the parameters in Table 2. The calculation has been performed considering the whole volume of the crystallization zone (V = πr2h = 1.91 cm3). The HEN diagram (dashed line) is somewhat speculative because parameters such as the concentration of active nucleation sites (N0) and Φ are difficult to estimate. Here, we have chosen N0 = 2.5 × 107 cm−3 from the data obtained by Garten and Head (1973) for baryte nucleation. Figure 5 shows the crystal counts obtained by these authors as a function of supersaturation. The strong increase above a certain supersaturation level is representative of HON, whereas at a lower supersaturation, the number of nuclei is virtually independent of supersaturation and can be considered representative of the concentration of active HEN sites.

The diagram in Fig. 4 can be used to assess whether nucleation will or will not occur in an evolutionary system and at what degree of supersaturation. However, the diagram does not seem to work with gel experiments. As shown in Fig. 6, the experimental Ωth data points correspond to longer times in comparison with the theoretical nucleation curves. The experimental values of Ωth have been plotted on the ordinate against tR on the abscissa, i.e. considering a stationary supersaturation rate. The open circles show the induction times to be expected at equivalent supersaturation in 1.91 cm3 of ‘free’ solution. As observed, ti and tR differ by orders of magnitude.

The pore size effect

As seen in equation 8, the induction time depends on the solution volume that we consider: the larger V, the shorter ti, in agreement with the probabilistic nature of nucleation. A problem arises because choosing the ‘relevant’ volume for nucleation-rate calculations is not always easy. The diagram in Fig. 4 was calculated by considering the whole volume of the crystallization zone, which does not seem to be the right choice. The fact that porous media are highly effective at suppressing nucleation lies at the core of this problem. During cluster formation, the solution vicinity becomes poor in solute, and the disappearance of a subcritical cluster is ensured if no exchange brings new growth units into its vicinity (Prieto et al., 2002). This effect is particularly important in fine porous media, where the solution is trapped in small pores connected by tortuous routes, and solute transport occurs by diffusion. Under these conditions, each single pore can be envisioned as a crystallization chamber, and we can use the pore size as the ‘relevant’ volume in equation 8.

Figure 7 shows two S-N-T diagrams calculated in this way for pore sizes of 1 μm and 0.1 μm (see Table 2). The experimental results plot between both curves, which demonstrate that the relevant volume in determining the nucleation behaviour is probably the pore volume and that the largest pores must be the preferential nucleation places. Although the HEN curve has not been represented, the high threshold supersaturation points towards an HON mechanism. The only exception occurs in experiment BRT-7. In this case, Ωth and RΩ are comparatively very small and a HEN mechanism could be expected. In counter-diffusion experiments the gel is said to be inert because it does not take part in the precipitation reaction. Nevertheless, in the same way as in free solutions, in real systems nucleation may occur on foreign nanoparticles or active centres present in either the polymeric framework or the interstitial solution. Anyway, Fig. 7 illustrates that S-N-T diagrams can be a suitable tool to assess HEN conditions in rocks, sediments, soils, etc., where the prediction of precipitation reactions is complicated not only by the porous nature of these media but also by their interaction with the precipitate (Stack et al., 2014).

Porous media are known to affect the transport of solutes through a diversity of mechanisms such as adsorption, absorption and precipitation (Brown and Calas, 2013). However, the existence of a pore-size dependence of the metastability of supersaturated solutions is less known and not often considered in reactive-transport modeling. Such dependence usually results in inhibition of precipitation in nanopores and preference for the precipitation in macropores. Stack et al. (2014) studied the kinetics of precipitation of calcium carbonate in an amorphous-silica medium that contained two categories of pore sizes, macropores (Ø >30 μm) and nanopores (Ø = 8–30 nm). Those authors observed precipitation exclusively in the macropores and demonstrated how the interaction between substrate and precipitate is a controlling factor, with more “favorable” interactions allowing precipitation to occur in smaller pores. In fact, precipitation on nanopores can be enhanced using additives that increase the “favourability” of the substrate–precipitate interaction (Stack et al., 2014).

In the CNT framework, such ‘favourability’ can be expressed in terms of the HEN coefficient, Φ. The smaller Φ is, the smaller the effective interfacial tension, σeff (see equation 7) and the larger the HEN probability. The effect is shown in Figure 8a for the nucleation of baryte in 1 μm-sized pores. Figure 8b illustrates the influence of the pore size in the case of HON in a porous medium with two categories (100 μm and 0.1 μm) of pores. As can be seen, for Ω = 6000 the induction time differs by nine orders of magnitude! Obviously, nucleation would occur in the 100 μm pores. Only in case that the supply of reactants continued and Ω increased or kept constant for a long time, precipitation could eventually occur in the nanopores. These two examples show separately the main factors involved in the precipitation behaviour in porous media. In practice, the scenario will be a complex combination of both examples in which the interplay of nucleation mechanisms, pore-size effects, and supersaturation evolution will determine the final outcome.

Conclusions

S-N-T diagrams are shown to be a useful tool for interpreting and predicting nucleation behaviour in porous media. Such diagrams can be determined experimentally or estimated from CNT-based calculations. With this aim, the present author has adopted a pragmatic understanding of CNT in which the nucleation rate expression is imagined as a species of fitting function with the interfacial tension being the main fitting parameter. The method can be extended to a diversity of scenarios. The correct choice of the relevant volume for the studied system becomes clearer when one considers that nucleation in porous rocks is favoured in fractures and open spaces, where a faster supply of growth units is ensured.

Acknowledgements

This work was performed in the framework of the Marie Curie Network ‘Geologic Carbon Storage’ (European Commission, FP7-People-ITN-CO2-REACT-317235). The author acknowledges Dr Cristina Ruiz-Agudo and two anonymous reviewers for their insightful comments. He also thanks the editors for their work and suggestions.

This paper is published as part of a special issue in Mineralogical Magazine, Vol. 78(6), 2014 entitled ‘Mineral–fluid interactions: scaling, surface reactivity and natural systems’.
P1Freely available online through the publisher-supported open access option.