Ice accumulation on solid surfaces is a severe problem for safety and functioning of a large variety of engineering systems, and its control is an enormous challenge that influences the safety and reliability of many technological applications. The use of molecular dynamics (MD) simulations is popular, but as ice nucleation is a rare event when compared to simulation timescales, the simulations need to be accelerated to force ice to form on a surface, which affects the accuracy and/or applicability of the results obtained. Here, we present an alternative seeded MD simulation approach, which reduces the computational cost while still ensuring accurate simulations of ice growth on surfaces. In addition, this approach enables, for the first time, brute-force all-atom water simulations of ice growth on surfaces unfavorable for nucleation within MD timescales. Using this approach, we investigate the effect of surface wettability and structure on ice growth in the crucial surface–ice interfacial region. Our main findings are that the surface structure can induce a flat or buckled overlayer to form within the liquid, and this transition is mediated by surface wettability. The first overlayer and the bulk ice compete to structure the intermediate water layers between them, the relative influence of which is traced using density heat maps and diffusivity measurements. This work provides new understanding on the role of the surface properties on the structure and dynamics of ice growth, and we also present a useful framework for future research on surface icing simulations.

Ice growth on surfaces can prevent the operation of—and cause damage to—a broad spectrum of man-made structures and devices, such as aircraft, ships, wind turbines, photovoltaic devices, heat exchangers, and telecommunications equipment.1–3 In each of these examples, liquid water comes into contact with an ice nucleating agent, which then initiates the process of freezing under the right temperature conditions. This heterogeneous ice nucleation process is the dominant mechanism of ice formation. Liquid water can also freeze in the absence of nucleating agents via homogeneous nucleation, although this is rare. In fact, pure liquid water has been shown to exist in a metastable supercooled state down to −37.5 °C4 and can be stored for months in a laboratory at −20 °C.5 Freezing of relevance to scientific and engineering applications typically occurs at higher temperatures, which means it occurs via heterogeneous nucleation. The wide scope of applications listed above has motivated the study of mechanisms of heterogeneous icing and development of techniques, which either prevent/delay icing or ensure that formed ice can be easily removed. Previous experiments have highlighted the importance of surface wettability and topography to surface icing.6–8 In addition, many research papers indicate how a small variation of the surface roughness could lead to a significant change of the icing mechanism,9,10 demonstrating the potential of surface-driven control of icing and the necessity for its systematic analysis. This raises an important question that is yet to be decisively answered: What is the role of the surface structure and chemistry in promoting/inhibiting icing?

There are two main challenges faced by experimentalists when trying to understand the mechanisms of surface icing. First, the analysis of experimental results is limited by the stochastic character of icing and the difficulties associated with the presence of water impurities and surface heterogeneities. Second, the mechanisms of icing occur in nano-to-microsecond timescales and length scales of a few hundred to thousand molecular diameters, which are very hard to observe experimentally. An alternative to running experiments is to use highly accurate molecular simulations using a well-validated tool such as molecular dynamics (MD). Since ice nucleation is fundamental to the initiation of icing, the majority of MD simulation studies have focused on the complex problem of the spatial and temporal evolution of heterogeneous ice nucleation. However, MD has its own set of challenges as individual atoms need to be tracked and interatomic interactions need to be computed at every time step, making it extremely computationally expensive to simulate very small systems (∼104 atoms) for very small periods of time (∼10−6 s). Additionally, given ice nucleation is a rare event, there is no guarantee that ice formation would occur in any given simulation, and many hundreds of simulations may be necessary, making this approach impractical. While there are some papers that successfully conducted brute-force simulations of surface icing using all-atom water models, they did so by focusing on specific types of surfaces that favored nucleation. This was achieved by either modifying the surface structure11,12 or adding surface charge13,14 in a way that induced nucleation to occur within the timescales of an MD simulation. For surfaces that are not necessarily tuned to promote ice nucleation in this manner, which we call unfavorable surfaces, running successful MD simulations of icing within the allowable simulation timescales remains a challenge.

There are two ways in which the MD community has bypassed these issues: using coarse-grained models (which reduce the interatomic interactions to be computed) or using enhanced sampling methods (which accelerate the freezing process). The majority of the MD literature uses coarse-grained models, such as the mW model,15 wherein the water molecule is approximated using a single interaction point. The use of the mW model, which is an order of magnitude faster than all-atom water models, such as TIP4P/Ice,16 has produced an extensive literature on the various factors affecting heterogeneous nucleation, such as surface wettability, surface crystal morphology, surface roughness, and interfacial water layering.17–23 Alternatively, observable nucleation for all-atom water models can be forced to occur within reasonable simulation timescales using enhanced sampling methods,24 including free energy methods25–27 and path sampling methods.28–30 

However, coarse-grained models and enhanced sampling techniques31,32 have their own drawbacks. Compared to all-atom models, the coarse-grained mW model provides incorrect values of interfacial free energy and nucleation rate and produces inaccurate ice growth rates that are four orders of magnitude higher.33 With enhanced sampling methods, the simulations are still computationally expensive and are therefore usually performed at high supercooling; this is suboptimal, as these low temperatures are not always representative of real-world applications.34 Additionally, the vast majority of enhanced simulation methods have only been used for homogeneous ice nucleation,35–37 and the literature on heterogeneous ice nucleation using enhanced sampling remains sparse. The use of favorable surfaces, the reduced accuracy, the reliance on artificially low temperatures, and the focus on homogeneous nucleation limit the real-world relevance of MD simulation results; as a consequence, the experimental and MD communities working on surface icing have largely operated independently of each other. However, given the important insights that MD simulations can provide to designing new surfaces that combat icing, it is important to be able to run accurate icing simulations on realistic surfaces not artificially designed to promote ice nucleation. This leads to a second important question: Can highly accurate all-atom water models be used to run MD simulations of icing on unfavorable surfaces within reasonable timescales?

Finally, in recent years, the structure and dynamics of the interfacial region, i.e., the exact location where the ice meets the surface, has gained importance in the heterogeneous icing community. Recent experiments have shown that a nanometer-thick quasi-liquid layer (QLL) can form at the surface–ice interfacial region.38 Other studies indicate that the presence of QLL leads to extremely low ice adhesion strength values, much lower than what surface-energy models predict.39,40 The structure and dynamics of QLL at the surface–ice interface is therefore important because it can influence ice formation and growth.41–44 Contrary to the studies of QLL at the ice–vapor interface, the studies of QLzL at the surface–ice interface are exclusively experimental. There have been no MD studies yet that investigate QLL formation and dynamics at the surface–ice interface. This leads to the third and final important question that we focus on here: What is the molecular influence of wettability and surface structure on the QLL at the surface–ice interface?

This paper presents a methodology to address the first question raised above in a way that avoids the difficulties encountered in conventional heterogeneous ice nucleation MD simulations. The proposed approach enables monitoring of the phase transition of liquid-water molecules into ordered ice structures with a brute-force approach, without the need for enhanced sampling, coarse-grained models, or surfaces favoring ice nucleation. This enables investigation into parameters that influence formation of ice on any surface for which accurate intermolecular potentials exist. Following that, a study into the effect of the surface wettability and lattice structure on ice growth is conducted, shedding light on the second question raised above. Finally, the approach can also be used to get molecular insight into the QLL structure and dynamics whenever such structures are observed in the interfacial region. This can be then linked back to the surface type and wettability, enabling insight into the third question raised above. This paper represents the first attempt at investigating these three open questions in surface icing using brute-force MD simulations. Note that, as with the alternatives discussed previously, there are limitations with this presented approach as well, such as the fact that the surfaces presented here are very simple in structure (and therefore not representative of the real world). However, it is hoped that this approach presents a useful alternative to the existing methods when studying heterogeneous ice growth using MD.

To overcome the limitations associated with coarse-graining or enhanced sampling, an alternative approach of “seeding” is used here. Seeding methods involve inserting an ice nucleus, or “seed,” into the supercooled liquid water at the start of the simulation.45 If the seed is ensured to be larger than the critical nucleus at that temperature, it will grow and ice growth will proceed as would be the case following homogeneous/heterogeneous nucleation.45 This approach is computationally very efficient and, unlike the alternative approaches, can produce ice at mild supercooling comparable to experiments. Note that seeded MD does not provide insight or information into the mechanism of nucleation itself, as it assumed that the water molecules have already formed a critical nucleus at the start of the simulation. Instead, the focus is to study the influence of the surface on formed ice adjacent to it, in particular, the first few layers of ice molecules that comprise the interfacial region.

