Interface potential for water on ice

Most experiments in the literature report premelting layer thicknesses as a function of temperature. However, premelting can also be understood as the condensation of water vapor onto the bulk ice surface. Viewed as an adsorption problem, one sees that the layer thickness is both a function of temperature and vapor pressure²⁵. Strictly, ice in contact with water vapor can only be in equilibrium along the sublimation line. It follows that the premelting thickness away from the sublimation line can only be meaningfully characterized for small deviations away from coexistence, where vapor condensation and freezing occur at exactly the same rate. Ice can then be out of equilibrium, while the premelting film remains in a stationary state of constant thickness³⁸. The failure to recognize this important point is the source of much confusion in the literature and largely explains why results for the premelting layer thickness differ by orders of magnitude close to the triple point.

Here, we show that an analysis of equilibrium surface fluctuations of ice along the sublimation line can be exploited to calculate an approximate interface potential for the premelting film. Input in a suitable theory of crystal growth dynamics, this allows us to characterize the premelting layer thickness at arbitrary temperature and pressure.

To see this, we write the effective surface free energy per unit surface area at solid/vapor coexistence as ω(h;T) = g(h;T) − Δp_lv(T)h, where Δp_lv(T) is the pressure difference between the liquid and vapor bulk phases at the solid/vapor coexistence chemical potential. The free energy ω(h;T) may be calculated over a limited range of h, by simulating at solid/vapor equilibrium. During the course of the simulation, the film thickness fluctuates according to P(h;T), a probability distribution which can be easily collected. This can be used to obtain the free energy from the standard fluctuation formula \(A\,\omega (h;\!T)=-{k}_{{\rm{B}}}T\,{\mathrm{ln}}\,P(h;\!T)\), where k_B is Boltzmann’s constant and A is the surface area^39,40. On the other hand, Δp_lv(T) is a purely bulk property and can be readily calculated by thermodynamic integration from available data (see “Methods” and ref. ⁴¹). With both ω(h; T) and Δp_lv(T) at hand, a batch of simulations along the simulation line can provide g(h;T) = ω(h;T) + Δp_lv(T)h for a set of temperatures over a range of overlapping film thicknesses. Since the interface potential is expected to exhibit only a small temperature dependence, the set of piecewise functions g(h;T) at different temperatures may be combined to build a master curve g(h) over the whole range of film thicknesses spanned in the temperature interval of the simulations (see “Methods” and Supplementary Note 1).

In principle, computer simulations of ice premelting are extremely challenging. The environment of a given molecule changes from solid to liquid and then to vapor over the scale of just one nanometer or less. The local polarization changes significantly across the interface, and therefore the average many-body forces differ greatly depending on the local position. Such a complicated situation is best described with electronic quantum-mechanical calculations, or explicit many-body potentials^42,43. Unfortunately, simulations with this level of detail for system sizes as large as required here appear unfeasible. Therefore, we employ the TIP4P/Ice model⁴⁴. Although this is a point-charge non-polarizable potential, it predicts accurately both the solid/liquid and liquid/vapor surface tensions⁴⁵. Furthermore, in the range between 210 and 271 K, it produces film thicknesses that lie between 3 and 10 Å, consistent with a growing body of evidence from experimental probes^17,18,46.

The results obtained with the TIP4P/Ice model for thicknesses up to one nanometer are analyzed as described above to produce the interface potential shown in Fig. 1.

In practice, the equilibrium film thickness can grow well beyond one nanometer as the triple point is approached, so that a complete model of the interface potential requires additional input from theory and experiment.

Mean-field liquid state theory shows that a short-range contribution of the interface potential originating from molecular correlations in the adsorbed film obeys the following equation^47,48,49,50:

where C_i are positive constants, κ₁ and κ₂ are inverse decay lengths (whichever is shorter is the inverse bulk correlation length), and \({q}_{0}\approx 2\pi \ {d}_{0}^{-1}\), where d₀ is the molecular diameter.

