Surface phase transitions and crystal habits of ice in the atmosphere

See allHide authors and affiliations

Science Advances  20 May 2020:
Vol. 6, no. 21, eaay9322
DOI: 10.1126/sciadv.aay9322


With climate modeling predicting a raise of at least 2°C by year 2100, the fate of ice has become a serious concern, but we still do not understand how ice grows (or melts). In the atmosphere, crystal growth rates of basal and prism facets exhibit an enigmatic temperature dependence and crossover up to three times in a range between 0° and −40°. Here, we use large-scale computer simulations to characterize the ice surface and identify a sequence of previously unidentified phase transitions on the main facets of ice crystallites. Unexpectedly, we find that as temperature is increased, the crystal surface transforms from a disordered phase with proliferation of steps to a smooth phase with small step density. This causes the anomalous increase of step free energies and provides the long sought explanation for the enigmatic crossover of snow crystal growth rates found in the atmosphere.


The Nakaya diagram documents the hidden mystery of snow crystal growth (1): As temperature is cooled down from 0°C to −40°C, ice crystals in the atmosphere change their habit from plates, to columns, to plates and back to columns in a puzzling and unexplained sequence that holds the key to our understanding of the ice surface (Fig. 1).

Fig. 1 Growth of ice crystals at low supersaturation.

(Middle) As temperature ΔT = TTt decreases below the triple point, Tt, the habit of hexagonal ice prisms grown in the atmosphere changes sharply from platelike to columnar at ca. −4 K, from columnar to platelike at −10 K, and somewhat less sharply from platelike to columnar at −40 K (8). The facet which grows faster as indicated by the arrows, dictates the prevalence of plates or columns. The change of crystal habits results from a crossover in the growth rates of the basal and prism facets. (Top and bottom) Sketch of surface structural evolution with temperature as found in our work. Blue lines represent the i/f surface, and red lines represent the f/v surface. The basal surface (top row) is a high-temperature reconstructed flat (HT-RF) phase from ΔT = 0 to −4 K. It becomes a DOF phase in the range between ca. −4 and −47 K and is transformed into an OF phase at lower temperatures. In this phase, surface disorder resulting from patches of liquid-like molecules remains. The prism surface (bottom row) is an HT-RF phase all the way from 0 to −17 K but is very close to the roughening transition at ΔT > − 2 K in our model. In the range from −17 to −57 K, it is a DOF phase and becomes an OF phase below −57 K. At the transition from DOF-like to HT-RF phases, step free energies increase anomalously and result in the crossover of crystal growth rates. The shaded areas illustrate the temperature range where melting of full bilayers has been accomplished. Blue, less than one full bilayer; white, less than two full bilayers; yellow, more than two full bilayers.

After heterogeneous nucleation (2, 3), there is still a long way before ice crystallites adopt the micrometer size found in cirrus clouds (4). Whether the ice embryos transform into mature columnar or platelike hexagonal prisms depends on the relative growth rate of basal and prism facets (5, 6). But what drives the crossover of basal and prism growth rates, and how such changes are related to the ice surface structure, remains completely unknown to date (7, 8).

Kuroda and Lacmann (9) speculated that the occurrence of surface phase transitions could result in sudden changes of the crystal growth rate and provide an explanation for the observed crystal habits, but unfortunately, the experimental verification of this hypothesis has remained elusive to date (1014) and our current understanding is still far from providing a molecular explanation (1417).

Actually, for a long time, experiments and simulations have been aimed at measuring the onset of premelting and the premelting layer thickness (10, 11, 18), but the results have been often contradictory and showed a continuous increase of the premelting film, with little signs of a phase transition (14). Recent experiments with increased sensitivity show evidence of layering phenomena on the basal facet (13), while simulations of the monoatomi Water model (mW) emphasize also the importance of the heterogeneous in-plane structure of ice’s premelting film (19). This is a very important observation that highlights the need to probe the structure of ice’s surface not only in the perpendicular direction (layer thickness) but also, in the horizontal direction, because often parallel correlations of the surface structure hold the key to understanding surface properties and crystal growth (2023).

One of the main achievements of simple solid on solid (SOS) crystal growth models has been the identification of the roughening transition, which signals the onset of diverging parallel correlations of the surface height fluctuations of a crystal (20). Below the roughening temperature, TR, at an ordered flat (OF) phase, crystal surfaces exhibit only limited disorder, with surface height fluctuations that remain finite and congruent with the underlying bulk lattice (as in the usual Kossel crystal). At TR, a roughening transition of the Kosterlitz-Thouless type occurs, where the crystal surface unbinds from the underlying crystal mold and exhibits diverging height fluctuations, which do not differ qualitatively from those found at a fluid-fluid interface. Below TR, at the OF phase, crystal growth is a slow activated process that occurs via two-dimensional (2D) nucleation and is completely suppressed at low saturation. On the contrary, at a rough phase above TR, crystal growth is no longer activated and depends linearly on saturation (20).

