Ice formation on kaolinite: Insights from molecular dynamics simulations

From the immense extent of glaciers to the microscopic length scale of living cells, ice shapes life as we know it.¹ For instance, the formation of clouds² and the weathering of rocks³ originate from water freezing in the atmosphere and on earth, respectively. Moreover, technologies such as cryotherapy and cryopreservation⁴ are greatly influenced by the microscopic details of ice formation. It is surprising to discover that pure water freezes only at very strong supercooling, i.e., when it is brought to temperatures lower than –30 °C.⁵ Water must, therefore, be freezing heterogeneously, as in a world where water would only freeze homogeneously, the Arctic Ocean would hardly turn into the icy playground of polar bears.^6,7 Furthermore, if water only froze homogeneously, too much solar radiation would reach us as there would be no screening by ice-rich clouds.^8,9 Luckily, ice can form at mild supercooling (i.e., at few degrees only below 0 °C) heterogeneously,¹⁰ with the aid of foreign substances which lower the free energy cost needed to form a sufficiently large (or critical) nucleus of crystalline ice within supercooled liquid water. The nature of these impurities is astonishingly diverse:¹⁰ for example, bacterial fragments, soot, pollen, volcanic ashes, and mineral dust have all been shown to boost the rate of ice nucleation. Hence the question: what makes these very different substrates so effective in promoting ice formation? Surprisingly, a conclusive answer is yet to be found, mainly because the microscopic details of heterogeneous ice nucleation are still largely unresolved.¹¹ 

Experiments can quantify the ability of a given substance to promote the formation of ice: for instance, a common approach consists of comparing the fraction of water droplets that freeze in a given time interval at a certain temperature with or without the presence of foreign particles.^12–14 However, the early stages and the atomistic mechanism of nucleation remain exceedingly challenging to probe experimentally. This is in large part because once the critical size (typically of the order of nanometers) of the ice nucleus has been reached, nucleation proceeds on very fast time scales (pico- or nanoseconds). This is why atomistic simulations can provide valuable insight, and complement experimental evidence by unraveling the details of the nucleation mechanism on different substrates at the molecular level. Having said that, however, nucleation is a rare event and seconds can pass before the spontaneous formation of a nucleus of critical size occurs. Running molecular simulations of these lengths is simply not tractable. This is why in the last few years substantial effort has been devoted to developing enhanced sampling techniques capable of tackling this time scale problem,^15,16 such as umbrella sampling,^17,18 forward flux sampling (FFS),^19,20 and transition path sampling.^21,22 Moreover, equally serious issues often go unacknowledged, such as the ability of the interatomic potentials of choice to represent the water, the substrate, and the interaction between the two, or the extent to which simulations of ice nucleation are affected by specific computational details, such as the flexibility of the substrate.

Probing the importance of such aspects would require to investigate heterogeneous ice formation across a collection of different substrates: this is currently possible only by taking advantage of the computational speed granted by the coarse grained mW model²³ for water. This approach led to many important findings, such as the influence of the hydrophobicity and/or lattice mismatch of different crystalline surfaces in promoting ice formation²⁴ and the characterization of ice formation on carbonaceous particles.^25–29 However, fully atomistic models are needed to deal with water at complex interfaces, such as crystalline surfaces of organic crystals or mineral dust particles. In these cases, it remains to be seen whether the description of the surface, and most importantly of the water-surface interaction, is accurate enough to allow for reliable results to be obtained.

Only a few works have probed ice formation on crystalline substrates by means of all-atom models.^30–34 Here, we study heterogeneous ice nucleation at strong supercooling (⁠$ΔT=Tm−T$ = 42 K, where T_m is the melting temperature) on the (001) surface of kaolinite using molecular dynamics (MD). Kaolinite is of great relevance in atmospheric science, and its surface structure and ice nucleation ability have been extensively investigated in both experiments^10,35–38 and simulations.^33,39–41

The building blocks of kaolinite are aluminosilicate layers, i.e., stacking of tetrahedral silicate sheets and octahedral hydroxide sheets.⁴² In the atmosphere, airborne kaolinite particles are found as tiny plate-like particles (with an in-plane size of the order of μm).⁴³ Because of its layered structure, facile cleavage of a kaolinite crystal leads to surfaces exposing the (001) basal planes (hence the plate-like morphology), which in turn can present either the siloxane or the hydroxylated face (see Fig. 1(a)). While the siloxane surface (KAO_Si) is basically hydrophobic and in general thought not to be such an effective ice nucleating agent, first-principles simulations (see, e.g., Ref. 44) have shown that the hydroxylated surface (KAO_OH) is able to stabilize an ice-like bi-layer thanks to its amphoteric nature, characterized by the ability to both accept and donate hydrogen bonds.

