Microscopic Mechanism and Kinetics of Ice Formation at Complex Interfaces: Zooming in on Kaolinite.

Most ice in nature forms because of impurities which boost the exceedingly low nucleation rate of pure supercooled water. However, the microscopic details of ice nucleation on these substances remain largely unknown. Here, we have unraveled the molecular mechanism and the kinetics of ice formation on kaolinite, a clay mineral playing a key role in climate science. We find that the formation of ice at strong supercooling in the presence of this clay is about 20 orders of magnitude faster than homogeneous freezing. The critical nucleus is substantially smaller than that found for homogeneous nucleation and, in contrast to the predictions of classical nucleation theory (CNT), it has a strong two-dimensional character. Nonetheless, we show that CNT describes correctly the formation of ice at this complex interface. Kaolinite also promotes the exclusive nucleation of hexagonal ice, as opposed to homogeneous freezing where a mixture of cubic and hexagonal polytypes is observed.

The formation of ice is at the heart of intracellular freezing,[1] stratospheric ozone chemistry,[2] cloud dynamics,[3] rock weathering,[4] and hydrate formation.[5] Because ice nucleation within pure supercooled liquid water is amazingly rare in nature, most of the ice on earth forms heterogeneously, in the presence of foreign particles which boost the ice nucleation rate.[6] These substances, which can be as diverse as soot,[7] bacterial fragments,[8] or mineral dust,[9] lower the free-energy barrier for nucleation and make ice formation possible even at a few degrees of supercooling. However, the microscopic details of heterogeneous ice nucleation are still poorly understood. State-of-the-art experimental techniques can establish whether a certain material is efficient in promoting heterogeneous ice nucleation, but it is very challenging to achieve the temporal and spatial resolutions required to characterize the process at the molecular level. On the other hand, spontaneous fluctuations that produce nuclei of critical size are rare events. They thus happen on time scales (e.g. seconds) that are far beyond the reach of classical molecular dynamics simulations. This is why, to our knowledge, quantitative simulations of heterogeneous ice nucleation have been successful only when using the coarse-grained mW model for water.[7,10,11] Such simulations have gone a long way toward improving our fundamental understanding of heterogeneous ice nucleation, but coarse-grained models are not appropriate for many of the more complex and interesting ice-nucleating substrates.

A representative example is the formation of ice on clay minerals, a phenomenon critical to cloud formation and dynamics.[9,12] For instance, the heterogeneous ice nucleation probability in the presence of kaolinite, a clay mineral well studied by both experiments[6,13−16] and simulations,[17−20] seems to be related to its surface area,[13] but how exactly this material facilitates the formation of the ice nuclei is largely unestablished. Classical molecular dynamics simulations have recently succeeded in simulating ice nucleation on kaolinite.[19,20] However, finite size effects[19] and rigid models of the surface[20] prevented the extraction of quantitative results. In fact, it is exceedingly challenging to compute via atomistic simulations ice nucleation rates, which have been inferred (for homogeneous freezing only) along a wide range of temperatures[21] and recently computed directly at strong supercooling (ΔT = 42 K) for the fully atomistic TIP4P/Ice model of water.[22]

In this work, we have computed the rate and unraveled the mechanism at the all-atom level of the heterogeneous nucleation of ice on a substrate of practical relevance. We have considered the hydroxylated (001) surface of kaolinite as a prototypical material capable of promoting ice formation. We quantify the efficiency of kaolinite in boosting ice nucleation and find that this mineral alters the ice polytype that would form homogeneously at the same conditions. We also observe that ice nuclei grow in an anisotropic fashion, in contrast with the predictions of classical nucleation theory (CNT) which nonetheless we demonstrate is reliable in describing quantitatively the heterogeneous nucleation process.

Kaolinite (Al2Si2O5(OH)4) is a layered aluminosilicate, in which each layer contains a tetrahedral silica sheet and an octahedral alumina sheet, in turn terminated with hydroxyl groups. Facile cleavage along the (100) basal plane parallel to the layers results in surfaces exposing either the silica-terminated face or the hydroxyl-terminated one. The latter is believed to be the most effective in promoting ice nucleation, as the hydroxyl groups form a hexagonal arrangement that possibly templates ice formation.[19,23] Here we considered a single slab of kaolinite cleaved along the (100) plane so that it exposes the hydroxyl-terminated surface, while water molecules have been represented by the fully atomistic TIP4P/Ice model.[24] Further details about the structure of the water–kaolinite interface and the computational setup can be found in refs (19 and 25) and in the Supporting Information.