This basic framework of crystal growth theory can become far more complicated in extended SOS models, where one incorporates molecular interactions at different length scales (2123). Particularly, a preroughening transition at a temperature Tpr has been documented for some models, where the OF phase “roughens,” meaning that it exhibits up and down correlated crystal step proliferation, but it is “flat” in the sense that the crystal height fluctuations remain finite (21). Dynamically, this disordered flat (DOF) phase is thought to exhibit crystal growth behavior intermediate between an OF phase and a rough phase, because growth remains activated in principle, but the step free energies are expected to decrease much as a signature of step proliferation (21, 24). The occurrence of a DOF phase is thought to be accompanied partly with layering (25), as observed in experiments on the ice basal facet (13), and with strong parallel heterogeneity, as reported for the mW model of water (19).

Here, we show that, similar to bulk ice polymorphism, the facets of ice exhibit a number of surface polymorphs as temperature changes. Particularly, we identify a transition from an OF to a DOF phase that is akin to microphase separation of terraces on the crystal surface (23, 21). Upon increasing temperature, the premelting film thickness grows and the surface structure flattens again, resulting in the unexpected increase of nucleation step free energies and the crossover of crystal growth rates exactly as found in laboratory and atmospheric field studies (8, 17).


We perform large-scale computer simulations of the TIP4P/Ice point charge model of water (26) on elongated rectangular surfaces meant to characterize large wavelength correlations along one direction. A bulk ice sample is placed in vacuum at temperatures below the model’s triple point, Tt = 272 K. After a few nanoseconds, the surface spontaneously develops a layer of quasi-liquid disordered molecules that can be readily distinguished from the underlying bulk crystal network (Fig. 2). Averaging the position of the outermost solid and liquid-like atoms of the premelting layer about points r on the plane of the interface, we are able to identify distinct ice/film (i/f) zif(r) and film/vapor (f/v) zfv(r) surfaces (27), which separate the premelting film from the bulk solid and vapor, respectively (Fig. 2 and Materials and Methods). Using these definitions, we can monitor the local thickness of the premelting film as h(r) = ∣ zfv(r) − zif(r)∣ and perform a detailed study of in-plane surface structure as a function of undercooling, ΔT = TTt. The thermally averaged thickness of the premelting layer, h, grows from about 3 Å at ΔT = − 82 K to 9 Å at ΔT = − 2 K with little measurable anisotropy (fig. S1) (18). However, as recently observed for the basal plane in the mW model (19), this thin disordered layer is laterally inhomogeneous up to about ΔT = − 9 K (Fig. 3B and fig. S2). Our study reveals that the heterogeneity is also found to a similar extent on the prism plane, with little qualitative differences (Fig. 4B and fig. S3).

Fig. 2 Characterization of the premelting layer.

After placing a slab of perfect ice in vacuum, a premelting layer of disordered water molecules is formed spontaneously in a few nanoseconds as shown in the snapshot. Using a suitable order parameter, it is possible to distinguish liquid-like (red) from solid-like (blue) molecules. The state of the premelting film may be described in terms of two different surfaces, zif(r) and zfv(r), separating the film from bulk ice and vapor phases, respectively. A film thickness, h(r), may be defined as the difference between ice/film and film/vapor surfaces such that h(r) =zfv(r) zif(r)∣.

Fig. 3 Surface fluctuations on the basal facet.

(A) Probability distribution of ice/film (blue) and film/vapor (red) surface fluctuations, as measured by the deviations of the interface position zα(x) about the average surface z¯α, for α = {if, fv}. Results are shown for different temperatures as indicated in the legend. (B) Snapshots of the basal surface at the same four temperatures. Red lines show the connected hydrogen bond network of all solid-like and liquid-like water molecules. The violet patches represent disordered liquid-like molecules. At low temperature, the surface is mainly a regular hexagonal honeycomb with a few patches of liquid-like molecules sitting on interstitial positions (as indicated by the yellow circles). The extent of filled interstitial positions increases as the premelting layer covers the surface. (C) Plot showing a snapshot of local surface height fluctuations δzif(r) (bottom, blue frame) and δzfv(r) (top, red frame) on the basal ice face. Notice the emergence of large-scale correlated patches for the DOF phase in the temperature range ΔT = −32 to −9 K (see movies S1 and S2). The patches disappear at high temperature as the surface flattens again. (D) Wave vector–dependent stiffness coefficients, as obtained from the inverse surface structure factor for i/f correlations (blue), f/v correlations (red), and crossed i/f -f/v correlations (green). Crosses are results from simulation; full lines are fits to the SG-CW model. The results show that all surfaces are smooth, as indicated by the divergence of Γ(qx) as qx → 0. Note that the sharp minimum appears at intermediate length scales in the temperature range ΔT = −32 to −9 K, where the DOF phase is present.

Fig. 4 Surface fluctuations on the prism facet.