In particular, Cox and co-workers³³ performed brute-force MD simulations of ice nucleation on the KAO_OH face of kaolinite using small (10² molecules) models of the clay-water interface. Despite the substantial finite-size effects affecting these simulations, the results suggested that non-basal faces of ice, specifically the primary prism face of the hexagonal ice polytype, can nucleate on the hydroxylated basal face of kaolinite. This evidence has also been observed in brute-force MD simulations⁴¹ employing much larger (10⁴ molecules) models. However, in that case nucleation was observed when the kaolinite surface was almost entirely kept frozen, i.e., atoms have been kept fixed at the experimental atomic positions of the bulk phase. In fact, it has been suggested⁴¹ that the flexibility of the kaolinite surface can substantially affect the time scale over which heterogeneous ice nucleation occurs. Moreover, we have recently succeeded in elucidating the kinetics of ice formation on this mineral³⁴ by means of FFS simulations,¹⁹ an accurate path sampling technique which has been successfully employed to investigate crystal nucleation and growth in different systems.^20,45,46

In this work, we investigate: (i) the type of ice (cubic, I_c, or hexagonal, I_h) that forms on the hydroxylated basal face of kaolinite; (ii) how the surface relaxation of the clay alters ice formation at the water-kaolinite interface; (iii) some aspects of how the force fields used perform for this system.

We demonstrate by seeded MD simulations that the hexagonal polytype of ice is likely to be the only one involved in the nucleation process, due to the favorable interaction between the hydroxylated (001) surface of kaolinite and the prism face of hexagonal ice. In fact, while large (∼450 molecules) I_c seeds dissolve into liquid water, I_h nuclei of the same size or smaller (∼250 molecules) grow within a wide temperature range. In addition, we show by means of very long (up to 2 μs) unbiased MD simulations that nuclei of hexagonal ice exposing the prism face to the clay surface spontaneously occur as natural fluctuations of the water network. These findings are consistent with previous computational studies,^33,34,41 and provide conclusive evidence of the dominant role of the hexagonal polytype in the early stages of ice nucleation on the hydroxylated (001) surface of kaolinite.

We have also addressed whether the flexibility of the substrate, in this case the kaolinite surface, can affect the kinetics of ice formation within MD simulations. We find that equally flexible (or non-flexible) kaolinite surfaces can lead to very different outcomes in terms of ice formation. In particular, it seems that small structural changes can significantly alter the nucleation ability of this mineral. Thus, we assess the sensitivity of the nucleation mechanism to the microscopic structural features of the water-kaolinite interface. We find that surface relaxation, however small, can substantially alter the templating effect of kaolinite and that these effects can facilitate the heterogeneous formation of ice. Specifically, small structural changes in the arrangement of the hydroxyl groups at the surface affect the free energy cost needed to form a hexagonal motif of water molecules within the first overlayer: this templating structure is in turn very effective in promoting the formation of ice.

Finally, we briefly investigate whether the TIP4P/Ice and CLAY_FF force fields are capable to describe supercooled water and kaolinite, respectively. We find that the CLAY_FF force field seems to provide a reliable description of ice nucleation on the (001) hydroxylated surface of kaolinite. However, the surface relaxation of the siloxane (001) surface of kaolinite predicted by the CLAY_FF is not consistent with first principles calculation results, thus putting into question the ability of this force field to deal with ice nucleation and growth on the siloxane face of kaolinite.

MD simulations have been performed using the GROMACS simulation package.⁴⁷ The CLAY_FF⁴⁸ and the TIP4P/Ice⁴⁹ force fields were used to model kaolinite and water, respectively. So as to address the question of surface flexibility and relaxation, we have considered three different water-kaolinite interfaces, where water molecules are in contact with the hydroxylated (001) face of kaolinite (KAO_OH):

• (S1) Frozen surface, unrelaxed: All the kaolinite atoms are kept fixed during the MD simulations at the experimental positions of the bulk system, except for the hydrogen atoms of the hydroxyl groups on the outer layer of the slab. These hydrogen atoms are bonded to the corresponding oxygen atoms via a harmonic constraint characterized by a spring constant of 2.3185$⋅103$ kJ/mol Å⁻² acting on the O–H bond length (1.0 Å), as required by the CLAY_FF force field.⁴⁸ About 6000 water molecules have been placed between two kaolinite slabs mirroring each other. This interface is identical to that reported in Ref. 41.

