Roles of Surface Energy and Temperature in Heterogeneous Ice

Article pubs.acs.org/JPCC

Roles of Surface Energy and Temperature in Heterogeneous Ice Nucleation Chu Li, Xiang Gao, and Zhigang Li* Department of Mechanical and Aerospace Engineering, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong

Downloaded via UNIV OF SOUTH DAKOTA on September 23, 2018 at 07:05:39 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Heterogeneous ice nucleation (HIN) is strongly related to the dynamics of hydrogen bonds in water at an interface. In this work, we investigate the microscopic kinetics of HIN through molecular dynamics simulations. The dynamics of hydrogen bond network (HBN) at interfaces is studied under the coupled effects of thermal fluctuation and water−surface molecular interactions. It is revealed that the lasting time of the HBN at the interface is critical to HIN. Under comparable thermal and surface effects, which result in a proper lasting time of the HBN, HIN is promoted. However, if the thermal effect or the surface effect dominates over the other, the lasting time of the HBN at the interface would be either too long or too short, leading to the failure of HIN. By varying the water−surface interaction strength, i.e., binding energy, and temperature, a diagram of HIN events is presented, which shows that HIN is only favored in certain temperature and surface energy ranges.

■

INTRODUCTION Ice can form either homogeneously in bulk water or heterogeneously at the surface of a foreign material. Heterogeneous ice nucleation (HIN) is the dominant mode of water crystallization in nature and is of great importance in a variety of areas, such as atmospheric physics,1,2 cryobiology,3 transportation,4 and the operation of common infrastructures.5 In many applications, a good control of ice nucleation under the influences of different factors (e.g., surface properties and temperature) is highly desired.6−11 Although HIN at the surface of different materials has been extensively explored in the literature,12−21 the roles of various parameters in HIN are still unclear. Previous experimental and numerical studies show that surface defects, such as cracks and cavities,22−24 may induce HIN through increasing the local density of water molecules. HIN can also be enhanced through surface functionalization (e.g., hydroxyl groups).25−28 In addition, the lattice mismatch between ice and the crystalline structure of the surface29−31 and the surface wettability greatly influence HIN.27,31,32 Unfortunately, how these factors fundamentally affect HIN is not well understood and some of previous results are even controversial. For instance, some work suggests that HIN can be boosted by increasing the surface wettability.33 Other studies, however, show that there is no direct relationship between surface wettability and ice nucleation ability.27,32 In addition, the findings about the role of ice-surface lattice mismatch in HIN are also inconsistent.29−31,34,35 Furthermore, the temperature effect, which is expected to be critical, and how temperature is coupled with other factors have not been systematically examined. Nevertheless, a clear picture about the kinetics of HIN requires not only the knowledge of classical thermody© 2017 American Chemical Society

namics but also a deep, microscopic understanding of ice nucleation processes at water−solid interfaces under the coupled effects of different factors. It is known that water molecules interact weakly through hydrogen bonds, which form an unstable hydrogen bond network (HBN)36 in the liquid state. Under the thermal fluctuations of water molecules, i.e., the effect of temperature, the HBN breaks and reforms constantly. This makes the average lasting time of the HBN quite short in bulk water. With the presence of a foreign material (e.g., a solid surface), the water−surface intermolecular interactions could greatly affect the dynamics and the stability of the HBN, in addition to the thermal effect. Under proper conditions, certain temperature and water−surface interaction strength, water molecules at the interface may reorganize and form hexagonal structure, and the consequent HBN there becomes stable, leading to the formation of ice. Therefore, the dynamics of the HBN at the interface, under the coupled effects of water−surface interaction and temperature, fundamentally determines water crystallization, which, to our knowledge, has not been fully studied. In this work, we investigate HIN at the molecular level and probe the kinetics of ice nucleation at interfaces through molecular dynamics (MD) simulations. Our attention is on how the coupling of surface energy (or wettability), which is described by the water−surface molecular binding energy εws, and temperature T affects the stability of the HBN at the interface and consequently influences HIN. It is found the Received: March 26, 2017 Revised: May 7, 2017 Published: May 10, 2017 11552

DOI: 10.1021/acs.jpcc.7b02848 J. Phys. Chem. C 2017, 121, 11552−11559

Article

The Journal of Physical Chemistry C

Figure 1. Simulation system and a case of ice nucleation. (a) A snapshot of a MD simulation system. Surface, oxygen, and hydrogen atoms are represented by gray, purple, and green spheres, respectively. (b) Illustration of the dipole vector δv⃗ of a water molecule. (c) Time variation of water potential energy and the number of ice molecules (the inset shows the distribution of the starting time of ice nucleation). (d−k) Snapshots showing the evolution of ice clusters at different times (ice molecules are denoted by purple spheres and hydrogen bonds in water are represented by pink dash lines).