Nada and Furukawa46,47 were the first to use seeded MD as described here but for ice–water interfacial simulations without a surface. They did not provide a complete description of the crystallization mechanism, since equilibrium was not reached due to the computational cost. The first systematic simulation of ice growth in the presence of water and vacuum was presented by Carignano, Shepson, and Szleifer.48 They placed an ice layer in contact with a water layer and observed the evolution of the system until equilibrium. A similar approach was also used recently to investigate the effect of freezing on a nanoparticle deposited on a surface.49 However, in contrast to these aforementioned studies where the objective was to simulate crystal growth in a setup with supercooled water placed in contact with ice, this work uses seeded MD to study, for the first time, ice growth on unfavorable surfaces. In this work, the supercooled liquid is sandwiched between the ice slab and the surface, which enables the study of ice growth at the surface, as shown in Fig. 1(a), a procedure we call “slab-seeding.”

FIG. 1.

(a) The MD setup for all slab-seeded ice growth simulations in this work. (b)–(d) MD time-evolution of the ice growth process near a solid FCC surface. (e) The mean molecular potential energy of water molecules during phase transition from liquid to solid. A plateau in the potential energy is reached when freezing is completed and steady state measurements are taken within the plateau. (f) Macroscopic contact angles of water droplets on different surfaces (FCC, BCC, and HCP) when varying the solid-water interaction potential ɛso. (g) A schematic representation of a flat overlayer and a buckled bilayer.

FIG. 1.

(a) The MD setup for all slab-seeded ice growth simulations in this work. (b)–(d) MD time-evolution of the ice growth process near a solid FCC surface. (e) The mean molecular potential energy of water molecules during phase transition from liquid to solid. A plateau in the potential energy is reached when freezing is completed and steady state measurements are taken within the plateau. (f) Macroscopic contact angles of water droplets on different surfaces (FCC, BCC, and HCP) when varying the solid-water interaction potential ɛso. (g) A schematic representation of a flat overlayer and a buckled bilayer.

Close modal

The size of the simulation box in the x, y, and z directions is 45.215, 31.326, and 90 Å, respectively, with periodic boundary conditions applied in all three directions. A vacant region of height 24 Å above the ice slab (not shown in the figure) is introduced in the computational domain in order to avoid the interaction of the ice slab with the surface through the top periodic boundary. This top region will be referred to as the “vapor” region in this paper. The ice slab consists of 80 hexagonal ice unit cells, with each unit cell composed of eight water molecules totaling 640 molecules (1920 atoms). The approximate thickness of the ice slab is 14.7 Å. The liquid region initially comprises 1632 water molecules (4896 atoms) with height 40 Å. This setup is designed to be large enough to avoid finite size effects while still being computationally tractable. Note that separate tests were run where the thickness of the ice slab and the vapor regions was increased, with a negligible impact on the observed results (details available in the supplementary material).

All ice growth simulations presented in this paper use the rigid TIP4P/Ice model, which is designed to reproduce the solid phase properties of water and has been shown to accurately predict phase diagrams and other physiochemical properties.16 The TIP4P/Ice model consists of four interaction sites, three of which are located at the oxygen (O) and hydrogen (H) atom positions of the standard H2O molecule. The fourth site, which is a massless site denoted by the symbol M, is coplanar with the oxygen and hydrogen atom positions and is placed at the bisector of the H–O–H angle. As with earlier TIP4P models,16 the O–H distance and H–O–H angle are fixed to the experimental values 0.9572 Å and 104.52°, respectively. The SHAKE algorithm50 is used to constrain the geometry of the water molecules. All intermolecular interactions between pairs of molecules for O, H, M, and the surface sites are based on the combined Lennard-Jones (LJ) and Coulomb potential. The total potential energy between two interaction sites i and j, Uij, is given by the sum of these two interactions,
Uij=4εijσijr12σijr6+14πε0qiqjr,
(1)
where rij is the distance between two interacting particles and ɛij, σij are the LJ parameters representing the depth of the potential well and the distance at which the potential energy is zero (often referred to as the diameter of the particle), respectively. In addition to the surface–surface and surface–oxygen interactions, a harmonic potential is applied to the surface to ensure that the surface atoms are constrained to their original positions.

Hexagonal ice is an anisotropic crystal, consisting of stacked hexagonal rings in the basal planes and the primary and secondary prismatic faces perpendicular to it. Ice growth in the basal plane occurs layer-by-layer, which is not the case for the prismatic plane where less organized ice growth is observed.46,47 To simplify ice growth in our initial tests of this slab-seeded approach, in this work, we only pick ice growth from the basal planes of the ice seed, as this allows uni-directional ice growth. Growth from other planes of the ice crystal can be tested in future work. Three different types of surfaces are used, i.e., FCC, BCC, and HCP, respectively. The lateral dimensions of the surfaces are equal to that of the ice slab and their thicknesses are determined by their unit cells, with two stacked layers used for all lattices. This results in surfaces comprising 704 (FCC), 308 (BCC), and 352 (HCP) atoms, respectively. Note that the planes of the surfaces exposed to the liquid are the FCC(001), BCC(001), and HCP(0001), respectively. The effect of varying the crystal plane facet was not explored in this study, and all results presented are for the facets listed above. As the crystal orientation is fixed for all surfaces, we use the shorthand FCC, BCC, and HCP in this paper when referring to the surfaces studied.

The freezing process is initiated by an ice seeding slab placed directly in contact with the supercooled water below it [Fig. 1(a)], and Figs. 1(b)1(d) illustrate how it proceeds on an FCC surface. This results in ice growth which approaches the surface over time [Figs. 1(b)1(d)]. The mean per-molecule potential energy of the water molecules during the freezing process is also tracked [Fig. 1(e)], and the drop in density of the water as it freezes results in a corresponding decrease in the mean potential energy of the system, until it reaches a final steady state value at roughly t = 80 ns. This corresponds to the completion of the freezing process as all available water molecules have been transformed into a hexagonal ice structure [see image at t = 120 ns in Fig. 1(d)]. Once steady state is reached, the simulation is run for an additional 80 ns, and all steady-state properties presented in this paper are obtained by averaging across this time period.

The wettability of the surface was varied by selecting four values of the surface–surface interaction parameter ɛss, which, in turn, changes the surface–oxygen interaction parameter ɛso and therefore the contact angle, as shown in Fig. 1(f). Consistent with previous MD studies, the calculation of the contact angle in MD simulations was designed to mimic experimental procedures. 2D axisymmetric slices of the droplet were taken, in which density was measured. This density was evaluated from the coordinates of all oxygen atoms, assuming axisymmetry around a centroidal axis normal to the solid surface. The whole domain was divided in cells, each of size 1 Å, and the local number density of every cell was calculated. The density of the cell varies as they go from bulk of the droplet to vapor, and isocontour lines were fitted to the density field. Then, the isocontour with 50% of the bulk liquid density was chosen, and using a circular fit, the tangent near the wall surface was drawn giving the contact angle. This process was then repeated for three droplet sizes, comprising 5000, 8000, and 11 000 molecules, respectively. Negligible variation in measured values of θ was observed for all droplets in all cases, confirming that the measurements reported here are independent of the droplet size.

The linear fit in Fig. 1(f) shows that the rate of decrease in θ is steeper for FCC when compared to BCC and HCP (which have similar slopes). This indicates that the wettability of the surfaces depends not only on the interaction strength but is also affected by the structure. In fact, for ɛso = 0.44 kcal/mol, the FCC surface becomes hydrophilic (θ ∼46°, i.e., much lower than 90°), while the BCC and HCP remain conventionally hydrophobic with θ > 90°. These differences likely arise from the fact that for the same structure, the FCC surface is denser with more atoms than the other two surfaces. As a result, the same value of ɛso results in a greater attractive force for the droplet placed on the FCC surface when compared to BCC/HCP surfaces, thus reducing the contact angle [Fig. 1(f)]. Note that σss = 2.573 Å was fixed for all cases, a value within the range typical of metals on which the surfaces in this study are based.