Content displayed as in Fig. 3. (A) A bimodal distribution in this facet is only visible at temperature ΔT = −32 K. (B) Here, the snapshots show the characteristic rectangular mesh of the prism facet. At low temperature, the liquid-like molecules form patches on the solid surface as in the basal face. (C) Emergence of large correlated domains signals a DOF phase that is clearly visible at ΔT = −32 K and vanishes at notably lower temperatures than for the basal facet (see movies S4 and S5). (D) Likewise, the sharp minimum of the stiffness coefficient is visible only at and below ΔT = − 32 K.

Where does the marked crystal rate anisotropy found in experiments come from? To shed light into this problem, we quantify surface fluctuations along the long direction, x, of the simulation box. This is done by monitoring deviations, δz(x), of the local surface position, z(x), from the instantaneous average surface, z¯, for either the ice/film δzif(x)=zif(x)z¯if or the film/vapor surface δzfv(x)=zfv(x)z¯fv (Materials and Methods). This analysis is essential to reveal differences between basal and prism planes and allows us to identify a number of phase transitions along the sublimation line, which cannot possibly be inferred from visual inspection of the snapshots (Figs. 3B and 4B and figs. S2 and S3) or the premelting layer thickness (fig. S1).

Structure of the basal facet

At low temperature (ΔT = −82 K), the basal plane consists of a relatively OF solid surface, as revealed by a singly peaked, close to Gaussian distribution of both i/f and f/v surface fluctuations (Fig. 3A and fig. S4A). From the snapshots (Fig. 3B), the solid surface is formed mainly of an oxygen-unreconstructed stack of chair hexagons (2830). Patches of disordered liquid-like molecules are found and often show a tendency to sit on interstitial positions at the center of the primary hexagonal mesh (as in the so-called Honeycomb Fletcher phase) (28). The distributions remain unimodal up to ΔT = −52 K, but somewhat broaden, revealing a large increase of disorder in this temperature interval, which is consistent with observations of sum frequency generation spectroscopy (10). At a temperature ΔT = −42 K close to the onset of substantial molecular mobility of the outermost surface layer (31), the distribution of δzfv(x) develops a distinct trimodal character (fig. S4A). A main peak is centered at the mean surface position, and two other peaks appear to the left and right. The central peak of the trimodal distribution for δzfv(x) gradually fades away into a bimodal, which persists up to ΔT = −6 K. In a narrow range between −9 and −6 K, the distributions of δzfv(x) and δzif(x) are both fully bimodal and congruent. Last, at the temperature of −2 K, the bimodal collapses sharply into one single unimodal distribution.

From our analysis, the outer f/v surface of the premelting film exhibits a bimodal distribution centered at the mean surface location all the way from −22 to −6 K. The onset of bimodality very much correlates with the vanishing of the (ppp-polarized) dangling OH bond stretch observed in sum-frequency generation (SFG) spectroscopy experiments (10). The separation between peaks in the bimodal is approximately 3.1 Å, somewhat smaller than c/2 = 3.65 Å, the distance between adjacent bilayers.

Furthermore, the bimodal evolves after the appearance at low temperature of a trimodal distribution with a main peak centered at the mean surface position, as if ice melting resulted in the formation of water-like molecules at half-integer lattice positions (fig. S4A).

Preroughening and smoothening transitions

A strongly disordered phase consisting of a smooth surface with large-scale step proliferation has been documented in the literature for SOS models and is known as a DOF phase (21). Exactly as suggested by the extended SOS model, we observe that the low-temperature smooth flat phase with unimodal probability distribution transforms into a locally rough interface with a trimodal probability distribution. This change is the hallmark of a preroughening transition (21, 25). Here, the low-temperature flat phase (OF) is transformed into a DOF phase at a preroughening temperature ΔTpr ≈ −47 K. However, contrary to the usual scenario, at high temperature, the DOF phase does not undergo roughening but rather transforms across an apparently first-order surface phase transition into a new flat phase with unimodal distribution at a smoothening temperature, ΔTs ≈ −4 K (Fig. 3A). This transition has not been previously reported but is consistent with existing surface phase diagrams for the extended SOS model (22, 25). Following the notation of SOS literature, we call this a high-temperature reconstructed flat (HT-RF) phase. In practice, the nature of the actual HT-RF phase found here for the basal facet can be inferred from the snapshots (Fig. 3B and fig. S2). A premelting film appears to be fully developed, and the underlying hexagonal mesh of solid molecules is covered with interstitial water molecules located at the center of the hexagons.

As expected for the preroughening scenario predicted in extended SOS models (23, 25), the transition to a DOF phase is strongly correlated with the growth of a premelting film in a loosely layer-wise fashion. This can be shown by inspecting the distribution of laterally averaged film thickness, P(h), in Fig. 5A (see also fig. S4B). We find that the stabilization of the bimodal DOF phase observed at about ΔT = −22 K results after the full formation of a second layer of premelted ice, as revealed by a shift in the maximum of the distribution from h = 5.0 Å at ΔT = −32 K to h = 7.4 Å at ΔT = −9 K. SFG experiments reveal a large increase in the premelting film thickness in this temperature range (13). Likewise, the transition from a DOF phase into the HT-RF phase is also accompanied by the formation of a third premelted layer, as revealed by the shift in the premelting layer thickness from h = 7.4 Å at ΔT = −9 to h = 9.8 Å at −2 K (Fig. 5A).

