Molecular simulations of heterogeneous ice nucleation. II. Peeling back the layers

Coarse grained molecular dynamics simulations are presented in which the sensitivity of the ice nucleation rate to the hydrophilicity of a graphene nanoflake is investigated. We find that an optimal interaction strength for promoting ice nucleation exists, which coincides with that found previously for a face centered cubic (111) surface. We further investigate the role that the layering of interfacial water plays in heterogeneous ice nucleation and demonstrate that the extent of layering is not a good indicator of ice nucleating ability for all surfaces. Our results suggest that to be an efficient ice nucleating agent, a surface should not bind water too strongly if it is able to accommodate high coverages of water.


I. INTRODUCTION
As liquid water is cooled below its melting point it crystallizes to solid ice.This familiar yet important process is not fully understood, especially at the molecular level.It is known that pure water can exist in the liquid state far below 0 • C and that the reason we see ice formation occur at temperatures above approximately −35 • C is due to the presence of impurity particles. 1This is known as heterogeneous nucleation.It is also known that different impurity particles aid ice formation with different efficiencies, for example, feldspar mineral particles have been found to be better ice nucleating agents than clay mineral particles. 2What is severely lacking, however, is a comprehensive understanding of heterogeneous ice nucleation: we simply do not understand which properties of a material affect its ability to nucleate ice.Given the ubiquity of ice formation, and its important role in the atmospheric, geological and biological sciences, as well as the problems it can cause in the food, transport and energy industries, acquiring a full understanding of heterogeneous ice nucleation remains a major challenge in urgent need of address. 31][12][13][14][15] Recent studies have used molecular dynamics simulation in combination with coarse grained models to probe the mechanisms of heterogeneous ice nucleation.In particular, Lupi et al. 13,14 have investigated the effect of graphitic surfaces on ice nucleation, and have found that the extent of layering of interfacial water correlates with the freezing temperature; our previous work, 12 on the other hand, focused a) Electronic mail: angelos.michaelides@ucl.ac.uk on how the hydrophilicity of an hexagonal surface that acts as a template for the basal face of ice affects the rate, and found that an optimal interaction strength between the water and the surface exists.In this article, we present further results from simulations of ice nucleation in the presence of 'graphitic' surfaces.Unlike Ref. 14, where the primary aim was to understand ice nucleation on soot particles, here we are not attempting to model actual graphitic surfaces.Rather, the aim of this study is to exploit the smoothness of the potential experienced by the water molecules at such surfaces, and compare to the results obtained in Ref. 12, where the hexagonal surface under investigation presented distinct adsorption sites for the interfacial water molecules.We therefore probe a far greater range of hydrophilicities of these 'graphitic' surfaces than previously sought by Lupi et al.In what follows, we will find that the graphitic surfaces also exhibit an optimal interaction strength with water for promoting ice nucleation, which coincides with the optimal interaction strength found for the hexagonal surfaces presented in Ref. 12. We will also see that the previously suggested layering mechanism 13,14 requires slight modification to be applied to strongly adsorbing surfaces and that the in-plane structure of the interfacial water molecules can affect the layering mechanism.This suggests that the layering mechanism cannot be used to explain the ice nucleating ability of surfaces in general.Finally, we discuss the origin of the observed optimal interaction strength and suggest a rule-of-thumb for relating the surface hydrophilicity to the ice nucleating efficiency.

A. Systems and force fields
We have investigated heterogeneous ice nucleation in the presence of a rigid graphene nanoflake (GNF) of varying hydrophilicity, which is totally immersed in water as shown in Fig. 1.The interaction of water with the GNF was modeled using the two-body part of the Stillinger-Weber (SW) potential in the same manner as Lupi et al. 13,14 The GNF consisted of 217 carbon atoms, with a carbon-carbon bond-length of 0.142 nm (perfect edges were assumed).As in Refs.13 and 14, we used σ SW = 0.32 nm to define the range of the water-carbon interaction (the functional form of the SW/mW potential is given elsewhere 16 ).The interaction strength was tuned by varying SW .We note here that for the graphitic surfaces in Ref. 14, values in the range 0.12 ≤ SW ≤ 0.2 kcal/mol were investigated; as this work is concerned with trying to obtain general understanding rather than modeling a specific system, we have broadened this range to 0.06 ≤ SW ≤ 1.5 kcal/mol.The total energy after geometry optimization of a single water molecule at the center of the GNF was used to define E ads , and no interaction was defined between the carbon atoms as their equations of motion were not integrated.