This resulted in a total of 12 different surfaces with surface properties listed in Table I. As ice growth is a stochastic process, four independent realizations are run for each surface, and all results presented in this work are obtained by averaging across the realizations. The equations of motion were integrated using the velocity Verlet algorithm, and the cut-off distances for all LJ and Coulombic interactions are 12 and 10 Å, respectively, and the skin distance for the neighboring lists is 2 Å. Long-range Coulombic interactions are computed using the PPPM method. System equilibration at 240 K was performed using the NVT ensemble via the Nosé–Hoover thermostat. After an equilibration time of 2 ns, the thermostats applied to the surface and the liquid region were turned off, while the thermostat applied to the seeding ice slab is kept for the remainder of the simulation. Given that the planar area of the ice interface is larger than ten molecular diameters in each direction, the simulations are not expected to have finite size effects on the interface properties.51 All simulations in this work were run using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) MD software.52 

TABLE I.

Table showing the different solid–solid (ss) and solid–water (so) Lennard-Jones interactions tested for all the three MD surfaces. σss (also σso) was fixed for all cases. Lorentz Berthelot mixing rules53 were used to calculate values of ɛso and σso.

ɛss (kcal/mol)ɛso (kcal/mol)σss (Å)σso (Å)
6.0 1.12 2.573 2.870 
0.91 0.44 2.573 2.870 
0.27 0.24 2.573 2.870 
0.017 0.06 2.573 2.870 
ɛss (kcal/mol)ɛso (kcal/mol)σss (Å)σso (Å)
6.0 1.12 2.573 2.870 
0.91 0.44 2.573 2.870 
0.27 0.24 2.573 2.870 
0.017 0.06 2.573 2.870 

1. 1D density plots

The major focus of our results will be on investigating the surface–ice interfacial region. In Fig. 2, we show the steady-state plots of 1D density for frozen water on all surfaces (FCC, BCC, and HCP) for all contact angles (going from superhydrophilic to superhydrophobic) considered in this work. Note that the black dashed lines, representing the density of the supercooled water prior to the initiation of the freezing process, are also provided for reference. The first few layers of water at the surface–ice interface are required to bind both the atoms forming the surface below and the bulk ice above. Inside bulk ice, water molecules typically are arranged into a hexagonal network, where molecules orient to optimize hydrogen bonding with layers above and below. This is commonly known as an ice bilayer when viewed from the basal face, as two peaks are observed in the density distribution [see Fig. 2(a-l), between ∼10 and ∼50 Å from the surface], each corresponding to the most probable location of the oxygen atoms (with the hydrogen atoms concentrated in the valley between the peaks).

FIG. 2.

Steady-state density profiles organized according to θ (rows) and surface topography (columns) after freezing for FCC [(a)–(d)], BCC [(e)–(h)], and HCP [(i)–(l)], respectively. The dashed black lines represent the density profiles of liquid water before crystallization and the solid lines represent the steady-state profiles. The numbers 1–6 present in each graph represent the layers closest to the surface (1–3) and vapor (4–6) interfaces, respectively, which are subject to investigation in this work. To compute density, the domain was divided along the z-axis into 1100 cells and the number of molecules in each cell was divided by the cell volume.

FIG. 2.

Steady-state density profiles organized according to θ (rows) and surface topography (columns) after freezing for FCC [(a)–(d)], BCC [(e)–(h)], and HCP [(i)–(l)], respectively. The dashed black lines represent the density profiles of liquid water before crystallization and the solid lines represent the steady-state profiles. The numbers 1–6 present in each graph represent the layers closest to the surface (1–3) and vapor (4–6) interfaces, respectively, which are subject to investigation in this work. To compute density, the domain was divided along the z-axis into 1100 cells and the number of molecules in each cell was divided by the cell volume.

Close modal

However, pristine ice bilayers showing well-developed double peaks are not observed in the interfacial region38 [see Fig. 2(a-l), between 0 and ∼10 Å from the surface]. Instead, depending on the type and temperature of the surface, water can form any number of one-dimensional chain networks (pentamers, hexamers, etc) or two-dimensional hydrogen bond networks (HBNs), commonly known as buckled structures. Typically, the molecular structure of the surface greatly influences the first layer of water adjacent to it, known as the first overlayer [see Fig. 1(g) for an illustration of a flat first layer and buckled bulk bilayers]. Between the first overlayer and bulk ice, there are usually one or more intermediate layers, which are closer in structure to the bulk ice bilayer but still distinct from it due to the influence of the surface and the first overlayer.54 

The complex balance between water–water and water–surface interactions produces diverse molecular arrangements in interfacial layers, which we will investigate in detail in this paper. For example, on surfaces to which water molecules bond very strongly, such as ruthenium (a reactive metal), experiments have observed 1D flat (or planar) hexamers in the first overlayer.55 These hexamers optimize the interaction between the water and the metal molecules, as all water molecules in the first overlayer are located as close to the surface as possible.55 In contrast, on metals for which the bond with water molecules is weaker, such as copper and silver, 2D buckled hexamers are favored. Here, buckled structures optimize the bonding within the water molecules.56 The transition between a flat (1D) to a buckled (2D) first overlayer results when water molecules favor bonding to each other over bonding to the surface.

To facilitate the discussion in this paper, the interfacial region at both the surface–ice and ice–vapor interfaces have been separated into three distinct layers each, which consist of the first overlayer and any intermediate layers. As can be seen in Fig. 2, they are labeled from 1–6, with 1 being the closest to the surface and 6 being closest to the vapor region. The importance of these layers arises because water molecules inside them interact with either vapor or the surface, which differs from bulk ice or bulk water. Additionally, in the absence of a surface, layers 1, 2, and 3 would be identical to layers 6, 5, and 4, respectively. Therefore, the differences between corresponding layers can also help shed light on the influence of the surface on the interfacial region.

Focusing on the FCC surface first [Figs. 2(a)2(d)], a few observations can be made as follows:

  1. In Figs. 2(a) and 2(b), for low values of the contact angle (θ ∼0° and θ ∼46°), i.e., for superhydrophilic and hydrophilic FCC surfaces, the peak in layer 1 is very sharp and appears at ∼2 Å from the surface, indicating a flat first overlayer where the water molecules bind strongly to the surface due to its high wettability (despite the surface not being reactive like Ruthenium). While the peak height varies in Figs. 2(a) and 2(b) as θ increases, in both cases, the unimodal peak corresponds to the first hydration layer. For increasing values of θ, the height of layer 1’s peak is strongly reduced. Thus, the density of the flat first overlayer decreases as the surface wettability decreases.

  2. In Fig. 2(c) (θ ∼108°) layer 1 peaks at ∼2 Å and a smaller shoulder appears at 3 Å, while in Fig. 2(d) (θ ∼161°), the inverse is observed, with a shoulder at ∼2 Å and layer 1 peak at ∼5 Å. This implies that the first overlayer has now transitioned from a flat to a buckled configuration, as hydrogen bonds between layers of water begin to be favored as wettability is reduced.

  3. The appearance of well-developed double peaks, representing the formation of an ice bilayer structure, occurs further away from the FCC surface, approximately at 10 Å in Figs. 2(a) and 2(b). These bilayers are observed at 12 Å from the surface in Figs. 2(c) and 2(d), implying that the ice bilayer starts further away from the surface as θ is increased.

  4. Comparing the density profiles of ice with that of supercooled water over the four FCC surfaces before freezing (dashed black line), it is observed that the density of water before and after ice growth are similar when very close to the surface and only start to differ at z > 4 Å. This indicates that the molecules closest to the surface are less influenced by the freezing process, and their structure is largely determined by the surface–ice interface regardless of whether they form part of supercooled water (before freezing) or ice (once freezing is completed).