Fig. 5 Fluctuations of premelting thickness.

Results are shown for basal (A and C) and prism (B and D) facets. (A and B) Distribution of average film thickness, h, as a function of undercooling, with ΔT = −82 K (blue), −32 K (green), −9 K (yellow), and −2 K (red). Dashed vertical lines represent full layers in units of the molecular diameter. The transition of premelting layer thickness across integer multiples of the molecular diameter occurs at broadly the same temperature for the basal and prism facets. (C and D) Surface plot of the instantaneous premelting thickness, h(r), for basal (C) and prism (D) surfaces (see movies S3 and S6). Notice the absence of large correlated domains at all temperatures, in marked contrast with the δzif(r) and δzfv(r) surfaces shown in Figs. 3 and 4.

The DOF phase is not only characterized by step proliferation and strong local disorder. Predictions from the extended SOS model also show that the steps are highly correlated, while, depending on the system, the preroughening transitions can be either second-order (accompanied by diverging parallel correlations) or first-order (with parallel correlation that remain finite) (25). We explore the extent of parallel correlations with a plot of the local surface height fluctuations of the ice/film δzif(r) and film/vapor δzfv(r) surfaces (Fig. 3C and fig. S6, A and B). At low temperature, the i/f surface exhibits small amplitude up and down domains, with small correlation lengths. At ΔT = −9 K, where the DOF phase is fully formed, we observe very large up and down domains of about 9 nm in length that remain correlated over the full simulation box, as is visible both in the figure and in the accompanying movies (movies S1 and movies S2). However, at ΔT = −2 K, just a few kelvin above, the large correlated domains vanish. Our results are consistent with glancing angle x-ray experiments, which reported the appearance of a large surface correlation length in the nanometer range, and the sharp disappearance of the long-range correlations close to the triple point (32).

The nature of these correlations can be quantified from the wave vector–dependent surface structure factor. Plots of the related effective stiffness Γαβ(qx)=kBTAzα(qx)zβ(qx)qx2 are shown in Fig. 3D (fig. S4C). The results confirm that both at the preroughening and at the smoothening transition, the parallel correlations remain finite at qx → 0, as can be inferred from the strong divergence of the effective stiffness coefficients (notice the logarithmic scale of the figure). Everywhere in the region where a DOF phase is present, however, the strongly correlated up and down domains are detectable as a sharp and deep minimum of the stiffness coefficients at intermediate wave vectors (i.e., as peaks in the surface structure factor). A full theoretical description of this strong enhancement of large but finite correlations seems difficult. However, we can definitively observe in our results how the location of the sharp minima at intermediate wave vectors decreases as the size of the correlated domains in Fig. 3C increases.

Structure of the prism facet

The study of surface fluctuations on the prism facet is substantially different. A low-temperature flat phase that preserves the expected low-temperature rectangular mesh is observed at ΔT = −82 K (Fig. 4B). At this temperature, the distribution of both zif(x) and zfv(x) is unimodal but becomes gradually broader and skewed to the left as temperature is increased (fig. S5A). At ΔT = −32 K, however, the distribution of zif(x) and zfv(x) becomes slightly bimodal. From ΔT = −12 K to ΔT = −2 K, the distribution becomes again completely unimodal and Gaussian like, but has broadened abruptly at ΔT = −2 K, signaling the approach of a roughening transition (20). Although the order parameter is not sharp enough to reveal a strong bimodality, a large number of steps are formed in the range between −62 and −22 K with distributions that seem to resemble a θDOF phase [similar to the DOF phase, but with a continuous change of the step coverage (25)]. The transition from a DOF phase to a new HT-RF phase is visible in the surface maps (Fig. 4C), where large-scale domains appear to end at −32 K. This is also visible in the surface structure factors depicted in Fig. 4D, which reveal again diverging stiffness coefficients at qx → 0 (whence, flat phase), and the complete disappearance of the sharp minimum beyond this temperature (whence, loss of long-range surface order). The transformation of the surface structure across the θDOF phase is also accompanied with continuous increase of the premelting film (Fig. 5B), but the transition from two to three layers and beyond appears to occur almost at the same temperature as in the basal facet. On the contrary, the transition from the θDOF phase to the HT-RF phases is distinctly different at the basal (ca. Ts = −4 K) and prism (ca. Ts = −27 K) facets.