B. Simulation settings
All simulations were performed using the LAMMPS simulation package 17 and the coarse grained mW model for water. 16The velocity Verlet algorithm was used to propagate the equations of motion of the water molecules, using a 10 fs time step.In all simulations, 2944 mW molecules were used and periodic boundary conditions were applied in all three dimensions.Temperature and pressure were maintained using the Nosé-Hoover thermostat and barostat (with a chain length of 10) with relaxation times of 1 ps and 2 ps respectively.A 100 ns trajectory was first performed at 290 K and 1 bar, from which initial configurations were drawn (different initial configurations were separated by at least 5 ns in the high temperature trajectory).At the start of the nucleation simulations, velocities for the water molecules were drawn randomly from a Maxwell-Boltzmann distribution to give an initial temperature of 205 K. Simulations were stopped after 500 ns if nucleation did not occur.To detect 'ice-like' molecules, we have used the CHILL algorithm of Moore et al. 18 Rates were extracted in the same manner as Ref. 12 and are directly comparable, since we have used the same simulation protocol.

C. Analysis
For each value of E ads , an extra simulation of 10 ns was performed at 215 K and 1 bar (following a 1 ns equilibration period from a 290 K configuration) such that sufficient statistics in the liquid state could be obtained.Similarly, a set of 10 ns simulations were performed at 225 K and 1 bar for the face centered cubic nanoparticle investigated in Ref. 12 (we refer to this nanoparticle as the 'FCC-111 NP').
The layering of interfacial water was computed as: where ρ(z) is the local water density at a height z above the surface (see Fig. 2), ρ bulk is the density of bulk liquid water at 215 K (or 225 K) and 1 bar (also obtained from a 10 ns simulation) and z bulk = 1.8 nm is a height at which ρ(z) → ρ bulk .The integration was performed using the Trapezium rule (Simpson's rule was also used, with a maximum discrepancy between the two methods of 3%, and agreement generally within 1%).In computing ρ(z), only water molecules in the column above the surface were considered. 19e probability density P (x, y) of water molecules in the plane of the surface was computed for the water molecules in the contact and second layers above the surface.A water molecule was defined as being in the contact layer if 0 ≤ z < 0.45 nm and in the second layer if 0.45 ≤ z < 0.8 nm as shown in Fig. 2 (a).A similar analysis was also performed for the FCC-111 NP, with water molecules defined as being in the contact layer if 0 ≤ z < 0.35 nm and in the second layer if 0.35 ≤ z < 0.7 nm, as shown in Fig. 2

A. Nucleation rates and the role of interfacial layering
In Fig. 3 we show how the nucleation rate R varies with the hydrophilicity of the GNF.Specifically, we have plotted log 10 (R/R hom ) against E ads /∆H vap , where R hom is the homogeneous nucleation rate and ∆H vap is the heat of vaporization of bulk mW water (10.65 kcal/mol at 298 K). 16 We can clearly see that for the weakest interaction strength, the GNF has little effect on the nucleation rate.As E ads increases, the rate rapidly increases to reach a maximum at E ads /∆H vap ≈ 0.3−0.4 that is approximately a factor 25 faster than homogeneous nucleation.Upon increasing E ads further, the rate steadily drops until E ads /∆H vap ≈ 1.0 when the rate begins to steadily increase again.For the most strongly interacting GNF investigated, the rate is a little over 10 times that of homogeneous nucleation.We also show the results from Ref. 12 of the nucleation rate in the presence of the FCC-111 NP.Both the GNF and the FCC-111 NP exhibit a maximum in rate at E ads /∆H vap ≈ 0.3−0.4,but differ in that at E ads /∆H vap ≈ 1.0, the FCC-111 NP crosses from promoting to inhibiting rather than exhibiting a local minimum like the GNF.This crossover from promoting to inhibiting can be explained by an excess of favorable adsorption sites at the FCC-111 NP (shown in Fig. 4 and discussed in detail in Ref. 12).
To explain the ice nucleating ability of the graphitic surfaces such as those considered here, Lupi et al. found that the layering of interfacial water L (see Eq. 1) correlates with the ice nucleating ability. 14In Fig. 5 (a), we show the dependence of L on E ads for the GNF.Unsurprisingly, L increases monotonically with E ads and we therefore cannot explain the observed non-monotonic dependence of R on E ads seen in Fig. 3 simply by the extent of layering.Figs. 6 (a) and 6 (b) provide some insight into why this is the case, where we show typical structures of water in contact with the GNF at 215 K for E ads /∆H vap ≈ 0.25 and E ads /∆H vap ≈ 1.9, respectively.For values of E ads that yield the highest rates, such as in Fig. 6 (a), water forms structures in the contact layer that  resemble the hexagonal structure of ice.Indeed, when ice formation is observed at 205 K, it appears to be driven by the formation of such hexagonal patches in the contact layer, consistent with previous studies. 13,14In the case of the strongly adsorbing GNF shown in Fig. 6 (b), it is clear that the number of water molecules in contact with the GNF has increased and that, rather than an hexagonal structure similar to ice, a structure consisting predominantly of smaller membered rings is now observed.This structure bears a strong resemblance to those seen in confined water layers under high pressures 20,21 and can be understood as water maximizing its interaction to the surface with only a slight cost in hydrogen bond energy. 10hen ice formation is observed at this surface, it does so in the layers of water above the contact layer, with the structure in the contact layer remaining unchanged.
Thus for high values of E ads , the contact layer is 'inactive' with respect to nucleation.The observation that the contact layer becomes inactive to ice nucleation for strong adsorption energies is enough to understand why we begin to see a decrease in the nucleation rate beyond E ads /∆H vap ≈ 0.3−0.4.To explain the increase in rate beyond E ads /∆H vap ≈ 1.0, however, requires further analysis.To this end, we have computed the layering of interfacial water with contributions from the first layer excluded: where z 0 = 0.45 nm.In Fig. 5 (a) we can see that L * , like L, also increases monotonically with E ads , but much more slowly.We can also see that the value of L * at E ads /∆H vap ≈ 1.9 is similar to the value of L at E ads /∆H vap ≈ 0.2 and from Fig. 3, that these two adsorption energies yield similar rates (both approxi-mately a factor 10 faster than homogeneous nucleation).
It therefore seems that beyond E ads /∆H vap ≈ 1.0, the extent of layering in the second layer of water and above becomes sufficient to promote ice nucleation.This can be seen more clearly in the inset of Fig. 5 (a); bearing in mind that the most weakly interacting GNF yields L ≈ 3 and does not promote ice nucleation, we can see that L *  GNF with E ads /∆Hvap ≈ 1.9.At this more strongly adsorbing surface, the structure in the contact layer consists predominantly of smaller membered rings.This persists even after ice nucleation at 205 K, which at this interaction strength, occurs in the water layers above.
only begins to exceed this value for E ads /∆H vap > 0.7.

B. The layering mechanism depends upon the in-plane structure of water
The previous section provides strong evidence in support of the layering mechanism, albeit with a slight modification to what was originally proposed. 13,14Conceptually, the layering mechanism is appealing and can perhaps be understood in terms of reducing the entropic barrier to nucleation: if water molecules are restricted to motion in a particular plane (with a density acceptable for ice nucleation), then the space that they can explore is effectively reduced by one dimension relative to the bulk liquid, which subsequently reduces the number of possible configurations that the water molecules can explore.This argument, however, implicitly assumes that the effective potential experienced by the water molecules within a layer is uniform.
For the graphitic surfaces investigated in this study, the assumption of a uniform effective potential within the layers is likely to be reasonable.Now consider the FCC-111 NP with E ads /∆H vap ≈ 1.9 which, as can be seen in Fig. 3, does not promote ice nucleation.The GNF with similar E ads , on the other hand, enhances ice nucleation by a factor ∼10 relative to homogeneous nucleation.Furthermore, the layering (excluding contributions from the contact layer) above both of these surfaces is similar, with L * ≈ 6 for the FCC-111 NP and L * ≈ 5 for the GNF, as shown in Fig. 5 (b) (a value of z 0 = 0.35 nm is used in Eq. 2 for the FCC-111 NP).The difference in rates can be understood in terms of in-plane structure, as seen in Fig. 7, where we show − ln[P (x, y)], where P (x, y) is the probability density of water molecules in the plane of the surface at 215 K, both for the contact layer and the second layer above the surface (see Sec. II C).At the FCC-111 NP, the water molecules bind at distinct adsorption sites (see Fig. 7 (b)) and do not diffuse over the 10 ns timescale of the simulation.Importantly, this impacts upon the structure of the water molecules in the second layer (see Fig. 7 (d)), which tend to be found directly above those in the contact layer.In contrast, the water molecules in contact with the GNF (see Fig. 7 (a)) do not adsorb at particular adsorption sites and are not immobile like at the FCC-111 NP, resulting in a smearingout of P (x, y).Accordingly, P (x, y) for the second layer above the GNF (see Fig. 7 (c)) is much smoother than at the FCC-111 NP.The smoothness of P (x, y) for the second layer above the GNF means that the water molecules can rearrange to form ice-like structures, whereas the corrugation in P (x, y) for the second layer above the FCC-111 NP appears to frustrate the water molecules, hindering ice formation.
In Refs.13 and 14, the general applicability of the layering mechanism was left as an open question.The results of our simulations show that the extent of layering can be important for ice nucleation, but that if the coverage of water at the surface is too high, then the contributions of the contact layer to the layering should be omitted.Furthermore, the importance of layering is dependent upon the water molecules experiencing a rather uniform effective potential within the layers.The extent of layering is therefore not a universal descriptor for the ice nucleating ability of surfaces.
C. The effect of surface hydrophilicity on heterogeneous ice nucleation In this section, we will discuss the observation that the GNF and the FCC-111 NP both have optimal interaction strengths at E ads /∆H vap ≈ 0.3−0.4 in more detail.If the interaction between the water and the surface is too weak, the induced structural differences from the bulk liquid are not significant enough to promote ice nucleation at either the GNF or the FCC-111 NP.What is more intriguing is why the ice nucleation rate at both surfaces decreases beyond E ads /∆H vap ≈ 0.3−0.4.Despite the differences in surface topography (the GNF can be considered a smooth surface whereas the FCC-111 NP presents distinct adsorption sites), both surfaces share a common feature: they can both accommodate water coverages that are higher than that when ice forms at the surface.This has been demonstrated qualitatively in Figs. 4 and 6.Fig. 8 provides quantitative evidence for this statement, where we show the lateral density of water molecules σ in the contact and second layers: where the integral runs over the layer of interest (see Sec. II C).Note that σ can also be computed explicitly by counting the average number of water molecules in a given layer: such an approach also gives information regarding the fluctuations of σ.For the GNF, shown in Fig. 8 (b), we can see that σ for the contact layer steadily increases with E ads whereas for the second layer, although increasing slightly initially, σ is essentially constant.In terms of layering, this means L * is increasing primarily through a narrowing of the second peak in ρ(z), whereas L also has a contribution from an increased number of water molecules at the surface.To a lesser extent the same is true for the FCC-111 NP, shown in Fig. 8 (b), although some variation in σ for the second layer is also observed.
We have also marked in Figs. 8 (a) and (b) the water coverage of ice that forms at the surface for the GNF and the FCC-111 NP, respectively. 22We can clearly see that for E ads /∆H vap ≈ 0.3−0.4,the value of σ in the contact layer is comparable to that of ice that forms at the surface.Thus, as E ads increases further, it becomes increasingly favorable for water to adsorb to the surface, to the detriment of ice nucleation occurring in the contact layer.We therefore suggest a rule-of-thumb for the role of surface hydrophilicity in ice nucleation: for surfaces that can accommodate water coverages higher than that required by ice, binding to the surface should not be too strong, with optimal adsorption energies approximately 30-40% the heat of vaporization of liquid water.We must stress that the importance of surface hydrophilicity will depend upon the water coverage that the surface can accommodate.As we have shown in Ref. 12, surfaces with E ads /∆H vap 0.4 can exhibit excellent ice nucleating efficiency, provided that the coverage of water at the surface does not exceed that of ice.

IV. CONCLUSIONS
We have presented computer simulations of heterogeneous ice nucleation in the presence of a graphene nanoflake immersed in water and investigated how its hydrophilicity affects the nucleation rate.The results of our simulations in part support the previously proposed layering mechanism of Lupi et al., 13,14 although we have seen that for strongly adsorbing surfaces, the increased layering, due to a higher coverage of water molecules, is detrimental to ice nucleation.Under such conditions, by excluding the contribution of the water molecules in contact with the surface, the extent of layering is, however, still found to correlate with the heterogeneous nucleation rate.It has also been demonstrated that the layering mechanism is not universal, as surfaces that exhibit similar degrees of interfacial layering can yield vastly different rates.We attribute this finding to the extent that the surface affects the in-plane structure of the water molecules: for surfaces where the water molecules move in a smooth effective potential, like the graphitic surfaces investigated in this article, the extent of layering describes the nucleation rate well; for surfaces that present distinct adsorption sites, such as the FCC (111) surface investigated in Ref. 12, the induced structure can frustrate ice nucleation not only in the contact layer, but also in the layer above.
We have observed that an optimal interaction strength between the graphene nanoflake and water for ice nucleation exists, and that this coincides with the optimal interaction strength found for the FCC (111) surface investigated previously. 12This behavior has been rationalized by noting that both surfaces are able to accommodate coverages of water that are higher than that when ice forms at these surfaces.We have proposed a ruleof-thumb, which states that in order to nucleate ice efficiently, a surface should not bind water too strongly if it can accommodate high coverages of water.Such insight may prove useful when trying to predict a material's ice nucleating ability, especially as the coverage of liquid water should be obtainable through e.g.surface X-ray diffraction experiments.

FIG. 1 .
FIG. 1. (color online) A typical simulation of heterogeneous ice nucleation.The GNF is shown by large silver spheres, and ice-like molecules are shown by blue lines.The remaining liquid-like water molecules are shown by small gray spheres.The box shows the boundary of the simulation cell and periodic boundary conditions are used throughout. (b).

3 )FIG. 2 .
FIG. 2. (color online) Density profile ρ(z) of water above the surface of the (a) GNF at 215 K and (b) the FCC-111 NP at 225 K, for different values of E ads /∆Hvap.At both surfaces water forms layers, with the intensity and sharpness of the layers increasing with E ads .The dark gray shaded region indicates water molecules defined as belonging to the first layer, and the light gray shaded region water molecules defined as belonging to the second layer (see Sec. II C).

