How ice grows from premelting films and water droplets

The study addresses how ice grows and melts near the triple point when its surface is coated by a quasi-liquid premelting film. Classical vapor-deposition models that rely on adatom migration on terraces and kinks neglect this liquid layer and cannot reconcile conflicting measurements of ice growth rates. Experiments and simulations show that the ice/vapor interface becomes increasingly disordered near the triple point, forming a nanometric liquid film whose thickness depends on temperature and vapor pressure. The central question is how the presence and dynamics of this premelting film mediate crystal growth mechanisms and rates under various saturations, including the role of condensed water droplets on the surface. The purpose is to develop a physics-based framework that links equilibrium wetting concepts to out-of-equilibrium crystal growth on a reactive substrate (ice), enabling predictions of growth regimes and rates as functions of ambient conditions, and clarifying how the buried solid dynamics are transmitted to the observable liquid-vapor surface.

Past work on snow crystal growth successfully used phenomenological models to reproduce complex habits, but mapping model parameters (e.g., kinetic anisotropy) to actual temperature and saturation remains unresolved. Measurements show a quasi-liquid premelting layer forms on ice near the triple point, with thickness that increases toward coexistence, supported by various spectroscopic and simulation studies. Classical terrace-ledge-kink models omit premelting films, and attempts to include them have seen limited success. Similar challenges appear in materials systems where a metastable liquid can wet or condense as droplets on a growing substrate, altering growth mechanisms. Wetting theory on inert substrates uses an interface (binding) potential g(h) of adsorbed films to predict equilibrium thickness, but extending this to reactive substrates (where the solid grows or melts) is nontrivial. Prior theoretical tools include mean-field liquid-state theory for short-range structural forces, Dzyaloshinskii–Lifshitz–Pitaevskii theory for long-range van der Waals forces, and coarse-grained surface Hamiltonians (sine-Gordon and capillary-wave) for crystal surfaces. These inform the present model that merges equilibrium interface potentials with nonequilibrium growth kinetics.

