## Abstract

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.

## INTRODUCTION

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).

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 (*10*–*14*) and our current understanding is still far from providing a molecular explanation (*14*–*17*).

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 (*20*–*23*).

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, *T*_{R}, 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 *T*_{R}, 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 *T*_{R}, 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 *T*_{R}, 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 (*21*–*23*). Particularly, a preroughening transition at a temperature *T*_{pr} 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*).

## RESULTS

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, *T*_{t} = 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) *z*_{if}(**r**) and film/vapor (f/v) *z*_{fv}(**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**) = ∣ *z*_{fv}(**r**) − *z*_{if}(**r**)∣ and perform a detailed study of in-plane surface structure as a function of undercooling, Δ*T* = *T* − *T*_{t}. 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).

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,

### 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 (*28*–*30*). 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 δ*z*_{fv}(*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 δ*z*_{fv}(*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 δ*z*_{fv}(*x*) and δ*z*_{if}(*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 Δ*T*_{pr} ≈ −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, Δ*T*_{s} ≈ −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).

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 δ*z*_{if}(**r**) and film/vapor δ*z*_{fv}(**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 *q _{x}* → 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 *z*_{if}(*x*) and *z*_{fv}(*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 *z*_{if}(*x*) and *z*_{fv}(*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 *q _{x}* → 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.

*T*

_{s}= −4 K) and prism (ca.

*T*

_{s}= −27 K) facets.

We find that, on average, for both the basal and prism facets, *z*_{fv}(**r**) follows broadly the corrugation imposed by the *z*_{if}(**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**) = *z*_{fv}(**r**) − *z*_{if}(**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 (*21*–*23*). 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 *z*_{fv}(**r**) and *z*_{if}(**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.

## DISCUSSION

### 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*)*A*(*T*) is a temperature-dependent prefactor, *S* = ln *p*/*p*_{c} is the saturation, *a* is the surface area per molecule, *g _{m}* is a numerical factor,

*k*

_{B}is Boltzmann’s constant,

*p*is the ambient water vapor pressure, and

*p*

_{c}is the ice/vapor saturation pressure.

For a mononuclear mechanism, the exponent *m* is unity, and *g _{m}* = π, whereas for polynuclear growth,

*m*= 2/3 and

*g*= π/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 (

_{m}*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.

## MATERIALS AND METHODS

### 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 γ

*= 75.7 mN/m from experiment) and solid-liquid surface tension [γ*

_{lv}*= 29.8 mN/m from (*

_{sl}*33*) compared with recommended results γ

*= 28 mN/m by Pruppacher and Klett (*

_{sl}*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 *dt*^{2}. 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 *39*). Water-like molecules are those with a *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 *z*_{if}(**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 *z*_{if}(**r**) and *z*_{fv}(**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 {*z*_{if}(**r**)} and {*z*_{fv}(**r**)} of points on the grid (*27*). The instantaneous mean position of the i/f surface *z*_{if}(**r**)} overall points on the grid. From this value, we obtain *z*_{if}(*x*) are obtained from δ*z*_{if}(**r**) upon averaging along points *y*. δ*z*_{fv}(**r**) and δ*z*_{fv}(*x*) are obtained likewise. Fourier transforms of δ*z*_{if}(*x*) and δ*z*_{fv}(*x*) are obtained by summing δ*z*_{α}(*x*)*e*^{−iqxx} 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**) = *z*_{if}(**r**) − *z*_{fv}(**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**) = *z*_{fv}(**r**) − *z*_{if}(**r**). The full Hamiltonian is given by* _{wv}* is the water-vapor surface tension, γ

*dictates the coupling of surface deformations,*

_{iv}*u*accounts for the cost of moving the surface

*z*

_{if}away from integer lattice spacing,

*g*(

*h*) is the interface potential dictating the equilibrium film thickness, and

*q*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}*27*)

*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

*g*′′ and γ

*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,*

_{iv}*, respectively (*

_{wv}*36*). As in extended capillary wave Hamiltonians, these are modeled as even polynomials of

*q*to order 4. Once

*are known, we fit the remaining parameters*

_{wv}*g*′′ and υ to match simultaneously the low wave-vector regime of

*q*

^{2}〈

*z*

_{if}(

*q*)

*z*

_{if}(

*q*)〉,

*q*

^{2}〈

*z*

_{fv}(

*q*)

*z*

_{fv}(

*q*)〉, and

*q*

^{2}〈

*z*

_{if}(

*q*)

*z*

_{fv}(

*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

*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

*+ β*

_{iw}*, which we use as an estimate for the step free energy of the ice/vapor interface in Fig. 6.*

_{wv}## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/21/eaay9322/DC1

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.

## REFERENCES AND NOTES

**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.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution License 4.0 (CC BY).