FIG. 3 .
FIG. 3. (color online)Dependence of the heterogeneous nucleation rate on surface hydrophilicity.As E ads increases so too does the hydrophilicity.The homogeneous and FCC-111 data are taken from Ref. 12. Like the FCC-111 NP, the GNF also exhibits a maximum rate at E ads /∆Hvap ≈ 0.3−0.4,but in contrast, exhibits a local minimum in the rate at E ads /∆Hvap ≈ 1.0.

FIG. 4 .
FIG. 4. (color online) Typical structures that form in the contact layer at the FCC-111 NP at 225 K.The FCC-111 NP is shown in silver and the water molecules in blue.(a) FCC-111 NP with E ads /∆Hvap ≈ 0.3.An hexagonal overlayer, which is commensurate with the surface, and resembles ice is observed and facilitates ice formation at 205 K. Unoccupied, or 'excess', adsorption sites (highlighted by yellow circles) are present when this structure forms.(b) FCC-111 NP with E ads /∆Hvap ≈ 1.9.For this stronger interaction with the surface, the excess sites are occupied and the hexagonal structure resembling ice is no longer observed.

FIG. 5 .
FIG. 5. (color online)Dependence of the extent of layering of interfacial water on E ads .The black squares show the layering over the whole density profile L as defined by Eq. 1, whereas the red circles show the layering excluding the contact layer L * as defined by Eq. 2. (a) Result for the GNF at 215 K.Both L and L * increase monotonically with E ads , but L does so much more rapidly.The inset contains the same data but with a smaller scale for the y-axis, which shows more clearly that L * > ∼ 3 only once E ads /∆Hvap > ∼ 0.7 (when E ads /∆Hvap ≈ 0.1, L ≈ 3 and nucleation is not enhanced).(b) Results from the FCC-111 NP at 225 K.Although L is much greater at the FCC-111 NP than at the GNF, the values of L * are comparable.