Figure 2(e-l) shows the density profiles over the other two surfaces, BCC and HCP, for different wettabilities. There is a key difference when compared to FCC, which is that the first overlayer is not flat for any of the wettabilities studied, and some type of buckling is always observed. Layer 1 on the superhydrophilic BCC and HCP surfaces is consequently thicker than that of the FCC with a first peak at ∼1.2 Å and a second peak almost double in height at ∼3 Å and the onset of double peaks starting at ∼6 Å. In Figs. 2(f) and 2(j), (θ ∼106° and θ ∼110°, respectively) layer 1 has a large peak at ∼3 Å, and layer 2 comprises a degenerated double peak consisting of a shoulder at ∼6 Å and a peak at ∼7.5 Å.

Well-developed double peaks (indicating ice bilayers) are observed on BCC/HCP surfaces at ∼12 Å and beyond. Comparing the profiles beyond 12 Å, a small shift toward the right of the profile in Fig. 2(j) is noted, but for the others, BCC and HCP show a similar profile to FCC with equal distance between the bilayers (∼3.8 Å) and similar bilayer height. Comparing the density profiles over BCC and HCP surfaces, for increasing values of θ, layer 1 is observed to move further away from the surface. Moreover, when compared to the density profile of supercooled water, differences are observed on superhydrophobic surfaces at z ∼3 Å, which is lower when compared to superhydrophilic surfaces. This can be understood from the waning influence of the surface on the molecules near it as wettability is reduced, and the intermolecular interactions between the surface and water become weaker.

In summary, the decrease in ɛso results in a decrease in the first overlayer density (lower layer 1 peak), as the surface is unable to accommodate the same number of liquid molecules. Additionally, the density of the surface also plays a role. As the FCC surface has a higher wettability than BCC and HCP for the same value of ɛso, there is a stronger resultant attractive force toward the water molecules producing a flat first overlayer for FCC, which eventually transitions to a buckled overlayer at lower wettabilities. Note that the differences observed in the steady-state density profiles in Fig. 2 are limited to the vicinity of the surface–ice interface. The ice bilayers and the layers 4–6 (which constitute the ice–vapor interface) are identical on all surfaces. This is expected since the influence of the surface wanes considerably as you get further away from it.

2. 2D density heat maps

To probe further into the behavior of layers 1–2, 2D density “heat maps” were created for each layer, wettability, and structure. This allows us to analyze the local distribution of water molecules in each layer parallel to the wall, which would not be possible using just the 1D distribution. It has been noted in Sec. III A 1 that the water molecules on an FCC structure for θ = 0° forms a flat first overlayer. This flat structure can be understood as water maximizing its interaction with the surface with only a slight cost in hydrogen bond energy. In Fig. 3, the heat map of the first overlayer is placed alongside the heat map of the FCC structure itself. It is well-known that the surface can play a key role in templating the structure of water layers adjacent to it, sometimes forcing the first layer to mimic its registry rather than that of ice.57, Figure 3 shows distinct similarities between the FCC surface heat map and that of the first overlayer with many FCC atoms binding adjacent water molecules, which indicate some surface templating in this layer. However, unlike the cubic structure of the surface, hexamer and pentamer rings are observed (a couple of which are highlighted using dotted white lines in Fig. 3), consistent with experimental observations for water overlayers on FCC platinum surfaces.58 

FIG. 3.

(a) Two-dimensional heat map of FCC surface wall atoms at 240 K and (b) two-dimensional heat map of water molecules in layer 1 on top of the FCC surface for θ = 0°. For all heat maps, the density was computed by dividing the domain of each layer along the x-axis and y-axis into 5670 cells, and the mean number of molecules in each cell was divided by the cell volume. Note that the thickness, and therefore the volume, of individual layers may vary depending on the surface type and wettability. Note that the layer dimensions are measured trough-to-trough, which is to say that the valleys on either side of the peak(s) form the boundaries of a particular layer in our measurement.

FIG. 3.

(a) Two-dimensional heat map of FCC surface wall atoms at 240 K and (b) two-dimensional heat map of water molecules in layer 1 on top of the FCC surface for θ = 0°. For all heat maps, the density was computed by dividing the domain of each layer along the x-axis and y-axis into 5670 cells, and the mean number of molecules in each cell was divided by the cell volume. Note that the thickness, and therefore the volume, of individual layers may vary depending on the surface type and wettability. Note that the layer dimensions are measured trough-to-trough, which is to say that the valleys on either side of the peak(s) form the boundaries of a particular layer in our measurement.

Close modal

In Sec. III A 1, we observed significant changes in layers 1 and 2 as we moved to higher values of θ (Fig. 2). These changes are studied further here using these heat maps. To isolate the effect of the surface on layers 1–2, heat maps of layers 5–6 have also been produced for comparison, which also constitute an interface but are not influenced by the surface (as it is sufficiently far away from it for intermolecular interactions to be negligible). As was the case for density, heat maps of layers 1, 2, and 3 would, in the absence of the surface, be identical to those of layers 6, 5, and 4, respectively. It can be seen in Fig. 2(g-l) that the intermediate layer 2 appears visually similar to layer 5 (for θ ≥ 138°), as both consist of poorly formed bilayer-like structures, with asymmetrical peaks. This implies that the surface likely plays a much weaker role in water ordering in layer 2 once θ ≥ 138°.

We confirm this observation in Figs. 4 and 5, which shows heat maps for layers 2 and 5 on BCC (for θ = 138°) and on FCC (for θ = 161°), respectively. All four heat maps broadly show similar structures. This similarity is observed despite the different crystal structure of the underlying surface, which indicates that the registry of the underlying surface is no longer discernible. In addition, hexagonal patterns reminiscent of bulk ice start forming, with areas of extremely low density in between (colored in white). Note that these patterns are not exact because some defects in the water-layer structure are still present.

FIG. 4.

Two-dimensional heat maps for layers 2 and 5 on the BCC surface when θ = 138°. For this hydrophobic surface, layers 2 and 5 for both BCC and HCP surfaces also show similarities in their density plot (see Fig. 2).

FIG. 4.

Two-dimensional heat maps for layers 2 and 5 on the BCC surface when θ = 138°. For this hydrophobic surface, layers 2 and 5 for both BCC and HCP surfaces also show similarities in their density plot (see Fig. 2).

Close modal
FIG. 5.

Two-dimensional heat maps for layers 2 and 5 on the FCC surface when θ = 161°. For this superhydrophobic surface, layers 2 and 5 for all surfaces also show similarities in their density plot (see Fig. 2).

FIG. 5.

Two-dimensional heat maps for layers 2 and 5 on the FCC surface when θ = 161°. For this superhydrophobic surface, layers 2 and 5 for all surfaces also show similarities in their density plot (see Fig. 2).

Close modal

When the contact angles are high enough for our surfaces to be characterized as superhydrophobic (θ ≥ 150°), layer 1 also starts to resemble layer 6, as seen in Figs. 2(d)2(l), which is confirmed in Fig. 6. No structure or patterning is observed, indicating that these layers feel neither the registry of the surface below it nor the ice above. This verifies that for superhydrophobic surfaces, the surface–ice interface is nearly identical to the ice–vapor interface.

FIG. 6.

Two-dimensional heat maps for layers 1 and 6 on the HCP surface when θ = 180°. For this and other superhydrophobic surfaces, layers 1 and 6 for all surfaces also show similarities in their density plots (see Fig. 2).

FIG. 6.

Two-dimensional heat maps for layers 1 and 6 on the HCP surface when θ = 180°. For this and other superhydrophobic surfaces, layers 1 and 6 for all surfaces also show similarities in their density plots (see Fig. 2).

Close modal

1. Distinguishing ice from non-ice molecules

While the density heat maps show the local concentration of molecules within layers, they do not directly reveal if they are liquid, ice, or neither. This information is obtained using the CHILL+ algorithm,59–61 which identifies ice and non-ice molecules using the HBN between them. Note that here the term “non-ice molecules” and not “liquid molecules” is used because not all non-ice molecules are necessarily in the liquid phase; CHILL+ also identifies any poorly formed bilayers or structures with pentamer/hexagonal rings that do not resemble ice as comprising non-ice molecules. Figure 7 illustrates how CHILL+ is applied to a typical case (FCC, θ = 46°). Figure 7(a) is the same density plot as before but with ice molecules separated from non-ice molecules using CHILL+, while Fig. 7(b) shows the distribution of ice and non-ice molecules in the domain. The minimal presence of ice molecules in the interfacial layers (between 0–10 Å) is expected. A similar zone, but reduced in size, is present in the ice–vapor interface in this instance (between 48–55 Å). CHILL+ also identifies a few non-ice molecules in the frozen bulk region (between 20–40 Å). This indicates the presence of defects in the formed ice, as the density plot shows that these water molecules are arranged in the form of imperfectly developed (slightly asymmetrical) bilayers.