The heterogeneous ice nucleation rate was obtained using the forward flux sampling (FFS) technique,[26] which has been successfully applied to homogeneous water freezing[27,28] and to diverse nucleation scenarios.[29−32] Within this approach, the path from liquid water to crystalline ice is described by an order parameter, λ. A set (Nλ) of discrete interfaces characterized by an increasing value of λ is identified along this order parameter. Here, we have chosen λ as the number of water molecules in the largest icelike cluster within the whole system plus its first coordination shell (see the Supporting Information). The natural fluctuations of the system at each interface, sampled by a collection of unbiased molecular dynamics simulations, are then exploited and the nucleation rate, J, is calculated usingwhere Φλ0 is the rate at which the system reaches the first interface λ0. The total probability P(λ|λ0) for a trajectory starting from λ0 to reach the ice basin is decomposed into the product of the crossing probabilities P(λ|λ) . The details of the algorithm are described in the Supporting Information.

To compare our results with the homogeneous data from ref (22), we have performed FFS simulations at the same temperature, T = 230 K, corresponding for the TIP4P/Ice model to ΔT = 42 K. The calculated growth probability, P(λ|λ0), as a function of λ is reported in Figure a. In contrast with the transition probability for homogeneous nucleation reported in ref (22), we do not observe any inflection region, i.e., a regime for which P (λ|λ0) decreases sharply (P(λ|λ) < P(λ|λ) for some i > j). This inflection is because in the early stages of homogeneous nucleation the largest nuclei are mostly made of hexagonal ice (Ih), which leads to rather aspherical nuclei that are very unlikely to survive and reach the later parts of the nucleation pathway. Within the inflection region the nuclei contain instead a substantial fraction of cubic ice (Ic). It seems that in forming this polytype the nuclei are able to adopt a more spherical shape and that this is essential for ultimately growing toward the critical nucleus size. In contrast, within this heterogeneous case, the presence of the surface allows this process of forming spherical Ic-rich crystallites to be bypassed. Here, nucleation proceeds exclusively heterogeneously at the kaolinite–water interface. During the early stages of the process, the fraction of ice nuclei on the surface (as defined in the Supporting Information) is only around 25%, as shown in Figure a, because at this strong ΔT natural fluctuations toward the ice phase are abundant and homogeneously distributed throughout the liquid. However, as nucleation proceeds, the nuclei within the bulk of the liquid slab become less favorable, until only nuclei at the water–kaolinite interface survive. From this evidence alone one can conclude that at this temperature kaolinite substantially promotes the formation of ice via heterogeneous nucleation.

Figure 1

(a) Calculated growth probability P(λ|λ0) and fraction of ice nuclei sitting on top of the kaolinite (001) hydroxylated surface as a function of λ. (b) Committor probability PC(λ) as a function of λ. The value of PC(λ) = 0.5 corresponds to the critical nucleus size NC = 225. A typical ice nucleus of critical size is shown in the insets.

Our FFS simulation results in a heterogeneous ice nucleation rate of JHetero = 1026±2 s–1 m–3, which can be compared with the homogeneous nucleation rate of JHomo = 105.9299±0.6538 s–1 m–3 reported in ref (22). The hydroxylated (001) surface of kaolinite thus enhances the homogeneous ice nucleation rate by about 20 orders of magnitude at ΔT = 42 K. This spectacular boost is similar to that reported for simulations of heterogeneous ice nucleation on graphitic surfaces[11] and on Lennard-Jones crystals[10] at similar ΔT using the coarse-grained mW model.

An estimate of the critical nucleus size, NC, can be obtained directly from the crossing probabilities assuming that λ is a good reaction coordinate for the nucleation process.[22] In this scenario, NC is the value for which the committor probability PC(λ) for the nuclei to proceed toward the ice phase instead of shrinking into the liquid is equal to 0.5. As shown in Figure b, PC(λ) = 0.5 corresponds in our case to a critical nucleus of 225 ± 25 water molecules. The estimate of the homogeneous critical nucleus size, obtained by means of the same approximate approach employed here, is NC = 540 ± 30 water molecules (as obtained by using the definition of λ employed in this work, see the Supporting Information), more than two times larger than our estimate for the heterogeneous case.