We find that, on average, for both the basal and prism facets, zfv(r) follows broadly the corrugation imposed by the zif(r), and the fluctuations of premelted film thickness appear as broad unimodal distributions with no sign of bimodality or diverging correlation lengths (Fig. 5, A and B). Accordingly, there are apparently no layering transitions along the sublimation line in the temperature range studied, at least in the thermodynamic sense. The lack of large correlation lengths or bimodality is very clearly observed in the surface plots of local film thickness h(r) = zfv(r) − zif(r) (Fig. 5, C and D; see also figs. S6C and S7C and movies S3 and S6) and is consistent with SFG experiments and preliminary indications for the mW model (12, 19).

Step free energies

The importance of DOF phases has been well documented (2123). The step proliferation is akin to a strong reduction of the step free energy and the sharp decrease of the threshold for linear growth (21, 22, 24). It is therefore expected that, as temperature rises across the smoothening transition from a DOF to HT-RF phase, the crystal growth rates will decrease anomalously.

We provide a quantitative test of this expectation using a model of coupled capillary wave and sine-Gordon Hamiltonians (SG-CW) for the spectrum of surface fluctuations discussed recently (27). The capillary wave Hamiltonian describes the f/v surface fluctuations, whereas the likelihood of step proliferation in the underlying i/f surface is described with a sine-Gordon Hamiltonian. Both models are coupled with an interface potential that sets the equilibrium film thickness as the difference between zfv(r) and zif(r). The model can be fit to the regular (Gaussian) part of the surface stiffness (Figs. 3D and 4D) and provides phenomenological coefficients for the surface tension and the step free energy, β (Materials and Methods). The nonmonotonous behavior observed for β correlates with the behavior measured experimentally for terrace spreading rates (5) and is qualitatively similar to step free energies obtained from growth measurements on snow crystals (Fig. 6) (17). For the prism facet, β increases sharply in the range ΔT = 20 to 30 K, somewhat below the estimated smoothening temperature. For the basal facet, on the other hand, β increases in the range between 9 and 6 K, very close to the estimated smoothening transition. Although the increase for the basal face is small, it is sufficient to trigger the crossover of step free energies at about ΔT = −6 K, in almost exact agreement with the crossover observed in the experimental results (Fig. 6). Accordingly, our results lend support to the suggested scenario of anomalous increase of β in the neighborhood of the smoothening transitions of basal and prism facets, respectively.

Fig. 6 Calculated step free energies.

Step free energies as obtained from a fit of the mean field SG-CW model to the regular (Gaussian) part of the stiffness coefficients in Figs. 3D and 4D for the basal (blue) and prism (red) facets. Results (squares with dashed lines) are compared to experimental data (circles with full lines) (17) and displayed on two different figures to avoid crowding. Notice that the scale in both figures is the same. The step free energies from the fit exhibit a crossover at ca. −6 K, 2 K above that found in the accepted Nakaya diagram, and display overall a trend similar to the experimental data.


Origin of the anomalous behavior

The amazing polymorphism of ice crystals results from several factors that are unique to water. (i) A subtle balance between proton ordering, which is usually favored enthalpically, and proton disordering, which is favored entropically; (ii) a large number of oxygen lattices consistent with hydrogen bonding; and (iii) the allowance of interpenetrated hydrogen bond networks. A combination of these factors appears to operate also at the ice surface, which, contrary to bulk solid phases, can also gain additional entropy by disruption of the ordered oxygen network.

Below 200 K, quantum mechanical calculations indicate that both basal and prism facets of pristine ice form stripes of alternating up and down dangling hydrogen bonds in a low-temperature highly ordered Fletcher phase (14, 29, 30). As temperature increases, entropy can be gained by disordering the hydrogen bond network while keeping oxygen positions fixed or, alternatively, by disrupting the oxygen framework altogether to gain translational entropy. The structural evolution of the ice surface as temperature increases results from an intricate interplay of these two factors, which appear to participate successively.

As the OF phase is heated, the surface free energy can be minimized by disordering the dangling hydrogen bond network, at the cost of breaking the stabilizing Fletcher stripes. This loosens the ice surface, resulting in the decrease of step free energies observed at low temperature (see Fig. 6). Once the surface hydrogen bond network has achieved full disorder, however, no further entropy gain can be achieved without much disruption of the oxygen lattice. This provokes the gradual melting of a bilayer and increases the premelting layer thickness (13, 18, 31). For thin films, this occurs with minimal enthalpy cost because the premelted water molecules can shift to interstitial sites corresponding to an interpenetrated crystal network (see Fig. 3C).

As large amounts of translational entropy are gained by the formation of a thin premelting layer, the surface enthalpy can be lowered by partially building up the hydrogen bond order again. This results in an increase of the local surface stiffness and the rise of the step free energy at ca. ΔT = −30 K for the basal face and ca. ΔT = −20 K for the prism face observed in Fig. 6.

With hydrogen bond order building up, it is expected that local domains of enhanced surface polarity can result (14). These domains are stabilized locally because of the hydrogen bond connectivity of the underlying bulk crystal network. However, at large distance, domains of equal polarity repel each other. The combination of short-range attractive and long-range repulsive interactions is precisely the main ingredient in the stabilization of a DOF phase (21, 22). The entropy cost of local ordering the hydrogen bond network is balanced by the presence of substantial translational disorder of the premelting film and the meandering entropy associated with the moving boundaries of the domains. The gain in entropy is sufficient to drive the line tension between domains effectively negative and promote the spontaneous formation of steps on the surface (19).