In practice, small amplitude capillary wave fluctuations at both the solid/liquid and liquid/vapor interfaces considerably wash away the oscillatory behavior and renormalize the mean-field coefficients. Our computer simulations for the interface potential of the basal face are consistent with this scenario: fits describe the simulations accurately up to 10 Å, and then predict a fast decay with very weak oscillations of the sinusoidal term (c.f.⁴¹).

In addition, there are algebraically decaying contributions to the interface potential which stem from the long-range van der Waals interactions. These forces originate from fluctuations of the electromagnetic field. Elbaum and Schick⁵¹ parameterized the dielectric response of ice and water to numerically calculate these contributions with Dzyaloshinskii–Lifshitz–Pitaivesky theory. Following ref. ⁵², we show that the resulting crossover of retarded to non-retarded interactions is given accurately as

where f is a parameter that accounts for the relative weight of infrared and ultraviolet contributions to the van der Waals forces, a is a wavenumber in the ultraviolet region, while b falls in the extreme ultraviolet and accounts for the suppression of high-frequency contributions (see Supplementary Note 2 for further details).

The algebraic decay of the van der Waals forces provides a negative contribution to the interface potential and produces an absolute minimum at finite thickness^41,51. This explains the observation of water droplets formed on the ice surface just a few Kelvin away from the triple point^15,34,53. The droplets observed in the experiment have a small contact angle of θ ∼ 2°, which imply a shallow primary minimum with energy \({\gamma }_{{\rm{lv}}}(\cos \theta -1) \sim -1{0}^{-5}\) J m⁻².

Combining all this information, we obtain g(h) = g_sr(h) + g_vdw(h) and fit our computer simulation results to this form, with C_i, κ_i, q₀, and α as fit parameters (Supplementary Table 1 and Supplementary Note 3). In fact, the simulation results can be fitted very accurately to g_sr(h) alone⁴¹, but the extrapolation of the simulation results to larger h is required to describe the behavior at saturation. Therefore, in the parameter search we impose that g(h) exhibits minima at energies ∼−10⁻⁵ J m⁻², as observed in experiment³⁴. The constrained fit yields an interface potential in good agreement with the available simulation data—see Fig. 1. Consistent with expectations from renormalization theory, the shallow minima in the interface potential are more widely spaced than one would expect from mean-field theory, located at h_α = 16.0 Å and h_β = 24.5 Å. We refer to these two as the α- and β-minima, respectively, and this interface potential provides a transition between a thin α and a thick β film at sufficiently large supersaturation as suggested in experiments of ice premelting in the basal facet^34,53.

Interface Hamiltonian

The interface potential is adequate for describing the equilibrium properties of homogeneous films. However, in order to account for droplets like that depicted in Fig. 2 and other such inhomogeneities, we must extend our description. Building on previous work^20,45, we begin by constructing a coarse-grained free energy (effective Hamiltonian) with all the required physics, consisting of a coupled sine-Gordon plus Capillary Wave (SG+CW) Hamiltonian with bulk fields,

The first two terms account for the free energy cost to increase the surface area of the solid/liquid and liquid/vapor surfaces in a long-wave approximation, where L_sl and L_lv are the height profiles of the two interfaces, defined as the distances from the solid–liquid and liquid–vapor interfaces to an arbitrary reference plane that is parallel to the plane of the average ice surface (c.f. Fig. 2). Furthermore, γ_sl and γ_lv are the solid/liquid interfacial stiffness coefficient and the surface tension, respectively. The cosine term accounts for the energy cost, u, to move the solid/liquid surface L_sl away from the equilibrium lattice spacing, as dictated by the wave-vector \({q}_{z}=2\pi \ {d}_{{\rm{B}}}^{-1}\), where d_B is the lattice spacing between ice bilayers at the basal face. This simple model is known to describe adequately nucleated, spiral, and linear growth^54,55,56,57. The interface potential coupling the two surfaces seeks to enforce the equilibrium thickness of the premelting film h = L_lv − L_sl. The last two terms account for the bulk energy of the system as measured relative to the (reservoir) vapor phase with fixed chemical potential μ, where Δp_sl = p_s(μ) − p_l(μ) is the pressure difference between the bulk solid and liquid phases, while Δp_lv = p_l(μ) − p_v(μ) is the pressure difference between the bulk liquid and vapor phases. These two terms account for the free energy change due to growth/melting of the solid phase at the expense of the premelting film, and exchange of matter between the latter and the vapor via condensation/evaporation.