• (S2) Frozen surface, relaxed: All the kaolinite atoms are kept fixed during the MD simulations at the average atomic positions of the system previously equilibrated in the absence of restraints at 230 K, except for the hydrogen atoms of the hydroxyl groups on the outer layer of the slab. Specifically, the atomic positions have been obtained via an average of 5 ns within the above mentioned equilibrium MD run at 230 K. As such, S2 does not represent a minimum energy configuration of the system. The O–H bonds are treated with the same harmonic constraint of S1. Again in this system about 6000 water molecules have been placed between two kaolinite slabs mirroring each other.

• (S3) Flexible surface: The positions of the silicon atoms within the kaolinite slab are restrained during the MD simulations by means of a harmonic potential characterized by a spring constant of 1.0$⋅103$ kJ/mol Å⁻². The O–H bonds are treated with the same harmonic constraint of S1 and S2. All the other kaolinite atoms are unrestrained. About 6000 water molecules have been placed on top of a single kaolinite slab, as depicted in Fig. 1(a).

Schematic representations of S1, S2, and S3 are shown in Fig. 1(b). Note that upon surface relaxation the arrangement of the hydroxyl groups becomes more symmetric in the xy plane and more corrugated with respect to the z axis (normal to the 001 plane). Additional details about S1, S2, and S3 as well as about additional models for the water-kaolinite interfaces are discussed in the supplementary material. The interaction parameters between the clay and the water were obtained using the standard Lorentz-Berthelot mixing rules,^50,51 which yield water-surface interaction energies in good agreement with high quality reference data⁵² from quantum Monte Carlo calculations.⁵³ The equations of motion were integrated using a leap-frog integrator with a time step of 2 fs. The van der Waals (non bonded) interactions were considered up to 10 Å, where a switching function was used to bring them to zero at 12 Å. Electrostatic interactions have been dealt with by means of an Ewald summation with a real space cutoff at 14 Å. The NVT ensemble was sampled at 230 K using a stochastic velocity rescaling thermostat⁵⁴ with a very weak coupling constant of 4 ps in order to avoid temperature inhomogeneities throughout the system. The geometry of the water molecules (TIP4P/Ice being a rigid model) was constrained using the SETTLE algorithm⁵⁵ while the P-LINCS algorithm⁵⁶ was used to constrain the O–H bonds within the clay. The system was equilibrated at 300 K for 10 ns, before being quenched to 230 K over 50 ns. For each interface S1, S2, and S3 (plus the other interfaces discussed in the supplementary material), 10 independent MD simulations have been performed to look for nucleation events and in order to extract the equilibrium canonical averages used to compute free energies. Details concerning the order parameter used to identify ice-like water molecules are reported in the supplementary material.

We have also performed density functional theory (DFT) calculations as part of this study: periodic kaolinite slab models were used within the plane wave pseudopotential approach, using both GGA (Generalized Gradient Approximation) and dispersion inclusive GGA exchange-correlation functionals. Full details are included in the supplementary material.

At strong supercooling (⁠$ΔT≲$40 K), homogeneous water freezing results in a mixture of I_h and I_c that is known as stacking disordered ice^57–61 (I_sd). However, previous computational studies^33,41 suggest that the formation of ice at the water-kaolinite interface proceeds via the nucleation of I_h only, possibly due to the favorable interaction between its prism face and the hydroxyl groups on the clay surface.