FIG. 7.

(a) The 1D density plot of non-ice (green) and ice (red) water molecules as a function of distance from the wall using the CHILL+ algorithm and (b) an MD snapshot of their distribution in the system. CHILL+ relies on the correlation of bond-order parameters to identify and count the number of staggered and eclipsed O–O bonds in the HBN, which distinguish hexagonal/cubic ice and non-ice molecules.

FIG. 7.

(a) The 1D density plot of non-ice (green) and ice (red) water molecules as a function of distance from the wall using the CHILL+ algorithm and (b) an MD snapshot of their distribution in the system. CHILL+ relies on the correlation of bond-order parameters to identify and count the number of staggered and eclipsed O–O bonds in the HBN, which distinguish hexagonal/cubic ice and non-ice molecules.

Close modal

Figure 8 plots the average number of non-ice molecules in layers 1–3, respectively, for all structures and wettabilities. Despite the significant difference highlighted previously in the first overlayer structure between FCC and BCC/HCP for θ = 0°, Fig. 8 shows that the number of non-ice molecules in layer 1 is similar (∼180) in all three surfaces, regardless of the flat/buckled configuration. Interestingly, it can be seen that the greatest variation in the number of non-ice molecules (with increasing wettability) is seen in layer 2, as the values can range between ∼90 (all surfaces, θ ∼180°) to ∼180 (FCC, θ ∼46°). Additionally, the number of non-ice molecules in layer 2 seems to vary in a non-monotonic manner with wettability. Layer 3 behaves similarly but the variation is smaller, as it consistently contains less than ∼80 non-ice molecules for all surfaces and wettabilities, with the number dropping to ∼20 once θ ≥ 130°. To explain these differences, we need to examine the dynamics of how water molecules behave inside these layers in more detail, which is the goal of Sec. III B 2.

FIG. 8.

Number of non-ice molecules in all three surfaces for all wettabilities, with indices inset to indicate the layer number as defined in Fig. 2.

FIG. 8.

Number of non-ice molecules in all three surfaces for all wettabilities, with indices inset to indicate the layer number as defined in Fig. 2.

Close modal

In Fig. 9, the water densities in Fig. 2 are replotted, but with ice and non-ice molecules separated using CHILL+ as shown previously in Fig. 7 (with non-ice density shown as a black line). It is apparent that layers 1–3 near the surface–ice interface and layers 5 and 6 at the ice–vapor interface are the only layers where the majority of the molecules are identified by CHILL+ to be non-ice. Note that for the superhydrophobic cases [Figs. 9(d)9(l)], layer 3 freezes as well, thus resembling the bulk ice layers. In fact, as θ increases, the surface–ice interface starts to resemble the ice–vapor interface, consistent with Sec. III B. Once we reach the superhydrophobic cases, we see that layer 1 resembles layer 6 in all three cases (recall Fig. 6), showing the expected absence of attractive forces from the surfaces to the liquid molecules. This is why distinguishing between ice and non-ice molecules is important, as comparing the surface–ice and ice–vapor interfaces now gives insight into the influence of the surface on the composition of water molecules (whether ice or non-ice) in the interfacial region, since ice/non-ice molecules can behave differently in their mobility and ability to form structures, which we will demonstrate in Secs. III B 3 and III B 4. For example, the heat maps for layers 2 and 5 of the BCC and FCC (for θ > 138°) are seen to be similar in Figs. 4 and 5, respectively. Here, it can be seen that the peaks forming the bilayers in layers 2 and 5 can be subdivided into a peak that is identified to be non-ice (close to either interface) and the peak identified as ice (closer to the bulk ice). This highlights both why the heat maps showed hexagonal structures and had defects, as these layers were only partly frozen, which means that they are partly liquid-like in structure. Such quasi-liquid layers (QLLs) have, in fact, been experimentally observed at the surface–ice interface.62–66 As mentioned earlier, the formation of QLLs is associated with low ice adhesion,39,40 which is desirable when designing novel icephobic surfaces. However, previous work using MD to study QLLs has focused primarily at the ice–vapor interface,67–68 and this paper represents the first attempt to study QLLs that form at the surface–ice interface, which is conducted next.

FIG. 9.

Steady-state density profiles similar to Fig. 2, but here the profiles are split into ice (color solid lines) and non-ice (black solid lines) molecules using the CHILL+ algorithm. See Fig. 2 for definition of indices 1–6.

FIG. 9.

Steady-state density profiles similar to Fig. 2, but here the profiles are split into ice (color solid lines) and non-ice (black solid lines) molecules using the CHILL+ algorithm. See Fig. 2 for definition of indices 1–6.

Close modal

2. Quantifying quasi-liquid layers

From an MD perspective, as long as the molecules in the interfacial region are not frozen (i.e., they are identified by CHILL+ to be non-ice), it can be included in our QLL thickness measurements using the following equation:67 
TQLL(Å)=Nnon-iceMρNAVLxLy1024,
(2)
where Nnon-ice is the average number of non-ice molecules, M is the water molecule’s mass in g/mol, NAV is Avogadro’s number, ρ is the bulk density of liquid water (g/cm3), and LxLy is the area of the interface in Å2. In cases where any molecules in any one of the three interfacial layers forms ice, it is excluded from our QLL measurements. This protocol enables the calculation of QLL thickness, providing a quantitative estimate of non-frozen molecules, and enables comparisons between the different simulated cases. Note that while we are not explicitly calculating a “thickness” in this instance (we are estimating an equivalent thickness based on the number of non-ice molecules that are present in the interfacial region), we did also quantify TQLL graphically (so by directly measuring thickness from the density plots). However, as both approaches produce reasonably similar results, we have only included this analytical calculation of TQLL in this paper.

Figure 10 plots TQLL at the surface–ice interface (markers) and ice–vapor interface (blue dashed line) for all surfaces and wettabilities. It shows that TQLL is higher for moderate values of θ (between 40° and 120°) and lower for both superhydrophobic and superhydrophilic surfaces. This is an interesting result that deserves further scrutiny. For low values of θ, the surface overlayer (be it flat or buckled) is strongly bonded to the surface molecules and therefore is identified by CHILL+ to be non-ice, which then means that it contributes to TQLL. Therefore, a thicker TQLL for low θ can be understood as resulting from the fact that the surfaces force molecules in layers 1–3 to adopt non-hexagonal (i.e., non-ice) configurations. What is less apparent is why TQLL peaks for moderate values of θ and then drops off at higher θ, in a similar pattern seen in Fig. 8 for layers 2 and 3. This will be investigated in Sec. III B 3. Second, TQLL is shown to approach the value of the ice–vapor interface at high values of θ. This is expected as the surface plays a minimal role at high θ, and this convergence can be considered to be a validation of these measurement techniques. Note that the TQLL values approach but do not actually achieve the ice–vapor interface values.

FIG. 10.

The thickness of the quasi-liquid layer, TQLL, is plotted using Eq. (2). The blue dashed line shows the corresponding value of TQLL at the ice–vapor interface.

FIG. 10.

The thickness of the quasi-liquid layer, TQLL, is plotted using Eq. (2). The blue dashed line shows the corresponding value of TQLL at the ice–vapor interface.

Close modal

3. Molecular mobility in every layer