lasting time of the HBN at the interface plays an important role in ice nucleation. To promote HIN, the lasting time needs to be proper, which can only be achieved at a certain range of temperatures and water−surface binding energies, out of which, the lasting time of the HBN is either too long or too short and the crystallization of water is hindered. By varying εws and T, a diagram for HIN events in the εws−T coordinates is provided.

potential is employed to consider nonelectrostatic forces between surface atoms and the oxygen atoms of water molecules,

METHODS Molecular Dynamics Simulation. Molecular dynamics simulations are performed using the LAMMPS packages37 with a time step of 2 fs. The simulation system contains a rigid surface of hexagonal wurtzite structure and a water film of 1560 molecules sitting on the prism face of the surface, as shown in Figure 1a. The lattice parameters for the surface are a = 4.519 Å and c = 7.357 Å. The TIP4P/Ice water model38 is used to describe molecular interactions among water molecules. The ∠HOH bond angle and OH bond length in water molecules are maintained by the SHAKE algorithm. The Lennard−Jones

where rij is the separation between oxygen atom i and surface atom j, rc = 12 Å is the cutoff distance, σ = 3.1668 Å and εws are the water−surface collision diameter and binding energy, respectively. The surface atoms are not charged and the surface energy or wettability are tuned by changing the value of εws.39,40 Simulations are performed in canonical, i.e. (N, V, T), ensembles. In all the simulations, the system is initially heated up to 300 K for 2 ns and then cooled down to a supercooled temperature in the following 1 ns. The temperature of the system is controlled by the Nosé−Hoover thermostat. The lengths of the system in the x, y, and z directions are 36.2, 36.8,

⎛⎛ ⎞12 ⎛ ⎞6 ⎞ σ σ Uij(rij) = 4εws⎜⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎟ , ⎜⎝ r ⎠ ⎝ rij ⎠ ⎟⎠ ⎝ ij

■

11553

rij < rc (1)


Article


Figure 2. Snapshots illustrating the dynamics of hydrogen bonds, which break and form continuously (gray and green spheres represent surface and hydrogen atoms, respectively. Oxygen atoms are purple/pink if they belong to water/ice molecules). (b) is 0.08 ps after (a). Hydrogen bonds are denoted by pink dash lines. Red arrows show the places where the hydrogen bonds form after 0.08 ps, while blue arrows indicate the locations of broken hydrogen bonds.

where δvi⃗ (t) is the unit dipole vector of water molecule i at time t, as shown in Figure 1b, t0 is a reference time, ⟨·⟩ denotes the ensemble average, and ψi(t0) is an indicator, which is equal to 1 or 0 depending on whether water molecule i is at the interface or not at t = t0 (the interface is defined as the layer of 4.3 Å in thickness above the surface). The stability of HBN is related to the decay rate of Cvv(t − t0), as will be discussed later.

and 70.0 Å, respectively, and periodic boundary conditions (PBCs) are applied in all the directions. A void space of ∼35 Å in the z direction above water is used to eliminate the unnecessary water−surface interactions caused by PBCs. The cutoff distance for the short-range Coulombic potential in the TIP4P/Ice water model is set to be 10 Å and the Particle− Particle Particle−Mesh (PPPM) method is applied to consider long-range interactions. The simulation time ranges from 0.3 to 0.6 μs depending on whether ice nucleation occurs or not. Totally, more than 40 μs simulations are conducted. Bond-Orientational Order Parameter. To differentiate ice from water, the bond-orientational order parameter41−43 q6(i, j) between a pair of molecules i and j is calculated, which is given as follows: q6(i , j) ≡

■

