How ice grows from premelting films and water droplets

Close to the triple point, the surface of ice is covered by a thin liquid layer (so-called quasi-liquid layer) which crucially impacts growth and melting rates. Experimental probes cannot observe the growth processes below this layer, and classical models of growth by vapor deposition do not account for the formation of premelting films. Here, we develop a mesoscopic model of liquid-film mediated ice growth, and identify the various resulting growth regimes. At low saturation, freezing proceeds by terrace spreading, but the motion of the buried solid is conveyed through the liquid to the outer liquid–vapor interface. At higher saturations water droplets condense, a large crater forms below, and freezing proceeds undetectably beneath the droplet. Our approach is a general framework that naturally models freezing close to three phase coexistence and provides a first principle theory of ice growth and melting which may prove useful in the geosciences. Melting and crystallization of ice close to the triple point is mediated by a thin liquid layer of water. Within an extended theoretical-numerical model, Sibley et al. capture both equilibrium properties and kinetic aspects of the film evolution to arrive at a broader perspective of how ice grows.

T he growth and melting of ice play a crucial role in numerous processes, from the precipitation of snowflakes 1 , to glacier dynamics 2 , scavenging of atmospheric gases 3 or climate change 4 . Yet, despite ice ubiquity both in large masses on the poles and as tiny crystals in the atmosphere, we still do not fully understand how ice actually grows (or melts) [5][6][7][8] .
Conflicting experimental measurements of ice growth rates [9][10][11][12][13] have been analyzed under a framework of classical crystal growth based on direct deposition from the vapor phase, followed by the subsequent two-dimensional migration of adatoms onto surface kinks 14 . However, the last two decades have witnessed great progress in the experimental characterization of the ice/vapor interface at equilibrium 7 . Results from different experimental techniques [15][16][17][18] , as well as computer simulations, confirm that the surface disorder of ice grows steadily as the triple point is approached, and what is sometimes referred to as a "quasi-liquid layer" of premelted ice is formed on its surface [19][20][21][22][23] . Unfortunately, classical growth models based on the terrace-ledge scenario do not account for the impact of premelting films at all, and attempts to incorporate this effect have met only limited success [24][25][26] .
Our current understanding of snow crystal growth illustrates this uncomfortable situation. The primary habit or aspect ratio of these familiar hexagonal crystallites can change dramatically with small changes in temperature and saturation, from extremely elongated needle-like crystals to almost flat plate-like dendrites 27 . But despite their variety and complexity, these shapes can be described using phenomenological models with amazing accuracy, based on just a number of parameters 28,29 . Particularly, the primary habit is dictated by a kinetic growth anisotropy factor, describing the ratio of horizontal to vertical growth rates 29 . Unfortunately, the mapping of this phenomenological parameter to the actual ambient conditions in the atmosphere, namely, temperature and water saturation, remains a long-standing topic in crystal growth science 12,13,25 . Accounting explicitly for the premelting layer appears an essential requisite to unveil the dependence of growth rates on ambient conditions. The difficulty to incorporate the role of premelting films on crystal growth theories is also encountered in many systems of interest in materials science [30][31][32] , where the partially stable liquid phase can even condense into liquid droplets on the growing substrate 15,[33][34][35] and change the mechanism of crystal growth substantially.
The problem is akin to one encountered in the theory of wetting, where one studies how a metastable liquid phase (say, water), adsorbs at the interface between a solid substrate (ice) in contact with a vapor (water vapor) as the liquid/vapor coexistence line is traversed 36 . For an inert substrate, wetting is very well understood in terms of the underlying interface potential g(h) that measures the free energy of the adsorbed film as a function of film thickness h 37 . Out of equilibrium, however, the substrate continually feeds from the adsorbed film at the expense of the mother phase, so it is debatable whether it is possible to define meaningfully a film thickness and corresponding interface potential.
Here, we combine state-of-the-art computer simulations, equilibrium wetting theory, and thin-film modeling to describe the kinetics of the ice surface in the vicinity of the triple point within a general framework for wetting on reactive substrates. Our results show that as the vapor saturation increases, ice first grows by terrace spreading below a premelting film with a well-defined stationary thickness. At higher saturations, however, the premelting layer thickness diverges, and growth actually proceeds from below a bulk water phase. In between these two regimes, at intermediate saturations, droplets condense on the ice surface, and growth proceeds mainly under the droplets. The different regimes are separated by well defined kinetic phase lines, whose location can be mapped to an underlying equilibrium interface potential.