Here, we have explicitly compared the preference of the hydroxylated (001) surface of kaolinite for nucleating either I_h or I_c by means of seeded MD simulations. Specifically, we have inserted several crystalline nuclei of either cubic or hexagonal ice into the system, and then subsequently observed at which temperature they shrink into the liquid phase and at which temperatures they proceed toward full crystallization. Seeded MD simulations are an efficient way of obtaining a qualitative picture of crystal nucleation and growth, having been successfully used recently to explore homogeneous water freezing.^62–64 In the case of heterogeneous ice nucleation, however, one serious issue with seeded MD simulations is the choice/construction of the crystalline seeds. In fact, it is: (i) rather difficult to guess a priori which crystallographic face—if any—of a certain ice polytype will form at the water-kaolinite interface and (ii) it is even more challenging to construct a feasible hydrogen bond network between the ice seed and the surface. In this work it is clear how to resolve problem (i) as we already know that the prism face of I_h and the basal face of I_c bind to kaolinite most strongly.³³ As for the hydrogen bond network, we have employed metadynamics simulations^65,66 to generate reasonable models of I_h and I_c spherical caps in contact with the kaolinite surface, as depicted in Fig. 2 and further detailed in the supplementary material. The ice nuclei obtained in this way (details are discussed in the supplementary material) have a very good match between the crystalline ice seeds and the kaolinite surface, that would be very difficult to obtain otherwise. We have considered ice nuclei containing ∼250 or 450 water molecules at a flexible interface (model S3, as discussed below). A criterion based on the q₃ Steinhardt order parameter⁶⁷ has been used to label each water molecule in the system as liquid, ice-like, and/or belonging to the cubic or hexagonal polytype, as described in the supplementary material. The starting configurations that were obtained from the metadynamics simulations have been equilibrated at 220 K for 200 ps, before collecting a series of unbiased MD runs at different temperatures (220, 225, 230, 235, and 240 K), that were initiated by randomizing the initial atomic velocities in a manner consistent with the corresponding Maxwell-Boltzmann distribution at the temperature of interest.

The results are summarized in Fig. 2, where we report the number of I_c and I_h molecules within the ice seeds as a function of simulation time. We find that seeds containing as many as 450 I_c molecules are not stabilized by the presence of the surface and dissolve at any temperature we have probed (see Figs. 2(a) and 2(b)). In contrast, 450-molecule I_h seeds clearly grow up to 235 K (see Fig. 2(c)), and even small 250-molecule I_h nuclei proceed toward crystallization below 235 K (see Fig. 2(d)). These findings indicate that the basal face of cubic ice is exceedingly unlikely to form on the KAO_OH surface, and that the heterogeneous critical nucleus size $N*$ for ice on top of this kaolinite surface is of the order of 250 molecules at 230 K. This is consistent with the estimate of $N*=$ 225 ± 25 molecules obtained via the FFS simulations reported in Ref. 34. Note that the homogeneous critical nucleus size at the same supercooling is more than two times larger: ∼540 molecules,⁶⁸ where the same order parameter has to be used to compare our results with those of Ref. 68, as discussed in the supplementary material. Thus, our findings suggest that ice nucleation on the KAO_OH surface most likely proceeds via the formation of the prism face of I_h, in agreement with the simulations of Cox et al.³³ and Zielke et al.⁴¹ as well as with our FFS simulations.³⁴ 

Before moving on we stress, however, one drawback of seeded MD simulation: they assume a priori the composition, the structure, the size, and the shape of the crystalline seeds. Thus, even if the latter do grow, one has to verify that such seeds are compatible with the spontaneous fluctuations of the system, in this case of the water network. To explore this, we have performed a very long (2 μs long) unbiased MD simulation, looking at the natural, pre-critical fluctuations of the water network toward the ice phase. In particular, we have determined for each configuration along the trajectory whether the largest ice nucleus in the system is predominantly⁶⁹ made of either I_h or I_c. The probability density of the distance of the center of mass of the ice nuclei from the kaolinite surface (along the z direction normal to the surface) is reported in Fig. 3 for both I_h and I_c nuclei. While a similar fraction of I_h and I_c nuclei can be observed within the bulk of the water slab, close to the water-kaolinite interface there is a very clear preference for I_h nuclei, while the probability for the I_c polytype drops sharply. Importantly, we find that of all the large (i.e., containing more than 60 water molecules) pre-critical ice nuclei, a substantial fraction (25%) still sit on top of the kaolinite surface (that is, at least one water molecule that belongs to the ice nucleus sits within the first overlayer on the clay). In addition, 98% of this subset consists of nuclei of I_h that expose the prism face to the kaolinite surface (as illustrated in the inset of Fig. 3). We also note that TIP4P water models predict only a negligible enthalpic difference of about 0.002 J/mol between I_h and I_c, while experiments and first principles results give an enthalpic difference that ranges between 35 and 160 J/mol.⁷⁰ This difference seems to be due to the absence of anharmonic effects in TIP4P water.⁷¹ Thus, the strong preference of the KAO_OH surface for I_h that we observe is even more significant, as in the homogeneous case the tiny enthalpic difference results in a competition between I_h and I_c that here is outshined by the most favorable interaction between the hydroxyl groups of the KAO_OH surface and the prismatic face of I_h. We can thus safely assert that the I_h nuclei we have investigated with seeded MD can indeed form on the KAO_OH face within spontaneous, pre-critical fluctuations of the water network, and that I_h nuclei exposing the prism face to the kaolinite surface are the most likely to nucleate at this supercooling.