RESULTS AND DISCUSSION It is known that water can be cooled down to 235 K at 1 bar without freezing,45,46 which roughly offers a reference temperature for water crystallization. For εws = 4.1575 kJ mol−1, at which the water contact angle is about 25° at 300 K, ice nucleation is observed at T = 230 K. The evolutions of the number of ice molecules based on the bond-orientational order parameter and the water potential energy are depicted in Figure 1c. It is seen that the potential energy remains almost constant initially, for t < 150 ns, and significant ice crystallization is not observed. Meanwhile, incipient ice-like clusters of hexagonal structure at the interface form and break constantly, and the size of the ice clusters fluctuates, as indicated by the number change of ice molecules (blue curve) in Figure 1c and the snapshots of the ice clusters at different times in Figures 1d−h and Figures 1a−e (top view) in the Supporting Information. After this “quiescent” period, for t > 150 ns, the potential energy drops, which corresponds to the continuous growth of ice, as indicated by the number of ice molecules in Figure 1c and snapshots in Figures 1i−k and S1f. To confirm that the potential drop is associated with ice nucleation instead of relaxation, a series of independent simulations for εws = 4.1575 kJ mol−1 and T = 230 K have been conducted. As shown in the inset of Figure 1c, the starting time of ice formation (the time when water potential starts to drop) roughly follows the Poisson distribution, which agrees with the stochastic process of ice nucleation. As discussed in the introduction, the formation of ice is largely determined by the dynamics of the metastable HBN at the interface under the influences of surface and temperature. The dynamics of the HBN is affected by the thermal fluctuation of water molecules and the water−surface intermolecular interactions. The thermal effect is straightforward. At a high temperature, the strong thermal fluctuation can break hydrogen bonds easily and the lasting time of the HBN is short, which is associated with a fast decay rate of the orientation correlation function Cvv. The effects of the surface are two folds. On the

6 ∑m =−6 Y6m(i)·Y 6*m(j)

|∑m Y6m(i)||∑m Y6m(j)|

(2)

where Y6m(i) is the spherical harmonics and Y*6m(i) is its complex conjugate. For molecule i having NC nearest neighbors, the average value of q6(i, j) at time t is calculated as follows: q6̅ (i , t ) =

1 Nc

j = Nc

∑ q6(i , j , t ) j=1

(3)

For bulk water molecules, Nc = 4. For water molecules at the interface, the neighboring molecules are searched in the water layer of 4.3 Å in thickness above the surface and Nc = 3 is used. Whether a molecule is counted as an ice or water molecule depends on the value of the following: I(i , t ) = Θ(q6̅ (i , t ) − q6̅ 0) ·Θ(q6̅ (i , t − Δt ) − q6̅ 0)

(4)

where Θ is the step function and q06̅ is a gate value. For bulk water molecules, q06̅ = 0.45,26 while q06̅ = 0.4 for water molecules at the interface. Δt = 0.4 ps is the time step used in the algorithm. If I(i,t) = 1, then molecule i is an ice molecule at time t. If I(i,t) = 0, then molecule i is a water molecule. Orientation Correlation Function. To probe the stability of the HBN at the interface, the orientation correlation function of water molecules, Cvv(t − t0),44 is computed for different times t0. Cvv(t − t0) is defined as follows: Cvv(t − t0) = ⟨δvi⃗(t ) ·δvi⃗(t0) ·ψi(t0)⟩

(5) 11554


Article


Figure 3. Orientation distribution of water molecules at the interface during different periods of time for εws = 4.1575 kJ mol−1 and T = 230 K. (θ and φ are the polar and azimuthal angles of δv⃗).

one hand, the water−surface interactions tend to break hydrogen bonds in water and make the HBN unstable. On the other hand, they can cause water molecules to reorganize at the interface through position and orientation changes, which may lead to the formation of orderly hydrogen bonds (e.g., hexagonal structure). Whether ice nucleation occurs or not depends on the competition between the thermal and surface effects. For the case in Figure 1c, the thermal effect is comparable to the surface effect. At the beginning (constant potential part), the hydrogen bonds at the interface are unstable such that they frequently switch between “formation” and “breakdown”, as shown in Figures 2a,b. Through the continuous breaking and formation of the HBN, water molecules adjust their orientations and reorganize at the interface. Gradually, ice clusters with hexagonal structure form at the interface, as indicated by the distribution of the polar and azimuthal angles, θ and φ, of the water dipole vectors δv⃗ in Figure 3. For t < 100 ns, the directions of δv⃗ are rather disordered. As time evolves, δv⃗ tends to assume specific directions as water molecules adjust their orientations, as shown in Figure 3 (155−160 ns), which correspond to a hexagonal ice structure. In this process, the stability or the lasting time of the HBN is critical, which can be measured by the decay rate of the orientation correlation function Cvv(t,−,t0) for the molecules at the interface. As shown in Figure 4, the decay rate of Cvv decreases with increasing time. This indicates that the HBN becomes more stable as water molecules at the interface reorganize and eventually form hexagonal structures. For t0 ≥ 175 ns, the decay rate remains

Figure 4. Orientation correlation function of interfacial water molecules for εws = 4.1575 kJ mol−1 and T = 230 K at different t0.