Note that the spectrum of equilibrium surface fluctuations of Eq. (3) can be obtained exactly up to Gaussian renormalization²⁰. Accordingly, the parameters required in the theory can be obtained in principle by requiring that the spectrum of fluctuations from the theory match the results from experiments or simulations^23,45. By virtue of this mapping, the input of Eq. (3) is averaged over fluctuations, so that Ω is to be interpreted as a renormalized free energy, which incorporates consistently all surface fluctuations in the scale of the parallel correlation length.

Gradient driven dynamics

The motion of the solid/vapor interface in the presence of a premelting film necessitates us to account explicitly for the different dynamical processes occurring at both the solid/liquid and liquid/vapor surfaces^24,25,26. On the one hand, L_sl evolves as a result of freezing/melting at the solid/liquid surface, and on the other hand, L_lv evolves as a result of both the condensation/evaporation at the liquid/vapor surface and freezing/melting at the solid/liquid surface. Finally, we must also account for advective fluxes of the premelting film over the surface. In practice, since we are concerned only with small deviations away from equilibrium, we can assume the dynamics is mainly driven by free energy gradients with respect to the relevant order parameters⁵⁸. Accordingly, we treat the freezing/melting and condensation/evaporation in terms of non-conserved gradient dynamics, and the advective fluid dynamics of the premelting film using a thin-film (lubrication) approximation, whence

where k_sl and k_lv are kinetic growth coefficients that determine the rate of crystallization and condensation at the solid/liquid and liquid/vapor surfaces, respectively, η is the viscosity in the liquid film and Δρ = ρ_s − ρ_l, where ρ_s and ρ_l are the densities of the solid and liquid, respectively. Models with some similar features were developed in ref. ^58,59,60.

Notice that the deterministic dynamics given by Eq. (4) is driven by the renormalized free energy, Eq. (3). Accordingly, the equation accounts for stochastic fluctuations implicitly, and it may be interpreted as dictating the evolution of the film profiles averaged over all possible random trajectories⁶¹. Alternatively, replacing the renormalized free energy by a mean-field Hamiltonian, one can assume the above result describes the most likely path of the system³². When the fluctuations are small, the coarse-grained Hamiltonian and the renormalized free energy do not differ significantly, and the evolution of the average trajectory becomes the same as the most likely path, as expected in mean-field theory. In Supplementary Note 4, we provide an extended discussion on this issue and show that Eq. (4) may be derived from a fully stochastic-driven dynamics of the mean-field Hamiltonian.

Kinetic phase diagram

The time evolution predicted by Eqs. (3) and (4) is extremely rich and varied and the full range can only be obtained numerically. However, if we assume that the surface is on average flat, then we obtain equations that enable us to predict the outcome of the numerical simulations and determine an accurate kinetic phase diagram. Coarse graining the evolution over the time period required to form a single new plane of the crystal, we replace the time derivatives of L_sl and L_lv by their average values (denoted as 〈⋅〉), yielding a rate law for continuous growth (Supplementary Note 5):

where w = q_zu, ϕ_sl = Δp_sl − Π, ϕ_lv = Δp_lv + Π and the disjoining pressure is defined as Π(h) = −∂_hg(h). In Eq. (5a), the plus sign corresponds to freezing (ϕ_sl > 0), while the minus sign corresponds to sublimation (ϕ_sl < 0).

Despite the coarse graining, Eq. (5) predict a complex dynamics in very good agreement with the numerical solutions of Eq. (4) (see below).