It has been suggested that the ability of the KAO_OH face to promote the formation of ice nuclei stems from the templating effect of the hexagonal arrangement of hydroxyl groups,^33,73 as depicted in Fig. 1. Moreover, recent MD simulations⁴¹ performed at strong supercooling (⁠$ΔT$ = 42 K) that employed the TIP4P/Ice water model indicate that the flexibility of these hydroxyl groups is crucial when it comes to promoting heterogeneous ice nucleation. Specifically, the formation of ice on the KAO_OH face happens spontaneously in unbiased MD simulations within ∼100 ns, provided that all the kaolinite atoms except for the hydrogen atoms of the hydroxyl groups at the water-kaolinite interface are frozen at the experimental positions of the bulk system. We refer to this setup as interface S1 (see Sec. II). On the other hand, restraining the dynamics of the oxygen atoms of the same hydroxyl groups using a harmonic potential, as opposed to completely freezing the atomic positions, prevented the formation of ice within the μs time scale.

In order to investigate how exactly the flexibility of the substrate influences the formation of ice in MD simulations, we have considered in addition to the S1 interface discussed in Ref. 41 two different water-kaolinite interfaces S2 and S3 (see Sec. II). In S2, kaolinite atoms are kept frozen during the MD simulations exactly as for S1, but the starting configuration has been obtained from a kaolinite slab previously equilibrated (without imposing any restraint) at 230 K. In S3, the kaolinite atoms are instead unconstrained and the surface is flexible. Upon running unbiased MD simulations on these three substrates, we found that no nucleation events occur for both S2 and S3 up to 2 μs, whereas on S1 ice formation occurs rapidly (within ∼100 ns). We expected this result for S3, as according to Ref. 41 the flexibility of the surface should prevent ice formation on this time scale. However, the fact that we do not observe ice nucleation on the S2 surface is surprising, as they are both frozen surfaces and the only difference between S1 and S2 is that the latter was relaxed prior to the MD simulation of the frozen substrate.

The structure of the KAO_OH face facilitates the formation of a hexagonal motif of water molecules within the first overlayer (depicted in the inset of Fig. 4(a)). This structure is in turn compatible with the prism face of hexagonal ice or the basal face of cubic ice.^33,41 The likelihood for formation of ice on the KAO_OH face is related to the free energy cost needed to form a hexagonal patch (or connected cluster of water molecule hexagons) containing N_hP water molecules at the water-clay interface. Here we have quantified this free energy cost for the three interfaces S1, S2, and S3 as follows: we start by pinpointing six-membered rings of oxygen atoms within the first water over layer on the surface. To do that, we took advantage of the R.I.N.G.S.⁷⁴ code, and analyzed our simulations with a nearest-neighbors distance set to 3.2 Å and King’s shortest path criterion.^75,76 We then selected only those rings for which each triplet i, j, k of adjacent oxygen atoms is characterized by an in plane angle $αh$ = 120 ° ± 20 °, where

Here, r_x(ij) is the x component of the distance vector between oxygen atoms i and j. We then calculate the number of water molecules N_hP within the largest connected clusters of hexagonal rings. We have collected about 10⁴ values of N_hP from the first 10 ns of 10 independent MD runs for each interface. As S1, S2, and S3 have a slightly different surface area S_A, we have normalized by the latter the value of N_hP. We subsequently evaluate the free energy profile as a function of $NhPSA$ as

where k_B is the Boltzmann constant and $P(NhPSA)$ is the equilibrium⁷⁷ probability density distribution for $NhPSA$⁠. The results, obtained by taking into account the first 10 ns of 10 independent simulations for each interface (in order to avoid the onset of ice formation for S1), are reported in Fig. 4(a). In the case of S1, the free energy profile is rather shallow: for instance, ∼1 k_BT is sufficient to produce a rather large hexagonal patch containing ∼90 water molecules. As a result, the whole first overlayer of water molecules relaxes into this hexagonal pattern within 50–100 ns, quickly triggering crystallization (as shown in Fig. 4(b)). This is why at this interface the formation of ice does not proceed via nucleation events, but instead through a relaxation process. In fact, the onset of crystallization is determined by the time needed for the first water overlayer to relax into the hexagonal template. We have verified that for ten independent simulations the induction times for “nucleation” at the S1 surface are all very similar, thus resulting in a survival probability for the liquid phase (reported in the supplementary material) which is typical of a relaxation process, as opposed to the stochastic nature of nucleation events. We also note that the kinetics of ice formation on S1 is nonphysical, being about six orders of magnitude faster than the nucleation rate we have obtained for S3 via FFS calculations.³⁴ 