almost constant, which corresponds to the complete formation of ice at the interface and the further growth of ice. The ice nucleation observed in the case of Figure 1c is the consequence of comparable thermal and surface effects. However, if the temperature T and the water−surface binding energy εws are changed, the competition between the two effects may be skewed and ice nucleation may not take place. By varying T and εws, systematic studies are performed and an ice nucleation diagram is given in the εws − T coordinates, as 11555


Article


Figure 5. Surface and temperature effects on HIN. (a) Diagram of ice nucleation in the εws − T coordinates. (b−e) Interfacial water structures in different regions of the diagram at t = 250 ns (b: εws = 4.1575 kJ mol−1 and T = 225 K in region II; (c) εws = 4.1575 kJ mol−1 and T = 250 K in region III; (d) εws = 2.4945 kJ mol−1 and T = 230 K in region IV; and (e) εws = 5.8205 kJ mol−1 and T = 230 K region V).

interactions, the thermal fluctuation is insufficient and cannot effectively break the HBN to allow hydrogen bonds to reorganize and form hexagonal structures. This is confirmed by the slow decay rate of Cvv as time evolves, as shown in Figure 6a. As a result, the stable while irregular structure (Figure 5b) of the interfacial layer, which produces large free energy barrier for HIN,47,48 is incapable of triggering the growth of hexagonal ice, as indicated by the orientation distribution of interfacial molecules in Figure 7a. In region III, however, the thermal effect is strong and can easily break the HBN at the interface, reducing the stability of the HBN, as indicated by the fast decay of Cvv in Figure 6b. In this case, water molecules at the interface are very dynamic, which makes it difficult to form ice, as indicated by the random water

depicted in Figure 5a, where each case represents at least three independent simulations. The purple dots in Figure 5a are the places HIN is observed and the corresponding shaded area is where ice nucleation is likely to occur. In the other areas, HIN is not observed. Nevertheless, the diagram can roughly be divided into five regions, as shown in Figure 5a. In region I (purple area), the thermal and temperature effects are comparable and the mechanism for ice nucleation is similar to the case in Figure 1c, as explained previously. Another example in this region, εws = 4.1575 kJ mol−1 and T = 240 K, is depicted in Figure S2. In region II, the temperature is low, and the thermal effect is relatively weak compared with the surface effect. As water molecules at the interface form HBN due to the water−surface 11556


Article


Figure 6. Orientation correlation function of interfacial water molecules at different water−surface binding energies and temperatures, where ice nucleation is not observed. (a) εws = 4.1575 kJ mol−1 and T = 225 K. (b) εws = 4.1575 kJ mol−1 and T = 250 K. (c) εws = 2.4945 kJ mol−1 and T = 230 K. (d) εws = 5.8205 kJ mol−1 and T = 230 K.

Figure 7. Orientation distribution of interfacial water molecules for the cases in Figure 6. (a) εws = 4.1575 kJ mol−1 and T = 225 K. (b) εws = 4.1575 kJ mol−1 and T = 250 K. (c) εws = 2.4945 kJ mol−1 and T = 230 K. (d) εws = 5.8205 kJ mol−1 and T = 230 K.

structure at the interface in Figure 5c and the water orientation distribution in Figure 7b. In region IV, the temperature is intermediate while the water−surface interaction is weak compared with the thermal effect. The surface attraction is not strong enough to capture and stabilize water molecules at the interface, as shown by the fast decay of Cvv in Figure 6c. Similar to region III, the HBN is unstable and is unable to form ice (Figures 5d and 7c). In region V, the surface effect dominates and exerts strong

constraints to water molecules at the interface. Water molecules cannot move freely and the HBN is stable once it forms, as suggested by Cvv in Figure 6d, where the decay rate of Cvv slows down with increasing time. Similar to region II, the stable HBN leads to the formation of mixed structures, which does not favor the formation of hexagonal ice,47,48 as shown in Figures 5e and 7d. In this region, the density of water at the interface is about 15% higher than that of hexagonal ice, which impedes ice formation.47,48 11557


Article