Similarly, in colloidal science, colloids with short-range attractive and long-range repulsive forces are known to microphase separate. Therefore, the transformation of the low-temperature stripe phase into a DOF phase with compact clusters of different height can be viewed as akin to the transition from lamellar to cluster phases in colloids with competing interactions.

However, the topology of the basal and prism surfaces dictates very different enthalpy costs for the boundary between different domains. For the stripe phase of the basal facet, grain boundaries between locally ordered domains have a very low energy penalty (28). For the prism face, on the other hand, the cost of surface grain boundaries between ordered domains is very high (28) and does not allow for the stabilization of the DOF phase in such a large temperature range.

Eventually, as the triple point is approached, a new bilayer melts and the meandering entropy of the basal face is traded for the translational entropy gained by the additional layer in the premelted film. In this situation, the ice/film surface is similar to the ice/water interface, and step free energies approach the value expected for ice freezing from water at the triple point (9). There, the basal face is faceted, and has a finite step free energy, but the prism face is rough and has a vanishing step free energy (27), leading to the observed crossover of step free energies at high temperature (see Fig. 6).

Relation to crystal growth

Ice crystallites in the atmosphere are expected to grow by a 2D nucleation process, with crystal growth rate, j, given by (9)j=A(T)Smexp (gmβ2akB2T2S)(1)where A(T) is a temperature-dependent prefactor, S = ln p/pc is the saturation, a is the surface area per molecule, gm is a numerical factor, kB is Boltzmann’s constant, p is the ambient water vapor pressure, and pc is the ice/vapor saturation pressure.

For a mononuclear mechanism, the exponent m is unity, and gm = π, whereas for polynuclear growth, m = 2/3 and gm = π/3. In the nucleated regime at constant ambient conditions, the crystal growth rate is therefore dictated essentially by β. If the step free energy of the basal facet is larger than that of the prism facet, the basal face grows slower than the prism face and platelike crystals form. On the contrary, if the basal step free energy is smaller, the basal face grows faster and columnar crystals form (7, 9, 17).

Close to the melting point, our results show that the prism facet has smaller step free energies than the basal facet (Fig. 6). In this range, accordingly, growth of plates rather than columns is predicted to occur, in agreement with experiments (see Fig. 1). As the temperature decreases, our calculated step free energies cross over twice, predicting correctly the sequence of transitions from plates to columns and back to columns that is the most striking feature in the Nakaya diagram. The first crossover, corresponding to the transition from plates to columns, takes place at ΔT = −6 K, which is 2 K away from the value of ΔT = −4 K observed in experiments. The next crossover is predicted to occur at the temperature ΔT = −26 K, which is shifted by about 15 K from the value of ΔT = −10 K in the Nakaya diagram.

This suggests that the precise location of the remaining crystal habit transitions could depend not only on step free energies but also on additional factors, such as the surface mobility of premelted molecules that is embodied in the prefactor, A(T), or the diffusion limited vapor field that surrounds the crystal (9).

In conclusion, we show that in the range of about 80 K below the melting point, the main facets of ice may be found in up to three different surface phases with varying degree of surface disorder, as postulated by Kuroda and Lacmann with amazing intuition almost 40 years ago (9). The accompanying phase transitions provide the mechanism for a nonmonotonous change of the relative step free energies for 2D nucleation. Most notably, we observe a premelting mediated process of surface smoothening, which results in the anomalous increase of step free energies. This results in the crossover of relative growth rates of basal and prism facets that is required to explain the Nakaya diagram.

The explanation of the long-standing problem of snow crystal shapes (Fig. 1) proofs that we have now a close to complete molecular description of the surface structure and crystal growth rates of ice in the atmosphere, with immediate implications in atmospheric sciences, glaciology, and climate modeling.


Force field

Our study is performed with the TIP4P/Ice model of water (26). This model was purposely designed to best describe the properties of ice. It predicts a melting point of T = 272 K, in excellent agreement with experiment, and reproduces the most relevant surface properties at this temperature, such as liquid-vapor surface tension (γlv = 82 mN/m calculated by ourselves, compared with γlv = 75.7 mN/m from experiment) and solid-liquid surface tension [γsl = 29.8 mN/m from (33) compared with recommended results γsl = 28 mN/m by Pruppacher and Klett (34)].

The precise location of the surface phase transitions observed here could somewhat change with the molecular model used, but we expect the generic features to be quite generally observed for other accurate intermolecular potentials.

Initial configurations

Initial configurations are prepared from a perfect unit cell in pseudo-orthorhombic arrangement, consisting of two layers of hexagonal rings perpendicular to the hexagonal c axis and a total of 16 water molecules. For the basal interface, we arrange a stack of 46 × 8 × 8 cells of 47,104 molecules, with the long direction, x, aligned along the b axis of the pseudo-orthorhombic cell, corresponding to the so-called (basal)[pII] surface arrangement described by Davidchack et al. and used in our previous work (27, 35, 36).