On the other hand, the occurrence of large (⁠$NhP≳$50) hexagonal patches for S2 and S3 is exceedingly rare compared to S1: for instance, the same free energy cost needed for S1 to form a hexagonal patch containing ∼90 water molecules results for both S2 and S3 in a patch about two times smaller. Hence, despite the fact that S2 is kept frozen while S3 is fully flexible, the two interfaces show very similar free energy profiles. Indeed the free energy cost needed for the interface to form a templating water overlayer is the same for S2 and S3, and as a consequence, no ice nucleation is observed for either interface on the μs time scale (as illustrated in Fig. 4(b)). We note that we observe a very similar scenario for S1 when we restrain the oxygen atoms of the hydroxyl groups using a harmonic potential instead of freezing them completely: consistent with Ref. 41, we do not observe ice formation. This happens because the harmonic restraints do not prevent surface relaxation from taking place. Consequently, this constrained S1 surface resembles S2. This relaxation, however small, is enough to alter the ability of KAO_OH to produce a large enough hexagonal patch. In fact, we have also verified that by just relaxing the S1 interface at zero K, without equilibrating the kaolinite surface at 230 K, and subsequently freezing the atomic position exactly as we did for S1, we obtain results very similar to what we observe for S2. In particular, ice does not form within the μs time scale, and the extent of surface relaxation and the free energy cost to create the templating hexagonal overlayer are comparable to the outcomes of the S2 scenario.

To understand these results, we have examined the structures of the various slabs with the CLAY_FF force field and also DFT. The changes in the structure of the KAO_OH (001) surface upon relaxation, particularly the arrangement of the hydroxyl groups, are summarized in Fig. 5. We have quantified the corrugation (or surface roughness) of the surface in terms of the deviation of the height $rOHz$ (along the z direction normal to the slab) of the oxygen atoms of the hydroxyl groups from their mean height $⟨rOHz⟩$⁠, as shown in Fig. 5(a). While S1 is basically flat, we observe a small degree of corrugation, up to 0.2 Å, for S2 and S3. Relaxing S1 at the DFT level leads to a very similar degree of corrugation. The in-plane arrangement of the hydroxyl groups at the surface is also affected by surface relaxation. Fig. 5(b) shows the probability density distribution for the deviation of the in-plane (xy plane parallel to the slab) nearest neighbor distance $dOH,ijxy$ of the oxygen atoms of the hydroxyl groups from their mean nearest neighbor distance $⟨dOH,ijxy⟩$ for S1, S2, and S3. Both S2 and S3 are characterized by a non negligible spread of $dOH,ijxy$⁠, which in turn leads to a more symmetric arrangement (see Fig. 5(b)). The same conclusion holds for the structure obtained upon DFT relaxation of S1, albeit the extent of surface relaxation appears to be less pronounced. Despite the fact that the overall extent of surface relaxation for the KAO_OH surface appears to be quite small (ca. 0.1 Å), these marginal structural changes can play a significant role when it comes to the formation of ice on this clay. This is not entirely unexpected, as we have recently shown²⁴ that the ice nucleation rate for the coarse grained mW model of water²³ on Lennard-Jones crystals (at $ΔT$ = 70 K) can change by several orders of magnitude just because of deviations of ca. 0.2 Å in the lattice parameter of the crystalline surface.

It is important to note that S2 differs from S3 in terms of dynamical properties. For instance, the structural relaxation time (defined and discussed in the supplementary material) of the water network within the first overlayer for S2 is two times larger than that obtained for S3. This means that water dynamics at the interface with the KAO_OH face is slower for S2 (and qualitatively for S1 as well, as discussed in the supplementary material), as the frozen surface interacts more strongly with the water molecules, in a manner which is consistent with early findings for Lennard-Jones interfaces.⁷⁸ The absence of nucleation for S2 and S3 on the same μs time scale indicates that the dynamics of water at the kaolinite-water interface has a lesser impact than the structure of the clay on the tendency for ice to form.