The approach combines molecular simulation, equilibrium wetting theory, and mesoscopic thin-film modeling. 1) Interface potential from simulation: Along the ice/vapor sublimation line, NVT simulations using the TIP4P/Ice water model (GROMACS 5.0.5) are performed in the 210–270 K range for an equilibrated ice slab. Film thickness fluctuations h are sampled to obtain the distribution P(h;T), converted to an effective surface free energy per area w(h;T) = −(kT/A) ln P(h;T). The bulk pressure offset Δp_lv(T) between liquid and vapor at solid/vapor chemical potential is calculated by thermodynamic integration. The interface (binding) potential is reconstructed via g(h;T) = w(h;T) + Δp_lv(T) h, and data across temperatures are stitched into a master g(h) over nanometric thicknesses. TIP4P/Ice accurately reproduces solid/liquid and liquid/vapor surface tensions and yields premelting thicknesses of 3–10 Å in the studied range. 2) Theoretical form of g(h): Short-range structural contributions are represented by a damped oscillatory term consistent with mean-field liquid-state theory with renormalization by capillary fluctuations, while long-range dispersion forces are added using a compact form for van der Waals contributions g_vdw(h) that captures the crossover between non-retarded and retarded regimes. The overall g(h) = g_s(h) + g_vdw(h) is fitted to the simulation data and constrained by experimental observations of shallow minima (contact angle ≈2°, energy ≈ −1×10⁻⁵ J m⁻²). The fitted g(h) displays two shallow minima at h_α ≈ 16.0 Å and h_β ≈ 24.5 Å, consistent with a thin-to-thick film transition under increasing supersaturation. 3) Coarse-grained interfacial Hamiltonian: A coupled sine-Gordon plus capillary-wave Hamiltonian describes two interfaces: the solid–liquid height L_sl(x,t) and the liquid–vapor height L_lv(x,t). The functional includes interfacial stiffness/tension terms, a lattice pinning term −u cos(q L_sl), the binding potential g(h) with h=L_lv−L_sl, and bulk terms involving ΔP_sl = P_s(μ) − P_l(μ) and ΔP_lv = P_l(μ) − P_v(μ), with μ set by the vapor reservoir. Parameters can be renormalized to match fluctuation spectra from simulations/experiments. 4) Gradient-dynamics evolution: Near coexistence, dynamics are driven by free-energy gradients. Non-conserved kinetics for crystallization and condensation/evaporation are coupled with conserved hydrodynamics for the premelting film via thin-film (lubrication) approximation. The evolution equations for L_sl and L_lv include kinetic coefficients k_sl and k_lv, viscosity η, and density difference Δρ, plus advective terms for the film. Stochastic effects are implicitly incorporated via the renormalized free energy; equivalently, the equations describe the most probable trajectory for small fluctuations. 5) Coarse-grained rate laws and stationarity: Averaging over one lattice-plane formation time yields rate expressions for continuous growth and introduces the disjoining pressure Π(h)=−∂g/∂h. In quasi-stationary states where the premelting thickness is constant, the stationarity condition is Π(h) = −Δp_k(p_v,T), where the kinetic overpressure Δp_k depends on ambient conditions and kinetic coefficients, reducing to the Derjaguin condition Π(h)=−Δp_v for an inert substrate (k_sl=0). This mapping enables prediction of non-equilibrium wetting behavior using an effective interface potential. 6) Kinetic phase diagram: Using estimated parameters (from gas-kinetic theory, literature surface tensions, viscosities, step free energies) and the fitted g(h), phase lines in p–T space are computed: kinetic coexistence (Δp_k=0), kinetic α→β transition (double-tangent construction on g), and kinetic spinodal (Δp_k=−Π_spin where L_lv becomes linearly unstable). 7) Numerical simulations: The coupled thin-film equations are solved numerically (Matlab) using a method of lines with periodic pseudospectral spatial derivatives and stiff ODE time stepping (ode15s). Nondimensionalization uses length x₁≈0.49 nm and time τ≈0.11 ns. Simulations explore dynamics across the phase diagram for various initial conditions (terraces, droplets), validating analytical predictions. Detailed simulation parameters (e.g., PME grids, thermostats, constraints) and numerical details are provided in the Methods and Supplementary Notes.

• Interface potential: A fitted g(h) consistent with simulations and experiments exhibits shallow α and β minima at h_α ≈ 16.0 Å and h_β ≈ 24.5 Å, with energies near −1×10⁻⁵ J m⁻², matching observed small droplet contact angles (~2°). Short-range structural forces are strongly renormalized by surface fluctuations, and long-range van der Waals forces stabilize finite-thickness minima. - Kinetic regimes and phase lines: The model predicts distinct growth regimes in p–T space separated by kinetic phase lines tied to g(h): (i) Activated terrace nucleation/spreading near the sublimation line where the surface is pinned (φ² < w²); (ii) Continuous (kinetically rough) growth with constant premelting thickness when φ² > w² and Π(h)=−Δp_k has a solution in the α minimum; (iii) Above the kinetic coexistence line (Δp_k>0), transient droplet condensation is possible; (iv) Above the kinetic α→β line, the β film becomes transiently preferred under droplets (“sunny side up” states); (v) Beyond the kinetic spinodal, the film thickness diverges and growth proceeds beneath a macroscopically thick liquid layer. - Surface dynamics: Simulations confirm that, at low saturation, a buried terrace on the solid–liquid interface induces a corresponding terrace at the liquid–vapor surface; growth occurs via terrace spreading while the premelting film remains thin (≈ two lattice spacings). Slightly higher saturation yields stepwise growth with height jumps of one lattice spacing. Below kinetic coexistence, quenched droplets are unstable: they trigger terrace formation at the rim and flatten as the solid fills beneath, leaving a flat α-thick film. Above kinetic coexistence, surface defects or fluctuations nucleate droplets; a crater forms beneath, and the crystal advances mainly under the droplet until a stationary film thickness is reattained. Above the α→β line, droplets relax to stand on a β-thick rim transiently before crater evolution and catch-up. In linearly unstable conditions, multiple satellite droplets can form before long-time divergence of h. - Quantitative observations: The model reproduces transient droplet stabilization well above liquid–vapor coexistence and explains experimental observations of droplet formation and α/β coexistence. In one simulation below kinetic coexistence, while a droplet relaxes, the ice surface rises by ~290 x₁ (≈142 nm) before inhomogeneities vanish, consistent with the coarse-grained growth law. - Mechanistic insight: The hidden solid–liquid dynamics are transmitted to the observable liquid–vapor surface, validating interpretations of terrace translation, spiral growth, and nucleation despite the presence of a nanometric premelting film. The transition from thin to thick film regimes reduces step formation energies under droplets, accelerating growth where droplets condense (e.g., at crystal corners).