For the prism facet (pI), the simulation box is prepared from a stack of 40 × 8 × 8 unit cells of N = 40,960 molecules, with the long direction, x, aligned along the a axis. This corresponds to the (pI)[basal] arrangement in our recent work. For each such arrangement, we prepare an independent hydrogen bond network as described in (37). Using this method, a hydrogen bond network satisfying the Bernal-Fowler rules is selected randomly on a subsystem. This arrangement is replicated four times to obtain the full system. Initial configurations are selected such that the net dipole moment is smaller than 0.4% of one single TIP4P/Ice molecular dipole moment. Notice that in a bulk simulation of ice, hydrogen bond networks are arrested and no further sampling is possible. However, for simulations where the bulk ice phase terminates at a free boundary, as here, ring rotations take place spontaneously during the molecular dynamics simulations and the sampling is performed over many realizations of the hydrogen bond disorder.

After forming the ice slab, we perform NpT simulations of the bulk solid at the desired temperature to obtain equilibrated unit cell dimensions. The solid is then scaled to the average equilibrium cell value, placed in vacuum, and equilibrated again in the canonical ensemble under periodic boundary conditions.

Computation details

Large-scale simulations are carried out on the MareNostrum 4 facility at Barcelona Supercomputing Center from the Spanish National Supercomputing Network. Classical molecular dynamics simulations are performed with the 2016.4 version of the GROMACS package. Trajectories are evolved with the GROMACS default leap-frog algorithm. Both the Lennard-Jones and Coulomb interactions are truncated at 0.9 nm. The electrostatic interactions are calculated using the particle mesh Ewald method. Simulations are thermostated with the velocity rescale algorithm (38) and the Berendsen barostat when required. A relaxation time of 2 ps is used for both the thermostat and barostat. The time step used is 3 fs. This is a somewhat larger time step than 2 fs usually used in thermostated dynamics of TIP4P models, but we have checked that energy is well conserved for time lengths sufficiently larger than the thermostat relaxation time of 2 ps. We believe that the need to perform reliable ensemble averages with correlation lengths in the scale of several many nanometers warrants an increase in the time step. For symplectic algorithms, it is known that this effectively provides the correct dynamics for an effective Hamiltonian with corrections of order dt2. For this reason, the transitions that we find could be shifted by 1 or 2 K from those obtained for the exact Hamiltonian. Three runs (300 ns long) are launched for each temperature and crystal facet starting from an independent hydrogen bond network. We discard the first 75 ns of each run and perform averages over the remaining 225 ns. Accordingly, our results report thermal averages collected over a total of 675 ns for each temperature and facet.

Surface analysis

Order parameter. Before determining the i/f and f/v surfaces, we label water molecules as either solid or liquid like, using the q¯6 parameter (39). Water-like molecules are those with a q¯6 parameter below a threshold q¯6*(T). To determine the threshold, we simulate the probability distribution of q¯6 at a number of temperatures in either bulk solid or liquid water. The threshold value q¯6*(T) is determined such that the number of mislabeled liquid molecules on the solid phase is equal to the number of mislabeled solid molecules on the liquid phase. Previously, we have shown that using these criteria, we are allocating ice allotropes (ice clathrate and ice Ic environments) as well as surface undercoordinated molecules (mainly interfacial ice Ih) that appear at the ice surface as belonging to the premelting layer. Whence, our ice-liquid surface separates the strict hexagonal ice template from other less ordered interfacial environments. We have found that the number of molecules allocated to ice clathrate, ice Ic, and interfacial ice Ih environments does not appear to change much in the temperature interval studied here so that the premelting layer fluctuations are hardly affected by this choice (40).

Surface location. A cluster analysis is performed to determine which molecules pertain to the condensed phase. Water molecules with oxygen atoms at a distance less than 3.5 Å belong to the same cluster. The i/f surface is determined from the positions of solid-like atoms in the largest solid cluster. We use the heights of the four topmost (or bottommost) solid atoms of the upper (lower) interface. At a given point r on the plane of the interface, we find all the solid-like atoms lying within a rectangular prism centered at r. The base of the prism is half the area of the pseudo-orthorhombic unit cell. The surface height zif(r) at that point is determined from the average location of the four uppermost solid-like atoms. At the same point, the liquid surface for the upper (lower) interface is determined by averaging the position of the uppermost (bottommost) four liquid-like molecules of the cluster of condensed molecules lying within a rectangular area of 3σ × 3σ Lennard-Jones molecular diameters. The surfaces zif(r) and zfv(r) are determined over points on a grid on the plane of the interface. The grid has twice as many points as unit cells along the x direction and just as many points as unit cells along the y direction. We perform the surface analysis from the set {zif(r)} and {zfv(r)} of points on the grid (27). The instantaneous mean position of the i/f surface z¯if is determined as the lateral average of {zif(r)} overall points on the grid. From this value, we obtain δzif(r)=zif(r)z¯if. The laterally averaged fluctuations δzif(x) are obtained from δzif(r) upon averaging along points y. δzfv(r) and δzfv(x) are obtained likewise. Fourier transforms of δzif(x) and δzfv(x) are obtained by summing δzα(x)eiqxx over all points along x. We have shown previously that the fluctuations of i/f and f/v surfaces so defined consistently provide the expected ice/water and water/vapor surface tensions as the premelting layer thickness grows on approaching the triple point (27).