First, for points in the temperature–pressure plane where \({\phi }_{{\rm{sl}}}^{2}\,<\,{w}^{2}\), the crystal surface is pinned by the bulk crystal field and remains smooth. Within this region, continuous growth cannot occur. Instead, the loci of points obeying \({\phi }_{{\rm{sl}}}^{2}={w}^{2}\) encloses a region of activated growth, where the crystal front advances via nucleation and spread of new terraces^55,56.

For state points where \({\phi }_{{\rm{sl}}}^{2}\,> \,{w}^{2}\), the thermodynamic driving force becomes larger than the pinning field. The surface then undergoes kinetic roughening, and growth can proceed continuously. The growth of the premelting film thickness may be found by subtracting the growth rate of 〈∂_tL_sl〉 from that of 〈∂_tL_lv〉, yielding:

In practice, we are interested in mapping the phase diagram for quasi-stationary states, where the solid and liquid phases grow at the same rate, so that the premelting film thickness remains constant, i.e., such that \(\langle \frac{\partial h}{\partial t}\rangle =0\)^25,26. Solving for this equality provides a condition for the film thickness as a function of pressure and temperature, which is conveniently written as:

where Π(h) is the disjoining pressure, while Δp_k(p_v, T) is a function of the ambient conditions as set by p_v, but depends parametrically also on the growth mechanism and rate constants (see Supplementary Note 5).

To illustrate the significance of this equation, consider the simple case of a rough surface, i.e., such that w = 0. Then, solving Eq. (6) for stationarity, readily yields Eq. (7), with the kinetic overpressure given in the simple form:

Notice that Δp_sl and Δp_lv are purely bulk quantities that only depend on the imposed thermodynamic conditions of the system, and convey the state-dependent information to the kinetic overpressure (Supplementary Table 2 and Supplementary Note 6). In the limiting case where the substrate is strictly inert, k_sl = 0, then Eq. (8) becomes Π(h) = −Δp_lv exactly, which is the Derjaguin condition for the equilibrium film thickness on inert substrates. This is very convenient, because we can then predict the outcome of the non-equilibrium dynamics by analogy with the behavior of equilibrium films on inert substrates, albeit with the effective overpressure Δp_k replacing Δp_lv. Likewise, one sees that an effective interface potential ω_k(h) = g(h) − Δp_kh determines the dynamics of the system in the quasi-stationary regime.

This allows us to determine the kinetic phase diagram, identifying the regions in (p, T) space where the different outcomes of the interfacial wetting dynamics is to be expected (Fig. 3). In particular, we identify three significant kinetic phase lines:

• The line of kinetic coexistence (dotted-red line in Fig. 3) occurs when Δp_k = 0. The location of this line can be obtained from Eq. (7), for the choice Π(h) = 0. States above this line are effectively oversaturated and have stationary film thicknesses consistent with Π(h) < 0 and are effectively oversaturated. Accordingly, the Laplace condition for droplet formation is met for the first time, and droplets can be stabilized transiently. However, this occurs well above the liquid–vapor coexistence line, and explains why droplets reported in experiment are formed only above the condensation line^34,53.

  

  Panel a shows the equilibrium phase diagram and kinetic phase lines. The red solid line is the sublimation line, whereas the dashed lines are metastable prolongations of the vaporization (blue) and melting (black) lines. The filled triangle (▴) indicates the triple point where these lines meet. The remaining features describe the outcome of the dynamics. The shaded area designates the region of activated growth. The dotted lines are kinetic phase lines corresponding to kinetic coexistence (red dotted), kinetic α → β transition (blue dotted) and kinetic spinodal (green dotted) lines as explained in the text. Panel b shows sketches with the dynamics observed in different points of the phase diagram, as indicated with the corresponding symbols. The colored lines describe the ice/liquid (red) and the liquid/vapor (blue) surfaces enclosing the premelting film. The black arrows show the direction of preferential growth. At the point marked by an asterisk (\(\ast\)), in the region of activated dynamics, growth proceeds by horizontal translation of nucleated terraces. At points such as that marked by a circle (\(\circ\)), above the region of activated dynamics, growth can occur continuously without activation in a steady state of constant film thickness. At points such as that marked by the open triangle (△), above the kinetic coexistence line, droplets can condense and are stabilized transiently with a crater growing inside. At points such as that marked by a square (□), beyond the α → β line, films in the β-thick state can be stabilized transiently and form at the rim of the droplet. At higher pressures, past the kinetic spinodal line (green dotted), such as the point marked with a lozenge (♢), the crystal growth rate can no longer match the condensation rate, and the film thickness diverges. The detailed dynamics corresponding to symbols in the phase diagram is illustrated in Figs. 4 and 5.