The results answer how premelting films and surface droplets mediate ice growth near the triple point by linking nonequilibrium dynamics to an equilibrium-derived binding potential. In steady states near sublimation, the premelting thickness is well defined by Π(h)=−Δp_k, enabling prediction of film thickness versus temperature and pressure for given kinetics. The coupled-interface model shows that the motion of the buried solid is communicated through the premelting film to the outer liquid–vapor interface, reconciling experimental observations of terrace motion with the existence of a nanometric liquid layer. The kinetic phase diagram provides a mechanistic map: thin-film regimes favor activated terrace nucleation and slow growth; crossing kinetic coexistence enables droplet-mediated growth with lower step energy costs, thereby increasing local growth rates; further saturation leads to β-thick films and ultimately instability with macroscopic wetting. These insights rationalize habit changes in snow crystals, explaining why corners (with higher local saturation) preferentially host droplets that can seed faster growth from corners to facet centers and why small crystallites (with higher vapor pressure) exhibit different basal growth mechanisms than larger ones. The framework generalizes to wetting on reactive substrates, suggesting a route to predict growth in other materials where metastable liquids condense during crystallization.

This work establishes a first-principles, mesoscopic framework for ice growth near three-phase coexistence that explicitly incorporates premelting films and droplet condensation. By extracting an interface potential g(h) from simulations and embedding it in a coupled sine-Gordon/capillary-wave Hamiltonian with gradient dynamics, the study identifies distinct kinetic regimes and phase lines that dictate whether growth occurs by terrace spreading under a thin film, by droplet-mediated mechanisms with transient β films, or beneath a macroscopically thick wetting layer. The model explains how buried solid dynamics manifest at the observable liquid–vapor interface and reconciles experimental observations of terrace motion and droplet behavior. These findings bridge microscopic theories and mesoscopic snowflake growth models and are broadly applicable to wetting on reactive substrates. Future work could test the predicted role of droplet condensation in tip splitting using advanced optical microscopy, extend the framework to other facets (e.g., prism plane above roughening), and incorporate impurities or more sophisticated water models to refine quantitative predictions.

• Molecular model limitations: TIP4P/Ice is a non-polarizable, point-charge potential; fully many-body or ab initio models would better capture interfacial polarization but are computationally prohibitive at required scales. - Equilibrium-to-kinetics mapping: The interface potential is inferred along the sublimation line and assumed weakly temperature dependent; extrapolation to broader conditions and construction of a master g(h) introduce approximations. - Coarse-grained dynamics: The model relies on renormalized free energy and thin-film (lubrication) approximations, assumes small deviations from equilibrium, and treats stochastic fluctuations implicitly; true stochastic dynamics and three-dimensional effects are not fully resolved. - Parameter estimation: Kinetic coefficients and some material parameters are estimated from theory/literature and adjusted to match experimental slopes, introducing uncertainty. - Observability: Experimental validation is limited by the difficulty of probing buried solid–liquid dynamics beneath the premelting film; comparisons rely on outer surface observations. - Geometry: Simulations and modeling are effectively two-dimensional with periodic boundaries; real crystal facets, defects, and corner geometries may introduce additional complexities.