At this supercooling, homogeneous water nucleates into stacking disordered ice (a mixture of Ih and Ic).[33−35] However, the presence of the clay leads to a very different outcome. To analyze the competition between Ih and Ic, we have adopted the topological criterion introduced in ref (22) (see the Supporting Information), pinpointing the building blocks of Ic (double-diamond cages, DDC) and Ih (hexagonal cages, HC) within the largest ice nuclei. The results are summarized in Figure : for ice nuclei in the bulk, a slightly larger fraction of HC with respect to DDC develops until they disappear because of the dominance of the much more favorable nuclei at the surface. In contrast, nuclei at the surface contain a large fraction of HC from the earliest stages of the nucleation, and they exclusively expose the prism face of Ih to the hexagonal arrangement of hydroxyl groups of the clay. This is consistent with what has been suggested previously by classical MD simulations[19,20] and demonstrates that at this supercooling heterogeneous nucleation takes place solely via the hexagonal ice polytype, in contrast with homogeneous nucleation. Experimental evidence[35] suggests that stacking disordered ice on kaolinite is likely to appear after the nucleation process because of the kinetics of crystal growth and the presence of surfaces other than the hydroxylated (001).

Figure 2

Average number of double-diamond cages ⟨DDC⟩Bulk and hexagonal cages ⟨HC⟩Bulk within the largest ice nuclei (identified according to the order parameter λ) in the bulk of the liquid slab only as a function of λ (nuclei in the bulk disappear beyond the value of λ marked by the vertical green line). Averages for the largest ice nuclei sitting on top of the kaolinite (001) hydroxylated surface (⟨DDC⟩Surf and ⟨HC⟩Surf) are also reported. The insets depict DDC and HC within an ice nucleus in the bulk at the early stages of nucleation (left) and a postcritical ice nucleus at the water–clay surface (right). Oxygen atoms belonging to the largest ice nucleus (hydrogens not shown) are depicted in blue (DDC), red (HC), and yellow (both DDC and HC). Atoms belonging to the largest ice nucleus but not involved in any DDC or HC are shown in gray.

In the homogeneous case, critical nuclei tend to be rather spherical even at this strong supercooling.[22] However, we see a very different behavior here. This is illustrated in Figure , where we show as a function of λ the asphericity parameter α (which is equal to zero for spherical objects and one for an infinitely elongated rod), for nuclei in the bulk and at the surface. Note that heterogeneous CNT predicts (on flat surfaces) critical nuclei in the form of spherical caps, the exact shape of which is dictated by the contact angle, θIce,Surf , between the nuclei and the surface.[11] For instance, α = 0.094 for a pristine hemispherical cap, corresponding to θIce,Surf = 90°. Also reported in Figure is the spatial extent Δz of the nuclei along the direction normal to the slab (the exact definitions of α and Δz are provided in the Supporting Information). Nuclei within the bulk tend to be rather spherical. A small increase in the asphericity is observed right before these nuclei disappear and are replaced with nuclei at the surface. This regime, in which the nuclei in the bulk grow substantially and become less spherical, possibly corresponds to the onset of the inflection region observed within the homogeneous case. However, here nucleation is dominated by the surface. While nuclei at the surface are initially quite similar to spherical caps, they tend to grow by expanding at the water–kaolinite interface because of the favorable templating effect of the hydroxyl groups, which favors the formation of the prism face of Ih.[19] This can clearly be seen by looking at the substantial increase in α for the nuclei at the surface, which is accompanied by a slight drop in Δz corresponding to an expansion of the nuclei in two dimensions. Once the nuclei have overcome the critical nucleus size, they tend to return to a more isotropic and compact form, while accumulating new ice layers along the normal to the surface. We note that due to the strong two-dimensional nature of the critical ice nuclei, special care has to be taken to avoid finite size effects. We have therefore used a simulation box with lateral dimensions of the order of 60 Å, which is large enough to prevent interactions between the ice nuclei and their periodic images (as discussed in the Supporting Information).

Figure 3

Asphericity parameter, α, and spatial extent of the ice nuclei along the direction normal to the clay slab, Δz, as a function of λ for ice nuclei in the bulk (αBulk and ΔzBulk). Nuclei in the bulk disappear beyond the value of λ marked by the vertical green line. Averages within the ice nuclei sitting on top of the kaolinite (001) hydroxylated surface (αSurf and ΔzSurf) are also reported. The insets correspond to typical ice nuclei containing about 105, 200, and 325 (from left to right) water molecules.