• The line of α → β kinetic transition (dotted-blue line in Fig. 3). At sufficiently high saturation, the linear term in ω_k(h) stabilizes the β state transiently, and it is possible to observe the coexistence between α and β states that has been reported in experiments^34,53. The line where the condition is met is obtained by solving a double tangent construct as in usual wetting phase diagrams (Supplementary Note 5).

• The kinetic spinodal line (dotted-green line in Fig. 3), which occurs when Δp_k = − Π_spin, with Π_spin the value at which the interface potential g(h) predicts that the liquid/vapor interface L_lv becomes linearly unstable, i.e., has a spinodal. This condition leads to a line p_spin(T) that can be obtained from Eq. (7), for the choice Π = Π_spin. Crossing this line signals the region of the p-T plane where ice crystal growth cannot match the rate of vapor condensation, and the premelting film thickness diverges.

The slope of the kinetic coexistence lines is dictated by the ratio of k_sl to k_lv, while the separation between kinetic phase lines is dictated by the depth and free energy separation between the minima.

Using gas kinetic theory, crystal growth theory, and literature data for water and ice, we estimate the model parameters k_sl, k_lv, w, η, γ_sl, γ_lv, Δp_sl and Δp_lv for the basal surface of ice (Supplementary Table 3 and Supplementary Note 7). These data, combined with the interface potential g(h) from computer simulations, allows us to draw the kinetic phase diagram depicted in Fig. 3. The shaded area surrounding the sublimation line is the region where crystal growth is a slow activated process, only proceeding via step nucleation and growth. In the absence of any impurities to speed up the nucleation, in this regime the substrate is effectively unreactive for time scales smaller than the inverse nucleation rate, and behaves as dictated by the equilibrium interface potential displayed in Fig. 4a. In practice, the experimental systems reported in ref. ³⁴ contain dislocations, so the crystal freezes by spiral growth and the region of unreactive wetting shown in Fig. 3 for the SG + CW model is not observed. The significance of this change in the growth mechanism can be illustrated by setting w = 0. In this case, the region of activated growth is removed altogether, growth proceeds continuously, and the kinetic phase lines all meet the solid/liquid coexistence line as they approach the triple point (Supplementary Fig. 1). This regime is also relevant for the prism plane above its roughening transition at about 269 K.

Panels a and f show the effective potentials for state points depicted as an asterisk (\(\ast\)) and a circle (\(\circ\)) in Fig. 3, respectively. Panels b–e and g–j show the corresponding solid/liquid and liquid/vapor surfaces at significant milestones in their time evolution (solid red and blue lines, respectively). The dashed red lines indicate the surface location for fully formed ice bilayers, and dashed blue lines show the heights of a premelting film at the α or β minima, as a guide to the eye. Panels a–e illustrate the evolution in the nucleated regime at (p, T) = (455 Pa, 269.5 K) (marked as an asterisk (\(\ast\)) in Fig. 3). Panel a shows the sine-Gordon and interface potential which dictates the surface dynamics. A small terrace nucleated on the solid/liquid surface (panel b) triggers the formation of a similar terrace on the liquid/vapor surface (panel c) and then spreads horizontally (panels d–e). Once the surface has flattened, further growth is not possible until a new terrace is nucleated (Supplementary Movie 1). Panels f–j illustrate the evolution of a droplet quenched to a pressure just below the kinetic liquid–vapor coexistence line at (p, T) = (514 Pa, 269.5 K) (shown as a circle (\(\circ\)) in Fig. 3). The effective free energy, ω_k(h), (panel f) inhibits the growth of liquid wetting films. A droplet (panel g) triggers the formation of a terrace at the rim, which then spreads inside (panel h) and grows to fill the droplet completely (panels i–j). Subsequent growth occurs in a quasi-stationary state of constant film thickness (Supplementary Movie 2).