There have been two key observations made thus far: (a) a flat-to-buckled transition in the first overlayer occurs on the FCC surface as θ is increased from 46° to 108° (Fig. 2), while the BCC and HCP surfaces show buckled overlayers throughout; and (b) the thickness of the interfacial region (Fig. 10) is largest for moderate values of θ (between 40° and 120°) and is lower for superhydrophobic (θ > 160°) and superhydrophilic (θ = 0°) surfaces. To investigate further what is going on in the interfacial region across different surface wettabilities, we now turn to the mobility of the molecules in each layer. This is achieved by estimating the self-diffusion coefficient D for each layer in steady state using the formula
D=12nx(t)x02t,
(3)
where n = 3 is the diffusion space dimensionality and
x(t)x02=1Ni=1Nxi(t)x0i2
(4)
is the mean squared displacement (MSD) measurement, where N is the number of molecules and x0i is the initial position of the ith molecule. The MSD was calculated separately for layers 1–3 by averaging over the molecules that belong to each layer. Note that molecules can freely travel between the layers during this measurement, but the entrance and exit times for each molecule as it moves between layers are tracked. This means that if a molecule leaves the layer, it stops being incorporated into the calculation at that exact time step. Instead, its velocity contributes to the MSD of the layer it has entered. If this molecule then re-enters the original layer at a future time step, its entry point is treated as the starting position for MSD calculations for that time step. Thus, all the molecules within a given layer are included in the MSD calculation.

For example, there are about 150–200 molecules in the first layer on all surfaces (see Fig. 8), the bulk of which contribute to the MSD calculations, which are performed for every nanosecond of the averaging time (the averaging time consists of 80 ns in total). Each ns-long calculation proceeds as follows: given a simulation time step of 2 fs, the calculations take place every 10 000 time steps (so every 20 ps), which gives a total of 50 values in 1 ns. Once the MSD is calculated separately for every 1 ns, the 80 values obtained are then averaged at the end of the averaging time. As the MSD is approximately linear over time when the system is in steady state, then D can be calculated from the slope of the MSD plot. Validation of the method described above was performed by calculating D for bulk supercooled water at 240 K, and the value obtained (D = 0.199 × 10−9 m2/s) is close to the experimental value at 239.8 K69 (Dexp = 0.204 × 10−9 m2/s).

A high self-diffusion coefficient D implies that there is significant molecular motion, whereas a very low value of D implies that molecules are locked into a crystalline structure (for example, in a solid). Figure 11 shows the values of D for all surfaces and wettabilities. Here, it is important to note that (a) for the purposes of this calculation, no distinction was made between ice and non-ice molecules and (b) ice molecules have very low D, given they are in the solid phase. This may mean that as the relative proportion of ice to non-ice molecules within a layer changes with θ, the value of D will depend on the largest occupant (ice/non-ice molecules), as it is a molecule-averaged measure of mobility. To avoid confusion, the values of D for layers 5 and 6 (at the ice–vapor interface) are also included for reference. Note that as layers 5 and 6 are unaffected by the surface, their respective D values are constant (D = 8.486 × 10−12 m2/s and D = 1.355 × 10−10 m2/s, respectively), irrespective of surface wettability. At the surface–ice interface, D for layer 3 is similarly almost unaffected by θ or surface type and fluctuates around D = 1.277 × 10−13 m2/s. On the other hand, for layers 1–2 (first overlayer and intermediate layer), the surface influences the self-diffusion characteristics significantly, and a strong dependence on θ is seen, which is explored below.

FIG. 11.

The self-diffusion coefficient D plotted for all values of θ for layers 1–3 on (a) FCC, (b) BCC, and (c) HCP surfaces. The D values of layers 5 and 6 are provided as dashed lines for comparison.

FIG. 11.

The self-diffusion coefficient D plotted for all values of θ for layers 1–3 on (a) FCC, (b) BCC, and (c) HCP surfaces. The D values of layers 5 and 6 are provided as dashed lines for comparison.

Close modal

4. Untangling the effects of surface templating

Focusing on layer 1 first (plotted as black lines in Fig. 11), the variation of D as surface wettability θ is varied is qualitatively similar for all three surfaces; D is negligible at low θ and rises with θ with its highest value of D for the superhydrophobic surfaces considered here. Essentially, the molecules in the first overlayer are solid-like at low θ (D ∼0) as they are strongly adsorbed onto the surface atoms. However, as the influence of the surface atoms wanes (as θ increases), these molecules are freed from their binding sites, gain mobility, and move more freely, resulting in a rising D. Thus, the higher the wettability of a surface, the lower the mobility the molecules will have in the first overlayer. However, there are quantitative differences between the FCC surface and the BCC/HCP surfaces. For the FCC surface, the slope of D with θ appears almost linear, whereas both the BCC and HCP surfaces show a power-law behavior, with a steep rise in D as θ is increased from 100° to 180°. For superhydrophobic cases, the surface–ice interface resembles the ice–vapor interface as previously mentioned, and this is reflected in the values of D in layer 1 as well.

Moving on to the intermediate layer 2 (plotted as red markers and lines in Fig. 11), the variation of D with θ is more complex and further resembles trends seen earlier in Fig. 10. D is almost negligible for both extremely low (0°) and extremely high (160° and above) values of θ. In between those extremes, a non-monotonic behavior is observed, analogous to trends from Figs. 8 and 10. A rise in D is seen as θ is increased from 0° to 100° and then a subsequent drop in D is seen beyond that when θ rises from 100° to 180°. Conceptually, this can be understood as resulting because of two competing mechanisms: (a) the surface and the first overlayer (collectively called the surface templating layers here) attempting to structure layer 2 from below and (b) the ice-like layer 3 and the bulk ice layers above it (collectively called the freezing layers here) attempting to freeze it from above. The value of θ determines the relative strength of these two effects. At very low values of θ, the layers are dominated by surface templating, resulting in a highly structured solid-like layer 2. Similarly, at very high values of θ, the surface has very little influence, and it is instead the freezing layer 3 and the bulk ice above that structures layer 2 into an ice-like structure. In either of these extreme cases, when the freezing layers or the templating layers dominate, the mobility of the water molecules is greatly hindered. This reduces the number of possible configurations the water molecules in layer 2 can explore, as they are constrained to specific binding sites. These sites correspond to either the location of surface/overlayer atoms (low θ) or a structured ice-like HBN (high θ). Hence, D → 0 when θ → 0° or θ → 180°.

However, for moderate values of θ, both the templating and the freezing layers are of similar importance, and as a consequence, they balance each other; layer 2 molecules are therefore not strongly manipulated by either effect. This means that the motion of these molecules is not hindered as significantly, and they are able to continuously move in and out of template-like and ice-like configurations. This is reflected in their higher mobility and therefore a higher D, which appears to peak at θ = 100°. This therefore explains one of the two key observations, which is why TQLL peaks at moderate values of θ. It is now clear that this change in TQLL is driven primarily by changes in layer 2. Note that this competition between the templating and the freezing layers also explains the non-monotonic variation in non-ice molecule count with wettability observed in layer 2 (Fig. 8). The non-ice molecule count peaks at moderate values of θ and drops significantly as θ → 0° or θ → 180°, when the templating and freezing layers gain importance, respectively. In contrast, the non-ice molecule count in layers 1 and 3 in Fig. 8 does not vary as much with θ as there is no surface/ice-mediated competition occurring in these layers as θ → 0° or θ → 180°.

5. Understanding layer 2 on the FCC surface

While the above explanation is sufficient to explain the qualitative nature of the three MSD plots, it still does not touch on the quantitative differences between the surfaces, as D is far higher in layer 2 for FCC when compared to the other surfaces. This can be understood by the structural uniqueness of the FCC surface, which produces a flat overlayer at low θ and a transition to a buckled overlayer at higher θ. It is crucial to appreciate that the dynamic behavior of water molecules in a flat overlayer is distinct from that in a buckled overlayer. In a flat 1D overlayer, the water molecules are restricted to motion in a particular plane, and the space that the molecules can explore is effectively reduced by one dimension relative to the layers above. This means that they bind more effectively to any adsorption sites available on the surface, which impacts the structure of the water molecules in the layer above it. Specifically, the molecules in layer 2 experience a much smoother energy landscape, and the locations of the surface atoms are more smeared out when compared to the BCC/HCP surfaces.