Results
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 25 . 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 38 . 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 ωðh;TÞ ¼ Àk B T 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. 41 ). 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 manybody forces differ greatly depending on the local position. Such a complicated situation is best described with electronic quantummechanical 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 44 . Although this is a point-charge non-polarizable potential, it predicts accurately both the solid/ liquid and liquid/vapor surface tensions 45 . 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, κ 1 and κ 2 are inverse decay lengths (whichever is shorter is the inverse bulk correlation length), and q 0 % 2π d À1 0 , where d 0 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. 41 ).
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 51 parameterized the dielectric response of ice and water to numerically calculate these contributions with Dzyaloshinskii-Lifshitz-Pitaivesky theory. Following ref. 52 , 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 highfrequency 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 γ lv ðcos θ À 1Þ $ À10 À5 J m −2 .
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 0 , 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 41 , 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 −5 J m −2 , as observed in experiment 34 . 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π d À1 B , 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 20 . 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 58 . 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, 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 61 . 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 32 . When the fluctuations are small, the coarsegrained 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 meanfield 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): . 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 ϕ 2 sl < 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 ϕ 2 sl ¼ 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 ϕ 2 sl > 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 〈∂ t L sl 〉 from that of 〈∂ t L 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 h ∂h ∂t i ¼ 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 k h 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 .
• 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. 34 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.
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 κ À1 1 % 0:49 nm for the length scale and τ ¼ 3η ðκ 1 γ lv Þ À1 % 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 NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20318-6 ARTICLE NATURE COMMUNICATIONS | (2021) 12:239 | https://doi.org/10.1038/s41467-020-20318-6 | www.nature.com/naturecommunications 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. 53 or movie S1 in 34 ). 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 26 . 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 stablesee 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 quasiequilibrium 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 34 . 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 κ À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. 34 . 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.
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. 34 . Eventually, the saturation is large enough that the β film metastable minimum is washed away by the linear term Δp k h 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 observationssee Supplementary Movies S1 and S2 from ref. 34 . 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.

Discussion
In our study, we have discussed ice premelting, but our results rationalize the behavior of out-of-equilibrium premelting films and wetting on reactive substrates quite generally. In particular, we see that for small deviations away from the sublimation line, freezing occurs in a steady state regime with constant film thickness. In this regime, the thickness of the premelting film is dictated by an equilibrium interface potential and the underlying growth mechanism. For a given growth mechanism, our results show that the outcome of the out of equilibrium dynamics may be predicted accurately from an underlying free energy functional in analogy with wetting on inert substrates. As long as the system remains in this steady-state, the premelting film thickness is well defined and depends both on temperature and pressure. We emphasize that it is not possible to interpret the dynamics of the quasi-liquid layer without taking into account the behavior of the underlying substrate. In particular, our results demonstrate that the complex dynamics of a buried solid surface can be conveyed to the experimentally accessible outer surface of the quasi-liquid film. We also confirm that observation of terrace translation, spiral growth, and nucleation observed in experiments is fully consistent with the existence of a nanometer-thick premelting film as observed in simulations 19,20,22,23 . Accordingly, the motion of the experimentally accessible outer surface may be used to interpret the hidden dynamics of the inner surface, very much in agreement with expectations of the Kuroda-Lacmann model 25 .
The change from a thin to a thick film regime that occurs across well-defined kinetic lines can result in a significant change in the mechanism for crystal growth. In the thin-film regime, the growth of steps is energetically expensive, because the nuclei are barely buried by the premelting film: steps formed feel a large inhomogeneity as the density changes from solid to vapor across a thin water film. As the kinetic coexistence line is traversed, however, liquid droplets condense on the ice surface. Steps formed below feel a much smaller tension, similar to that at the ice/water interface. Their free energy of formation is therefore much smaller, and leads to a significant increase in the growth rate at places where droplets have condensed. This has immediate implications for our understanding of ice crystal growth 12,29 . Since crystal corners have high local saturation, droplets are more likely to condense there, providing a source of water for the crystal to feed by growth and spread mechanism from corners towards facet centers as observed in experiments 9,12,64 . Furthermore, small crystallites with large vapor pressure are more likely to have droplets condense at their corners, explaining why the growth mechanism on a basal facet appears to be different in large and in small crystallites 64 . Interestingly, this suggests that droplet condensation could play a role in the tip splitting mechanism of ice grown from the vapor. Advanced optical microscopy appears a candidate technique for the verification of this hypothesis. In summary, we find a discontinuous change of crystal growth mechanisms with saturation. Combined with recent findings of non-monotonic temperature dependence of step-free energies 23,41 , our results could help fill the gap between microscopic theories and mesoscopic models of snowflake growth 29 .