The fact that the system reaches the critical nucleus size by expanding chiefly in two dimensions is in sharp contrast with the heterogeneous nucleation picture predicated by CNT. Hence, the question arises: Is CNT able to describe heterogeneous ice nucleation on a complex substrate at this strong supercooling? Strikingly, the answer is yes. To show this, we compare the shape factor for heterogeneous nucleation , customarily used in CNT[36] to quantify the net effect of the surface on the free-energy barrier for nucleation ΔGC, with the volumetric factor . Details of this comparison are included in the Supporting Information. Note that while different, equally valid ways of defining an icelike cluster can lead to different values of NC, there is no ambiguity in the estimate of and as long as the same order parameter is used to define both NHomoC and NHeteroC. Thus, we obtain , in very good agreement with . Heterogeneous CNT has already proven to be reliable in describing the crystallization of ice on graphitic surfaces,[11] a scenario very different from ice formation on kaolinite. In fact, while the size of the critical clusters reported in ref (11) is similar to what has been obtained here (a few hundred water molecules), critical ice nuclei of mW water on graphitic surfaces are shaped as spherical caps, in line with CNT assumptions. This is due to the fairly weak interaction between water and carbonaceous surfaces,[7] which results in a weak wetting of the ice phase on the substrate. In contrast, our results show that ice nuclei on kaolinite tend to wet substantially the substrate, leading to shapes very different from spherical caps. In this regime, where the nuclei are small and the ice–kaolinite contact angle, θIce,Surf , is also small, line tension at the water–ice–kaolinite interface could introduce a mismatch between and (see e.g. refs (37 and 38)). However, this is not the case, as CNT holds quantitatively for the formation of ice on kaolinite even at the strong supercooling probed in this work.

The value of JHomo reported in ref (22) is about 11 orders of magnitude smaller than the experimental value extrapolated from ref (39). In addition, at the strong supercooling of ΔT = 42 K, no direct measures of JHetero exist for kaolinite (nor indeed for the homogeneous case), as pure water freezes homogeneously at T < ΔT ∼ 38 K. Consequently, extrapolations are necessary, leading to experimental uncertainties as large as 6 orders of magnitude.[40] Nonetheless, our quantifies the relative ice nucleation ability of kaolinite with respect to the homogeneous case, which can thus be compared with experimental values. Estimates of from measurements of ice formation on kaolinite particles can vary from 0.23 to 0.69 according to the interpretation of the experimental data,[41] and the seminal work of Murray[13] suggests a value of 0.11 for the exclusive formation of Ih observed in this work. The variability of these experimental results stems mainly from the diversity of the kaolinite samples (in terms of, for example, shape, purity, and surfaces exposed, the latter still largely unknown) and the difficulty interpreting the experimental data using heterogeneous CNT, for which tiny changes in quantities such as the free-energy difference between water and ice lead to substantial discrepancies.[22] To date, experiments must address populations of uneven particles and different nucleation sites. Here we provide a value of for a perfectly flat, defect-free (001) hydroxylated surface of kaolinite, with the hope of aiding the experimental investigation of well-defined, clean kaolinite substrates in the near future. We also note that our simulations of crystal nucleation are the very edge of what molecular dynamics can presently achieve. However, there is still room for improvement. For instance, heterogeneous ice formation can be affected by the presence of electric fields,[42,43] and similarly, water dissociation is common on many reactive surfaces;[44] these effects cannot be accounted for at present with the traditional force fields employed here.

In summary, we have calculated the heterogeneous ice nucleation rate for a fully atomistic water model on a prototypical clay mineral of great importance to environmental science. We have demonstrated that the hydroxylated (001) surface of kaolinite boosts ice formation by 20 orders of magnitude with respect to homogeneous nucleation at the same supercooling. We have found that this particular kaolinite surface promotes the nucleation of the hexagonal ice polytype, which forms because of the interaction of the prism face with the templating arrangements of hydroxyl groups at the clay interface. We have also found that ice nuclei tend to expand on the clay surface in two dimensions until they reach the critical nucleus size. This is in contrast with the predictions of CNT, which however holds quantitatively for ice formation on kaolinite even at this strong supercooling. Finally, we provide a value of the heterogeneous shape factor for the defect-free surface considered here, in the first attempt to bring simulations of heterogeneous ice nucleation a step closer to experiments. It remains to be investigated to what extent different surface morphologies can in general affect nucleation rates or alter the ice polytypes which form.