Instantaneous local film heights are obtained as h(r) = zif(r) − zfv(r). The instantaneous average film thickness is obtained from the mean of h(r) over the points of the grid.

SG-CW model and fit

We describe the coupled i/f and f/v surface fluctuations with an extended sine-Gordon model for the i/f surface and the capillary wave model for the f/v surface (20). The two terms are coupled via the interface potential, g(h), with the premelting film thickness given by h(r) = zfv(r) − zif(r). The full Hamiltonian is given byH=dr[γ˜iw2(zif)2+γwv2(zfv)2+γivzifzfvucos (qzzif)+g(zfvzif)](2)where γiw is the bare stiffness coefficient, γwv is the water-vapor surface tension, γiv dictates the coupling of surface deformations, u accounts for the cost of moving the surface zif away from integer lattice spacing, g(h) is the interface potential dictating the equilibrium film thickness, and qz is the wave vector for a wavelength equal to the lattice spacing. This Hamiltonian can be expanded to quadratic order in deviations away from the mean surface positions z¯if and z¯fv and yields for the thermally averaged surface fluctuations (27)zif2(q)=kBTAg+γwvq2[υ+g+γ˜iwq2][g+γwvq2][g+Δgγivq2]2zfv2(q)=kBTAυ+g+γ˜iwq2[υ+g+γ˜iwq2][g+γwvq2][g+Δgγivq2]2zif(q)zfv*(q)=kBTAg+Δgγivq2[υ+g+γ˜iwq2][g+γwvq2][g+Δgγivq2]2(3)where g′′ is the second derivative of the interface potential at the equilibrium film thickness, Δg′′ accounts for enhanced coupled compression-expansion of the film thickness, and υ=qz2u. In practice, to avoid underdetermined fits to limited data, we set Δg′′ and γiv to zero. This model for the spectrum of surface fluctuations has small and large wave-vector regimes. At large wave vectors, the spectrum of fluctuations depends only on the stiffness and surface tension coefficients, γiw and γwv, respectively (36). As in extended capillary wave Hamiltonians, these are modeled as even polynomials of q to order 4. Once γiw and γwv are known, we fit the remaining parameters g′′ and υ to match simultaneously the low wave-vector regime of q2zif(q)zif(q)〉, q2zfv(q)zfv(q)〉, and q2zif(q)zfv(q)〉 as obtained from simulations. The model reproduces very accurately the stiffness coefficients at high temperature. In the region where the DOF phase appears, it does, however, not grasp the enhanced height fluctuations. However, quite generally, the correlation functions may be expressed as a spectral series. The mean field solution corresponds to the “trivial” eigenmode, which results for the quadratic expansion of Hamiltonian. Higher-order solutions yield additional eigenmodes. For strictly periodical height profiles, the spectrum of eigenmodes has a band structure. For the quasi-long range observed in the DOF phase, we expect the spectrum will not have continuous bands, but rather a region of close eigenvalues that are separated from the mean field mode. Accordingly, the effect of the inhomogeneity appears as an additive correction to the mean field solution. This allows us to obtain meaningful fits to the regular part of the spectrum merely by neglecting the strong enhancement at intermediate length scales. The step free energy for the uncoupled sine-Gordon model can be obtained from the parameters γiw and υ as βiw=(8/qz2)(γiwυ)1/2 (36, 41). Because we find that the premelting film thickness is unimodal, a step on the i/f surface will provoke a similar step on the f/v surface, with a step energy of βwv=(8/qz2)(γwvg)1/2. The cost of creating a step on the premelting film is therefore the sum β = βiw + βwv, which we use as an estimate for the step free energy of the ice/vapor interface in Fig. 6.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Acknowledgments: We thank E. Lomba for helpful discussions and J. L. F. Abascal for invaluable support. Funding: We acknowledge funds from the Spanish Agencia Estatal de Investigación under grant no. FIS2017-89361-C3-2-P. This project was made possible because of the use of the MareNostrum supercomputer and the technical support provided by Barcelona Supercomputing Center from the Spanish Network of Supercomputing (RES) under grants QCM-2017-2-0008 and QCM-2017-3-0034. Author contributions: P.L. performed programs and simulations, analyzed the data, and discussed results as part of his Ph.D. thesis. E.G.N. supervised and discussed results. L.G.M. conceived the project, interpreted results, and wrote the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article