Methods
Computer simulations. Simulations of an equilibrated ice slab in the NVT ensemble are performed in the temperature range 210-270 K for the TIP4P/Ice model 44 using GROMACS 5.0.5. The equations of motion are integrated using the Leap-Frog algorithm, with a time step of 3 fs. Bond and angle constraints are applied using the LINCS algorithm. The canonical ensemble is sampled using thermostated dynamics with the velocity rescale algorithm 65 . The Lennard-Jones interactions are truncated at a distance of 9 Å. Electrostatic interactions are evaluated using the Particle Mesh Ewald algorithm with the same real space cutoff. We calculate the reciprocal space term using a total of 80 × 64 × 160 vectors in the x, y, z reciprocal directions, respectively. We use a 0.1-nm grid spacing and fourth-order interpolation scheme for the charge structure factor. Simulations are carried out in systems consisting of 8 × 8 × 5 unit cells of pseudo-orthorhombic geometry, each containing 16 molecules. The initial configurations for the solid ice slab are prepared with a random realization of the hydrogen bond network, following ref. 66 . One such initial lattice is provided as Supplementary Data 1. This is then simulated at 1 bar to obtain the equilibrium lattice parameters and placed in vacuum for further equilibration in the NVT ensemble during 15 ns. Averages are collected on production runs 35-ns long. During the simulations, we identify structurally liquidlike molecules using the q 6 order parameter 67 . Once these molecules are identified, we determine the locations of the liquid-vapor and solid-liquid surfaces as explained in ref. 45 . From these two surfaces, we calculate the local film thickness as the difference between these, h(x) = L lv (x) − L sl (x). For the calculation of the interface potential, the local film thickness for a given configuration is laterally averaged, in order to obtain the average liquid film thickness. The set of global film thicknesses obtained are used to compute the probability histograms P(h), from which g(h) can be calculated as detailed in the Supplementary Note 1. The results for g(h) are fitted to the model described in the main text. Parameter values and further details are given in Supplementary Table 1 and Supplementary Note 3.
Gradient dynamics. Numerical computations of the dynamics of the thin-film equations are performed using the method of lines, similar to that used in ref. 68 , but with a periodic pseudospectral method for the spatial derivatives. The method is extended to evolve the two interfaces (solid-liquid, and liquid-vapor), with coupling terms involving mass transfer and the two interface potentials naturally included. For the evolution of the solid-liquid interface, a pinning effect in the horizontal direction can occur if too few mesh points are used. Consequently, rather than using an extremely large number of points in the finite difference scheme, we implement a periodic pseudospectral method which significantly increases the rate of numerical convergence. The numerical method uses discretization on a regular (periodic) grid and a band-limited interpolant derived using the discrete Fourier transform and its inverse to form the differentiation matrices which act in real space. The presence of the premelting film avoids the need to explicitly evolve the contact lines, in comparison to some of our previous work using pseudospectral discretisation 69,70 . For the time stepping, the ode15s Matlab variable-step, variable-order solver is used. Our numerical calculations are performed on the non-dimensionalised version of the model equations. We find that choosing κ À1 1 % 0:49 nm and τ ¼ 3η ðκ 1 γ lv Þ À1 % 0:11 ns as our units of length and time in the non-dimensionalisation works well. Further details of the method and initial conditions are given in Supplementary Note 8 and Supplementary References.
Model parameters. Phase coexistence data required to compute Δp sl , Δp lv , structural properties of ice, and surface tension coefficients are obtained from the literature as described in Supplementary Tables 2-3 and Supplementary Note 6. The kinetic growth coefficients k sl is estimated from the kinetic theory of gases, and k lv is chosen such that the kinetic coexistence line has a slope similar to experiments. The sine-Gordon coefficient u = 1.3 × 10 −4 J m −2 is chosen to match stepfree energies from the literature. The viscosity is taken from literature values of undercooled water. Further details of the choice of model parameters are given in Supplementary Note 7. The actual model parameters used in this work may be found in Supplementary Tables 1-3.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request. Source data are provided with this paper.