The Importance of Feldspar for Ice Nucleation by Mineral Dust in Mixed-Phase Clouds. Nature 2013, 498, 355−358. (2) Tabazadeh, A.; Djikaev, Y. S.; Reiss, H. Surface Crystallization of Supercooled Water in Clouds. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 15873−15878. (3) Pandey, R.; et al. Ice-Nucleating Bacteria Control the Order and Dynamics of Interfacial Water. Sci. Adv. 2016, 2, e1501630. (4) Gent, R. W.; Dart, N. P.; Cansdale, J. T. Aircraft Icing. Philos. Trans. R. Soc., A 2000, 358, 2873−2911. (5) Stone, H. A. Ice-Phobic Surfaces That Are Wet. ACS Nano 2012, 6, 6536−6540. (6) Nagy, Z. K.; Braatz, R. D. Advances and New Directions in Crystallization Control. Annu. Rev. Chem. Biomol. Eng. 2012, 3, 55−75. (7) Kiani, H.; Sun, D. Water Crystallization and Its Importance to Freezing of Foods: A Review. Trends Food Sci. Technol. 2011, 22, 407− 426. (8) Kreder, M. J.; Alvarenga, J.; Kim, P.; Aizenberg, J. Design of AntiIcing Surfaces: Smooth, Textured or Slippery? Nat. Rev. Mater. 2016, 1, 15003. (9) Schutzius, T. M.; Jung, S.; Maitra, T.; Eberle, P.; Antonini, C.; Stamatopoulos, C.; Poulikakos, D. Physics of Icing and Rational Design of Surfaces with Extraordinary Icephobicity. Langmuir 2015, 31, 4807−4821. (10) Ehre, D.; Lavert, E.; Lahav, M.; Lubomirsky, I. Water Freezes Differently on Positively and Negatively Charged Surfaces of Pyroelectric Materials. Science 2010, 327, 672−675. (11) Kim, P.; Wong, T.; Alvarenga, J.; Kreder, M. J.; AdornoMartinez, W. E.; Aizenberg, J. Liquid-Infused Nanostructured Surfaces with Extreme Anti-Ice and Anti-Frost Performance. ACS Nano 2012, 6, 6569−6577. (12) Freedman, M. A. Potential Sites for Ice Nucleation on Aluminosilicate Clay Minerals and Related Materials. J. Phys. Chem. Lett. 2015, 6, 3850−3858. (13) Sosso, G. C.; Chen, J.; Cox, S. J.; Fitzner, M.; Pedevilla, P.; Zen, A.; Michaelides, A. Crystal Nucleation in Liquids: Open Questions and Future Challenges in Molecular Dynamics Simulations. Chem. Rev. 2016, 116, 7078−7116. (14) Murray, B. J.; O’Sullivan, D.; Atkinson, J. D.; Webb, M. E. Ice Nucleation by Particles Immersed in Supercooled Cloud Droplets. Chem. Soc. Rev. 2012, 41, 6519−6554. (15) Cabriolu, R.; Li, T. Ice Nucleation on Carbon Surface Supports the Classical Theory for Heterogeneous Nucleation. Phys. Rev. E 2015, 91, 052402. (16) Fraux, G.; Doye, J. P. K. Note: Heterogeneous Ice Nucleation on Silver-Iodide-Like Surfaces. J. Chem. Phys. 2014, 141, 216101. (17) Zielke, S. A.; Bertram, A. K.; Patey, G. N. A Molecular Mechanism of Ice Nucleation on Model AgI Surfaces. J. Phys. Chem. B 2015, 119, 9049−9055. (18) Lupi, L.; Hudait, A.; Molinero, V. Heterogeneous Nucleation of Ice on Carbon Surfaces. J. Am. Chem. Soc. 2014, 136, 3156−3164. (19) Glatz, B.; Sarupria, S. The Surface Charge Distribution Affects the Ice Nucleating Efficiency of Silver Iodide. J. Chem. Phys. 2016, 145, 211924. (20) Zielke, S. A.; Bertram, A. K.; Patey, G. N. Simulations of Ice Nucleation by Model AgI Disks and Plates. J. Phys. Chem. B 2016, 120, 2291−2299. (21) Zielke, S. A.; Bertram, A. K.; Patey, G. N. Simulations of Ice Nucleation by Kaolinite (001) with Rigid and Flexible Surfaces. J. Phys. Chem. B 2016, 120, 1726−1734. (22) Marcolli, C. Deposition Nucleation Viewed as Homogeneous or Immersion Freezing in Pores and Cavities. Atmos. Chem. Phys. 2014, 14, 2071−2104. (23) Christenson, H. K. Two-Step Crystal Nucleation via Capillary Condensation. CrystEngComm 2013, 15, 2030. (24) Kiselev, A.; Bachmann, F.; Pedevilla, P.; Cox, S. J.; Michaelides, A.; Gerthsen, D.; Leisner, T. Active Sites in Heterogeneous Ice Nucleation-the Example of K-Rich Feldspars. Science 2017, 355, 367− 371.

As the dynamics of HBN is essential to HIN, comparative simulations for εws = 4.1575 kJ mol−1 are also conducted using a coarse-grained water model (mW model18,27), for which ice forms immediately (

Roles of Surface Energy and Temperature in Heterogeneous Ice

Recommend Documents