Interface dynamics

An extensive set of numerical simulations performed over a wide range of the p − T plane and initial conditions confirms that the outcome of the dynamics is in excellent agreement with expectations from the kinetic phase diagram of Fig. 3. Here, we report results performed for the basal surface at T = 269.5 K and varying vapor pressure. Results are reported in reduced units of the model parameters, with \({\kappa }_{1}^{-1}\approx 0.49\) nm for the length scale and \(\tau =3\eta \ {({\kappa }_{1}{\gamma }_{{\rm{lv}}})}^{-1}\approx 0.11\) ns for the time scale.

First, we consider a state very near the sublimation line, where the system is found in the region of activated growth, and water vapor freezes by a growth and spread mechanism. The premelting film here is virtually in equilibrium for the time scale of the simulation, and adopts a thickness of roughly two lattice spacings. In our simulations (Fig. 4b–e and Supplementary Movie 1), an initial terrace mimicking a local defect on the solid/liquid surface L_sl, not observable by optical means, triggers the formation of a corresponding terrace on the liquid/vapor surface L_lv, with a step height equal to the solid lattice spacing. Crystal growth then proceeds by the spreading of the terrace, and the horizontal motion of the solid phase is conveyed to the external liquid/vapor surface. This motion can be observed directly by confocal microscopy, but of course, does not imply the absence of a disordered premelting film (c.f. Fig. 5 in refs. ⁵³ or movie S1 in ³⁴). Once the new full crystal lattice plane is formed, growth becomes stuck again until a new critical nucleus is formed stochastically.

Crossing the line of nucleated growth toward higher saturation, such that ϕ_sl > w, the thermodynamic driving force is large enough to beat the bulk crystal field, and growth then occurs without activation, as in a kinetically rough surface^54,56. However, if ϕ_sl is only marginally larger than w, the process occurs in a stepwise fashion, occurring with large time intervals of no growth, followed by height increments equal to the lattice spacing d_B in a short time²⁶. On further increasing ϕ_sl, crystal growth then proceeds in a truly quasi-stationary manner while the premelting film thickness remains constant, consistent with Eq. (7).

Interestingly, traversing the metastable prolongation of the liquid–vapor coexistence line does not change the growth behavior in any significant way. Although Δp_lv is now positive, Δp_k is still negative, so the thickening of h is still uphill in the effective free energy ω_k(h): i.e., the system behaves as if it is effectively undersaturated with respect to liquid–vapor coexistence and the vapor/liquid interface cannot advance faster than the crystal/liquid interface (c.f. Fig. 4f). For a purely flat interface, the stationary film thickness here is therefore somewhat smaller than that found at the sublimation line, but still remains confined within the α state of the interface potential (see Fig. 4f). A liquid droplet quenched to this region of the kinetic phase diagram is never stable – see Fig. 4g–j and Supplementary Movie 2. Instead, at the contact line of the droplet, terrace formation on the ice is triggered by the action of the disjoining pressure. The crystal then grows and the droplet flattens out, in order to reach a quasi-equilibrium film thickness consistent with Eq. (7). As a transient during the process, the premelting film thickness h can be stable in the β film state, reminiscent of the “sunny side up” states observed in experiment³⁴. Subsequently, the droplet disappears, leaving an Aztec pyramid-shaped solid surface that is covered by an α-thick film. Finally, the inhomogeneity completely disappears, and growth proceeds in a strictly quasi-stationary manner with a flat surface. Notice that during the relaxation process, the droplet is lifted upwards, as a result of the continuous ice growth occurring below. Indeed, comparing Fig. 4g with Fig. 4j, we find that well before the inhomogeneity is washed out, the ice surface grows by about 290 \({\kappa }_{1}^{-1}\), at a rate consistent with Eq. (5). This shows that the relevant relaxation time for large inhomogeneities is far larger than the coarse-graining time scale used to obtain the average growth rate law.