This is evident in Fig. 12 where density heat maps of layer 2 of the interfacial region on the FCC surface are compared to those on the BCC and HCP surfaces. To remove any effect of wettability and ensure that the differences observed result from the surface structure, cases with very similar wettabilities are selected for all surfaces (θ between 106–110°). It is clear that the FCC surface produces far fewer binding sites for the water molecules, as both the BCC and HCP surfaces show some degree of water ordering in layer 2, whereas the FCC structure has minimal ordering. This means that the mobility for molecules in layer 2 of FCC is less hindered by the corrugations in the overlayer below, allowing for greater mobility (and thus a greater D). An important consequence of this result is that FCC surfaces are not representative of all crystalline structures; indeed, they behave very differently from the BCC and HCP surfaces studied in this paper. In particular, they demonstrate a flat-to-buckled transition in the first overlayer with θ and a comparative lack of ordering in the second layer allowing for greater D for middle values of θ. This is not observed for BCC and HCP surfaces.

FIG. 12.

Density heat maps for layer 2 on (a) the FCC surface with θ = 108°, (b) the BCC surface with θ = 106°, and (c) the HCP surface with θ = 110°.

FIG. 12.

Density heat maps for layer 2 on (a) the FCC surface with θ = 108°, (b) the BCC surface with θ = 106°, and (c) the HCP surface with θ = 110°.

Close modal

In this paper, ice formation on a selection of crystalline structures was achieved using a slab-seeding approach, which involved sandwiching supercooled water between a seeding ice slab and the surface. This enabled ice growth to be successfully observed on FCC(001), BCC(001), and HCP(0001) surfaces and investigate the physics of the surface–ice interfacial region in detail. Note that there are simplifications made in the proposed approach to make the simulations computationally tractable, such as the use of very simple surfaces and an operating temperature (240 K) lower than that in the case of real-world icing scenarios. These imposed restrictions can be relaxed in future work on this topic.

The effect of surface wettability was tested by varying the surface–water interatomic potentials ɛso. It was shown that the relationship between ɛso and the contact angle θ was surface-dependent; in particular, the FCC surface consistently showed lower contact angles for the same value of ɛso. This was attributed to the greater density of the FCC surface, which meant that the number of solid–liquid interactions are larger.

Looking at the densities of the water post-freezing on all surfaces, it became apparent that the FCC surface also differed in terms of the structure of the first overlayer. The first overlayer was always buckled for the BCC and HCP surfaces; in contrast, the FCC surface induced a flat first overlayer for high wettabilities (θ < 50°), and this transitioned to a buckled overlayer as the wettability was decreased (θ > 100°). We also found that the number of molecules contained in the first overlayer was similar for all three surfaces, but the FCC is able to suppress water–water hydrogen bond formation and accommodate them in one single flat layer at high wettabilities, whereas the other surfaces were unable to do so and formed buckled overlayers instead.

The molecules were separated into ice and non-ice using CHILL+, and the thickness of the interfacial region of non-ice molecules was quantified (TQLL). TQLL peaked for moderate wettabilities and decreased for both extremely high and extremely low wettabilities. The self-diffusion coefficient D of layers within the interfacial region showed a similar trend. The first overlayer and the surface compete with the third layer and the bulk ice to simultaneously apply their influence on the structure in layer 2, which was found to be important when wettability was modified. For low wettabilities, the influence of layer 3 dominated and layer 2 displayed ice-like structures, whereas for high wettabilities, the surface constrained the water molecules in layer 2, and they adopted positions consistent with the surface adsorption sites. In either case, D of the layer 2 molecules was lowered compared to when the freezing and surface templating layers were balanced at moderate values of θ. TQLL also peaked in the same range of moderate θ as the molecules were more likely to be identified to be non-ice when they were not constrained by the surface or the bulk ice.

Finally, while the qualitative behavior of D with θ was similar across surfaces, there were significant quantitative differences. Importantly, the values of D were higher in layer 2 on the FCC surface at moderate wettabilities. This was shown to result from the screening effect of the flat overlayer below, which prevented layer 2 water molecules from being bound to the surface adsorption sites. Consequently, this meant that the layer 2 molecules on the FCC surface were able to move more freely compared to their counterparts on the BCC and the HCP surfaces. This is relevant to the MD studies of heterogeneous nucleation/growth in the literature, which have been predominantly conducted on FCC surfaces. Our results show that FCC surfaces induce qualitative and quantitative differences in the interfacial region when compared to other common crystalline structures, which may be relevant depending on the application.

See the supplementary material for additional simulation results that were run to verify the validity of the setup used in this work. Additional data (including a sample LAMMPS input script and datafile) can be freely accessed at https://doi.org/10.7488/ds/3221.