FIG. 6 .
FIG. 6. (color online) Typical structures of water in the contact layer at the GNF at 215 K.The GNF is shown in silver and the water molecules in blue.(a) GNF with E ads /∆Hvap ≈ 0.25.Patches of ice-like hexagons readily form in the contact layer and facilitate ice formation at 205 K. (b)GNF with E ads /∆Hvap ≈ 1.9.At this more strongly adsorbing surface, the structure in the contact layer consists predominantly of smaller membered rings.This persists even after ice nucleation at 205 K, which at this interaction strength, occurs in the water layers above.

FIG. 7 .
FIG. 7. (color online) In-plane distribution of water molecules above the surface, plotted as − ln[P (x, y)] (see Sec. II C) with E ads /∆Hvap ≈ 1.9 at 215 K. (a) and (b) show the contact layer at the GNF and FCC-111 NP, respectively.(c) and (d)show the second layer above the GNF and FCC-111 NP, respectively.Unlike at the GNF, the water molecules in the contact layer with the FCC-111 NP are bound at specific adsorption sites and do not diffuse over the timescale of the simulation (10 ns).The water molecules in the second layer above the FCC-111 NP consequently exhibit greater structure than those in the second layer above the GNF.

2 )E
FIG. 8. (color online) Dependence of the lateral density of water molecules σ on E ads for (a) the GNF at 215 K and (b) the FCC-111 NP at 225 K.The blue dashed lines indicate the water coverage of ice that forms at both surfaces (the blue shaded area in (a) indicates the standard deviation of a sample average, taken over the range 0.16 ≤ E ads /∆Hvap ≤ 0.56).The black and red shaded areas indicate the standard deviation in the measured values of σ for the contact and second layers, respectively.When E ads /∆Hvap ≈ 0.3−0.4,the coverage in the contact layer is close to the water coverage of ice that forms.