The situation changes significantly when saturation is raised above the kinetic liquid–vapor coexistence line, where Δp_k > 0. For thick enough films, h can now move downhill in the effective surface free energy (Fig. 5a). In this regime, small fluctuations or crystal defects that locally increase the film thickness beyond the spinodal thickness of g(h) trigger the formation of large liquid droplets on top of the premelting film, as observed in experiments—see Fig. 5b–e and Supplementary Movie 3; c.f. Fig. 1-D from ref. ³⁴. Essentially, when Δp_k > 0, the liquid pressure is large enough to sustain the tension of the droplet surface. However, the droplet cannot be fully stable here, since the system is open. The fastest way to decrease the overall free energy while the solid phase grows is to form a large crater below the droplet and then for the two interfaces to separate. Likewise, a droplet quenched to this region behaves initially as described above for droplets below the kinetic liquid–vapor coexistence. The difference is that once a few terraces have been formed at the rim, the crystal grows thereon inside the droplet towards its center by creating a premelting film of α thickness, without the droplet curvature flattening out (Supplementary Fig. 2 and Supplementary Movie 5). As growth proceeds, the interface profiles take a transient shape like that of droplets on soft substrates^62,63, with the solid surface growing higher in the contact line region. A crater develops, but is then filled by the growing solid, before the droplet disappears.

Panels a and f show the effective free energies, ω_k(h) that drive the time evolution of the ice surface at state points depicted as a triangle (△) and a square (□) in Fig. 3, respectively. Panels b–e and g–j show solid/liquid and liquid/vapor surfaces at significant milestones of the dynamics as described in Fig. 4. Panels a–e display the evolution of a surface defect at state point (p, T) = (517.5 Pa, 269.5 K) (shown as a triangle (△) in Fig. 3). The growth of a thick wetting film is now favorable, as illustrated by the negative slope of the effective free energy in panel a. A defect on the solid/liquid surface (panel b) triggers the formation of a liquid droplet (panel c). Ice then grows inside the droplet, forming a large crater (panels d–e) which vanishes eventually when the ice surface catches up with the liquid droplet and attains a stationary premelting layer thickness (Supplementary Movie 3). Panels f–j display the evolution of a droplet at (p, T) = (518.5 Pa, 269.5 K) above the kinetic α → β transition line (shown as a square (□) in Fig. 3). Here, the β state has lower free energy than the α state, as illustrated in panel f. During the time evolution of a droplet (panel g), a thick film of β thickness forms at the rim transiently (panel h), then the droplet evolves a crater (panels i–j) as the ice surface catches up with the droplet (Supplementary Movie 4).

Increasing further the pressure above the kinetic α → β transition line, the free energy of the β film becomes less than that of the α film (Fig. 5f). Therefore, a droplet prepared on top of an α film relaxes to a state where it stands on top of the preferred β state. This corresponds to the “sunny side up” configuration found experimentally at sufficiently high saturation—see Fig. 5g–j and Supplementary Movie 4; c.f. Fig. 1-A from ref. ³⁴. Eventually, the saturation is large enough that the β film metastable minimum is washed away by the linear term Δp_kh in ω_k(h). In this case, the system becomes highly unstable (i.e., linearly unstable to perturbations), and small satellite droplets can form, either in the neighborhood of a larger droplet, or directly from a single local perturbation on the solid surface (Supplementary Fig. 3 and Supplementary Movie 6), a situation that very much resembles experimental observations – see Supplementary Movies S1 and S2 from ref. ³⁴. Eventually, in the long time limit the inhomogeneities disappear completely, and the premelting film thickness diverges. Crystal growth then proceeds below a macroscopically thick wetting film that feeds on the surrounding bulk vapor.