This observation that the degree of surface relaxation can strongly affect the nucleation dynamics leads one to ask whether other technicalities play a role. For instance, kaolinite slabs have a non-zero dipole moment, which in principle can affect the nucleation process: in fact, it has been reported that electric fields can affect the freezing of water.^30,31 Moreover, water-surface and/or water-vacuum interfaces introduce structural and dynamical fluctuations into the water network, which must decay within the thickness of the water film on top of the mineral. However, it seems that these issues do not affect the outcome of MD simulations of ice formation in the case of kaolinite. In fact, we have been able to reproduce the results of the MD simulations of ice formation on kaolinite reported in Ref. 41 (at the same supercooling and employing the same water model) using a number of different computational setups, as discussed in the supplementary material. In contrast, we argue that the tiniest structural details of the surface are crucial in determining the kinetics of ice formation on crystalline surfaces within atomistic simulations, and that surface relaxation can play an even more relevant role than the flexibility of the surface in promoting the formation of ice on this particular kaolinite surface. At this stage, it is reasonable to speculate that both flexibility and surface relaxation are likely to be a general issue when dealing with MD simulations of heterogeneous ice nucleation. It could be that the structural details of a particular substrate will determine which one of the two would be the dominant factor in ruling the kinetics of ice formation.

Before ending, we briefly comment on the force field models employed. Many options are available to simulate water.^16,23,79,80 The coarse grained mW model²³ is computationally very fast and as such it has been extensively used to model heterogeneous ice nucleation on, e.g., carbonaceous particles^25–27 and Lennard-Jones crystals.²⁴ However, with coarse grained approaches, a truly atomistic description of the nucleation mechanism cannot be achieved, and more importantly it is difficult to describe the interaction between water and a complex material. This is why here we have employed the atomistic TIP4P/Ice rigid model⁴⁹ for water in this work. This model reproduces many structural and dynamical properties of liquid water as well as of different ice phases correctly,⁴⁹ and it has been recently used to obtain an accurate reference for the thermodynamics and kinetics of homogeneous freezing⁸¹ at 230 K, corresponding to a supercooling $ΔT$ = 42 K (TIP4P/Ice water melts at 270 K⁸²). At this supercooling, the dynamics of the water network is far from being homogeneous^83,84 and special care has to be taken to correctly reproduce quantities like the self-diffusion coefficient and the structural relaxation time, as detailed in the supplementary material.

The CLAY_FF force field⁴⁸ is widely used to model clays as well as the interaction between clays and water, including swelling properties⁸⁵ and confinement effects.⁸⁶ As we are interested in having a reliable description of the water-surface interface, we have investigated the extent of surface relaxation for two kaolinite surfaces customarily investigated in the context of heterogeneous ice nucleation: the hydroxylated (KAO_OH) and siloxane (KAO_Si) (001) faces (see Fig. 1(a)). As shown in Fig. 6(a), the arrangement of oxygen atoms at the surface, which is critical in templating ice formation,^33,41,44 can be described by the O–O–O angles $ϕ$ (for KAO_OH), $α$ and $β$ (for KAO_Si). We compare in Fig. 6(b) the values of $ϕ$⁠, $α$⁠, and $β$ at the experimental atomic positions (Exp., as obtained upon cleavage of the bulk crystal) with those obtained for the relaxed configurations of KAO_OH and KAO_Si, calculated by DFT and CLAY_FF. In the case of the KAO_OH face, $ϕ$ changes by ∼3° for both DFT and CLAY_FF. However, DFT and CLAY_FF simulations give substantially different results for the KAO_Si face: DFT predicts a marginal increase of the asymmetric buckling between oxygen and silicon atoms (⁠$α≃ 140$ and $β≃ 100$⁠), as depicted in the inset of Fig. 6(b); in contrast, the CLAY_FF relaxed structure is almost perfectly symmetric (⁠$α≃β≃ 120$⁠), the buckling is absent and the atoms at the surface form regular hexagonal patterns. Note that this symmetric arrangement of oxygen atoms at the KAO_Si surface has also been predicted by the CLAY_FF force field for materials such as mica,⁸⁷ while DFT calculations on the same system⁸⁸ resulted in buckled arrangements, in line with what we have observed here. We remark that different starting structures and/or different DFT exchange-correlation functionals do not affect these findings, as reported in the supplementary material. Whether the actual structure of the KAO_Si face upon surface relaxation is closer to the CLAY_FF or the DFT prediction remains to be seen. This is one reason why here we have limited the discussion to the KAO_OH (001) surface.