This work was supported in the UK by the Engineering and Physical Sciences Research Council (EPSRC) under Grant Nos. EP/N016602/1 and EP/R007438/1. All MD simulations were run on ARCHER, the UK’s national supercomputing service.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
O.
Parent
and
A.
Ilinca
,
Cold Reg. Sci. Technol.
65
,
88
(
2011
).
2.
R. W.
Gent
,
N. P.
Dart
, and
J. T.
Cansdale
,
Philos. Trans. R. Soc. London, Ser. A
358
,
2873
(
2000
).
3.
A. R.
Solangi
, “
Icing effects on power lines and anti-icing and de-icing methods
,” M.S. thesis,
The Arctic University of Norway
,
2018
.
4.
J.
Tian
,
X.
Dong
,
B.
Xi
,
P.
Minnis
,
W. L.
Smith
, Jr.
,
S.
Sun‐Mack
,
M.
Thieman
, and
J.
Wang
,
J. Geophys. Res.: Atmos.
123
,
1708
, (
2018
).
5.
H.
Huang
,
M. L.
Yarmush
, and
O. B.
Usta
,
Nat. Commun.
9
,
3201
(
2018
).
6.
T.
Darmanin
and
F.
Guittard
,
Mater. Today
18
,
273
(
2015
).
7.
X.
Huang
,
N.
Tepylo
,
V.
Pommier-Budinger
,
M.
Budinger
,
E.
Bonaccurso
,
P.
Villedieu
, and
L.
Bennani
,
Prog. Aerosp. Sci.
105
,
74
(
2019
).
8.
J.
Jiang
,
G. X.
Li
,
Q.
Sheng
, and
G. H.
Tang
,
Appl. Surf. Sci.
510
,
145520
(
2020
).
9.
P.
Eberle
,
M. K.
Tiwari
,
T.
Maitra
, and
D.
Poulikakos
,
Nanoscale
6
,
4874
(
2014
).
10.
Q. T.
Fu
,
E. J.
Liu
,
P.
Wilson
, and
Z.
Chen
,
Phys. Chem. Chem. Phys.
17
,
21492
(
2015
).
11.
S. A.
Zielke
,
A. K.
Bertram
, and
G. N.
Patey
,
J. Phys. Chem. B
120
,
1726
(
2016
).
12.
C.
Li
,
X.
Gao
, and
Z.
Li
,
J. Phys. Chem. C
121
,
11552
(
2017
).
13.
B.
Glatz
and
S.
Sarupria
,
J. Chem. Phys.
145
,
211924
(
2016
).
14.
B.
Glatz
and
S.
Sarupria
,
Langmuir
34
,
1190
(
2018
).
15.
V.
Molinero
and
E. B.
Moore
,
J. Phys. Chem. B
113
,
4008
(
2009
).
16.
J. L. F.
Abascal
,
E.
Sanz
,
R.
García Fernández
, and
C.
Vega
,
J. Chem. Phys.
122
,
234511
(
2005
).
17.
M.
Fitzner
,
G. C.
Sosso
,
S. J.
Cox
, and
A.
Michaelides
,
J. Am. Chem. Soc.
137
,
13658
(
2015
).
18.
L.
Lupi
,
A.
Hudait
, and
V.
Molinero
,
J. Am. Chem. Soc.
136
,
3156
(
2014
).
19.
Y.
Bi
,
R.
Cabriolu
, and
T.
Li
,
J. Phys. Chem. C
120
,
1507
(
2016
).
20.
J. K.
Singh
and
F.
Müller-Plathe
,
Appl. Phys. Lett.
104
,
021603
(
2014
).
21.
L.
Lupi
and
V.
Molinero
,
J. Phys. Chem. A
118
,
7330
(
2014
).
22.
S. J.
Cox
,
S. M.
Kathmann
,
B.
Slater
, and
A.
Michaelides
,
J. Chem. Phys.
142
,
184705
(
2015
).
23.
M. H.
Factorovich
,
V.
Molinero
, and
D. A.
Scherlis
,
J. Am. Chem. Soc.
137
,
10618
(
2015
).
24.
G. C.
Sosso
,
T.
Li
,
D.
Donadio
,
G. A.
Tribello
, and
A.
Michaelides
,
J. Phys. Chem. Lett.
7
,
2350
(
2016
).
25.
A.
Laio
and
M.
Parrinello
,
Proc. Natl. Acad. Sci. U. S. A.
99
,
12562
(
2002
).
26.
A.
Laio
and
F. L.
Gervasio
,
Rep. Prog. Phys.
71
,
126601
(
2008
).
27.
C.
Abrams
and
G.
Bussi
,
Entropy
16
,
163
(
2014
).
28.
C.
Dellago
,
P. G.
Bolhuis
, and
D.
Chandler
,
J. Chem. Phys.
108
,
9236
(
1998
).
29.
C.
Dellago
and
P. G.
Bolhuis
,
Adv. Polym. Sci.
221
,
167
(
2009
).
30.
R. J.
Allen
,
C.
Valeriani
, and
P.
Rein Ten Wolde
,
J. Phys.: Condens. Matter
21
,
463102
(
2009
).
31.
G. M.
Torrie
and
J. P.
Valleau
,
J. Comput. Phys.
23
,
187
(
1977
).
32.
R. J.
Allen
,
D.
Frenkel
, and
P. R.
Ten Wolde
,
J. Chem. Phys.
124
,
024102
(
2006
).
33.
J. R.
Espinosa
,
C.
Navarro
,
E.
Sanz
,
C.
Valeriani
, and
C.
Vega
,
J. Chem. Phys.
145
,
211922
(
2016
).
34.
P.
Pedevilla
,
M.
Fitzner
,
G. C.
Sosso
, and
A.
Michaelides
,
J. Chem. Phys.
149
,
072327
(
2018
).
35.
A.
Reinhardt
,
J. P. K.
Doye
,
E. G.
Noya
, and
C.
Vega
,
J. Chem. Phys.
137
,
194504
(
2012
).
36.
R.
Radhakrishnan
and
B. L.
Trout
,
Phys. Rev. Lett.
90
,
158301
(
2003
).
37.
D.
Quigley
and
P. M.
Rodger
,
J. Chem. Phys.
128
,
154518
(
2008
).
38.
O.
Björneholm
,
M. H.
Hansen
,
A.
Hodgson
,
L.-M.
Liu
,
D. T.
Limmer
,
A.
Michaelides
,
P.
Pedevilla
,
J.
Rossmeisl
,
H.
Shen
,
G.
Tocci
,
E.
Tyrode
,
M.-M.
Walz
,
J.
Werner
, and
H.
Bluhm
,
Chem. Rev.
116
,
7698
(
2016
).
39.
J.
Chen
,
R.
Dou
,
D.
Cui
,
Q.
Zhang
,
Y.
Zhang
,
F.
Xu
,
X.
Zhou
,
J.
Wang
,
Y.
Song
, and
L.
Jiang
,
ACS Appl. Mater. Interfaces
5
,
4026
(
2013
).
40.
J.
Chen
,
K.
Li
,
S.
Wu
,
J.
Liu
,
K.
Liu
, and
Q.
Fan
,
ACS Omega
2
,
2047
(
2017
).
41.
A.
Amirfazli
and
C.
Antonini
, in
Non-Wettable Surfaces: Theory, Preparation and Applications
(
The Royal Society of Chemistry
,
2017
), Chap. 11, pp.
319
346
.
42.
K.
Jha
,
E.
Anim-Danso
,
S.
Bekele
,
G.
Eason
, and
M.
Tsige
,
Coatings
6
,
3
(
2016
).
43.
D.
Chen
,
M. D.
Gelenter
,
M.
Hong
,
R. E.
Cohen
, and
G. H.
McKinley
,
ACS Appl. Mater. Interfaces
9
,
4202
(
2017
).
44.
Q.
Zeng
and
K.
Li
,
Crystals
9
,
250
(
2019
).
45.
R. G.
Pereyra
,
I.
Szleifer
, and
M. A.
Carignano
,
J. Chem. Phys.
135
,
034508
(
2011
).
46.
H.
Nada
and
Y.
Furukawa
,
J. Cryst. Growth
169
,
587
(
1996
).
47.
H.
Nada
and
Y.
Furukawa
,
Appl. Surf. Sci.
121–122
,
445
(
1997
).
48.
M. A.
Carignano
,
P. B.
Shepson
, and
I.
Szleifer
,
Mol. Phys.
103
,
2957
(
2005
).
49.
S.
Uchida
,
K.
Fujiwara
, and
M.
Shibahara
,
Nanoscale Microscale Thermophys. Eng.
24
,
53
(
2020
).
50.
J.-P.
Ryckaert
,
G.
Ciccotti
, and
H. J. C.
Berendsen
,
J. Comput. Phys.
23
,
327
(
1977
).
51.
M. M.
Conde
,
C.
Vega
, and
A.
Patrykiejew
,
J. Chem. Phys.
129
,
014702
(
2008
).
52.
S.
Plimpton
,
J. Comput. Phys.
117
,
1
(
1995
).
53.
H. A.
Lorentz
,
Ann. Phys.
248
,
127
(
1881
).
54.
J.
Carrasco
,
A.
Hodgson
, and
A.
Michaelides
,
Nat. Mater.
11
,
667
(
2012
).
55.
S.
Haq
,
C.
Clay
,
G. R.
Darling
,
G.
Zimbitas
, and
A.
Hodgson
,
Phys. Rev. B
73
,
115414
(
2006
).
56.
A.
Michaelides
and
K.
Morgenstern
,
Nat. Mater.
6
,
597
(
2007
).
57.
N.
Gerrard
,
C.
Gattinoni
,
F.
McBride
,
A.
Michaelides
, and
A.
Hodgson
,
J. Am. Chem. Soc.
141
,
8599
(
2019
).
58.
S.
Nie
,
P. J.
Feibelman
,
N. C.
Bartelt
, and
K.
Thürmer
,
Phys. Rev. Lett.
105
,
26102
(
2010
).
59.
E. B.
Moore
,
E.
de la Llave
,
K.
Welke
,
D. A.
Scherlis
, and
V.
Molinero
,
Phys. Chem. Chem. Phys.
12
,
4124
(
2010
).
60.
T.
Sayer
and
S. J.
Cox
,
Phys. Chem. Chem. Phys.
21
,
14546
(
2019
).
61.
A. H.
Nguyen
and
V.
Molinero
,
J. Phys. Chem. B
119
,
9369
(
2015
).
62.
Y.
Furukawa
, in
Handbook of Crystal Growth
, 2nd ed. (
Elsevier
,
Boston
,
2015
), Chap. 25, pp.
1061
1112
.
63.
K.-i.
Murata
,
H.
Asakawa
,
K.
Nagashima
,
Y.
Furukawa
, and
G.
Sazaki
,
Proc. Natl. Acad. Sci. U. S. A.
113
,
E6741
(
2016
).
64.
J. G.
Dash
,
A. W.
Rempel
, and
J. S.
Wettlaufer
,
Rev. Mod. Phys.
78
,
695
(
2006
).
65.
D. T.
Limmer
,
Proc. Natl. Acad. Sci. U. S. A.
113
,
12347
(
2016
).
66.
K.
Mochizuki
and
V.
Molinero
,
J. Am. Chem. Soc.
140
,
4803
(
2018
).
67.
J.
Benet
,
P.
Llombart
,
E.
Sanz
, and
L. G.
MacDowell
,
Phys. Rev. Lett.
117
,
096101
(
2016
); arXiv:1609.01439.
68.
P.
Llombart
,
P.
Llombart
,
E. G.
Noya
, and
L. G.
MacDowell
,
Sci. Adv.
6
,
eaay9322
(
2020
).
69.
W. S.
Price
,
H.
Ide
, and
Y.
Arata
,
J. Phys. Chem. A
103
,
448
(
1999
).
Published open access through an agreement with The University of Edinburgh

Supplementary Material