We have investigated several aspects of atomistic simulations of heterogeneous ice nucleation on a realistic surface, choosing as an example the well-characterized (001) surface of kaolinite, a prototypical clay mineral of relevance to ice formation in the atmosphere.

Previous MD simulations^33,41 suggest that ice nucleation occurs on the hydroxylated (001) kaolinite surface via the formation of hexagonal ice, due to the favorable interaction between its prism face and the hexagonal arrangement of the hydroxyl groups at the surface of the clay. Here, we have established the preference of the surface for hexagonal ice over the cubic polytype by means of seeded MD simulations using a fully flexible model of kaolinite. We find that nuclei of cubic ice exposing the basal face to the clay are not stabilized by the presence of the surface and that these nuclei therefore tend to shrink back into the liquid phase. On the other hand, nuclei of hexagonal ice substantially wet the kaolinite surface and proceed to grow. We have estimated the critical nucleus size for these nuclei of hexagonal ice to be roughly two times smaller than what has been reported for homogeneous water freezing at the same supercooling.

We have also verified by looking at the natural fluctuations of the water network at a flexible water-kaolinite interface using a very long (μs) unbiased MD simulation that the overwhelming majority of pre-critical ice nuclei that form on top of the clay are indeed made of hexagonal ice exposing the prism face to the surface. This demonstrates that such nuclei spontaneously form at the surface, and that indeed the hexagonal polytype is the only one involved in the early stages of heterogeneous ice nucleation on this particular kaolinite surface. However, we note that a number of different nucleation sites such as surface defects or functional groups on the edges of the mineral surface can all, in principle, contribute to ice formation on kaolinite (and indeed on many other clay minerals and related materials). Hence the need for the experimental characterization of these sites where nucleation is promoted: such a knowledge will allow computer simulations to address comprehensive nucleation scenarios, eventually leading to more accurate prediction of the ice nucleating ability of these materials.

As discussed, it has recently been reported⁴¹ that the flexibility of the clay surface influences the rate of ice formation. In this work, we find that surface relaxation can be equally important. In particular, small changes in the structure of the hydroxylated (001) kaolinite face drastically alter the free energy cost needed to form an extended hexagonal motif of ice-like molecules at the water-kaolinite interface. The occurrence of this templating layer leads to the formation of ice within ∼100 ns if the atoms of the kaolinite surface are frozen at the experimental positions of the bulk phase. However, upon surface relaxation the free energy cost for creating such a template is much higher and nucleation does not take place on the μs time scale. Note that of all the water-kaolinite interfaces considered in this work, S3 is arguably the best representation of the system, as actual surfaces are not frozen and/or unrelaxed.

In addition, we note that the CLAY_FF force field, customarily used to model clays as well as water interacting with clays, seems to provide a reliable description of ice nucleation at the water-KAO_OH interface. However, this classical force field predicts, in the case of the siloxane (001) surface of kaolinite, a surface relaxation which is not consistent with the outcome of DFT calculations. Thus, we argue that at this stage it is not clear whether the formation of ice on this particular surface can be safely modeled using the CLAY_FF force field.

We also remark that the fact that the nucleation process is so sensitive to the structure of the interface strongly suggests that future efforts should be devoted to produce more accurate interatomic force fields for water at complex interfaces. In fact, we have seen that subtle effects such as surface relaxation can truly affect the nucleation kinetics to a point where it becomes rather difficult to benchmark computational results and most importantly to compare them with experimental data.

See supplementary material for simulations, dynamical properties of the supercooled water network, the DFT calculations of kaolinite surface relaxation, the computational geometries used to model the water-kaolinite interface and how they affect the formation of ice, the formulation of the order parameter used to pinpoint ice nuclei, the ice formation on the interface S1, where structural relaxation is faster than ice nucleation, and the details of the metadynamics simulations used to generate the ice nuclei for seeded molecular dynamics.

This work was supported by the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 616121 (HeteroIce project). A.M. is also supported by a Royal Society Wolfson Research Merit Award. The work of A.M. and A.Z. has also been sponsored by the Air Force Office of Scientific Research, Air Force Material Command, USAF, under Grant No. FA8655-12-1-2099. The authors acknowledge the use of the UCL Grace and Legion High Performance Computing Facilities, the use of Emerald, a GPU-accelerated High Performance Computer, made available by the Science & Engineering South Consortium operated in partnership with the STFC Rutherford-Appleton Laboratory and the use of ARCHER UK National Supercomputing Service (http://www.archer.ac.uk) through the Materials Chemistry Consortium through the EPSRC Grant No. EP/L000202.