Kinetic coefficient for ice–water interface from simulated non-equilibrium relaxation at coexistence

From a macroscopic viewpoint, a freeze-front interface moves at a velocity that balances the rate of latent heat generation and heat transfer. Accordingly, solidification processes are often modeled as ”Stefan problems” with a mathematically sharp moving interface.^1,2 However, at the atomistic scale, the solid–liquid interface has a finite thickness. The temperature may vary across this interface and (for moving interfaces) deviate from the coexistence temperature. For example, the interfacial temperature is slightly supercooled for an advancing freezing front.^3–8 In extreme cases, the ordering of atoms or molecules at the interface may be unable to keep up with heat removal, causing the onset of a glass transition.^9–11 In this work, motivated by the aspects of ice growth inhibition by antifreeze proteins,^12–14 we consider ice–water interfaces near coexistence.¹⁵ 

For a small supercooling ΔT, the growth velocity is proportional to ΔT, i.e., v = ΓΔT, where Γ is the kinetic coefficient. The magnitude and orientational dependence of the kinetic coefficient (Γ) are crucial to understanding growth rates from the melt. There are three prevalent methods¹⁶ to compute Γ: forced velocity simulations, free solidification methods, and fluctuation simulations. The forced velocity method¹⁷ and free solidification method¹⁸ are used to compute Γ from the velocity vs undercooling relationship, while the fluctuation method¹⁹ is used to extract Γ from the time-correlation function of interface fluctuations at coexistence. The latter is a Green–Kubo approach²⁰ based on equilibrium fluctuations. These calculations have primarily been performed for metals^9,11,21–23 and simple liquid models, such as hard spheres²⁴ and Lennard-Jones fluids²⁵ and water.²⁵ This methodology has also been used to study the kinetics of crystal growth in Janus colloids²⁶ and orientational ordering in dipolar particle systems.²⁷ The free solidification method has also been used for studying ice–water interfaces.^10,28,29

In addition to supercooling, interface curvature can also drive freeze fronts to advance or recede.^30,31 In the context of a model for antifreeze protein mechanisms,¹⁴ Sander and Tkachenko noted that the same kinetic coefficient that controls the ice interface velocity under supercooling must also determine the ice interface velocity due to a curvature driving force.¹⁴ Their insightful model, now widely used to interpret data on antifreeze protein efficacy,³² accounts for both supercooling and local curvature effects on ice growth in the moments before and after an antifreeze protein engulfment event, respectively. However, their model uses the kinetic coefficient as an adjustable parameter.

This article presents an eigenfunction-based formulation to estimate Γ at coexistence without applying undercooling. We introduce a sinusoidal perturbation in the ice–water interface and track the evolution of the interface with time. We use the magnitude of this first Fourier mode as a proxy for quantifying interface height. At the start of simulation, the amplitude will be equal to the magnitude of the prepared sinusoidal perturbation. As the interface flattens, the amplitude becomes a small value close to zero. We compute the kinetic coefficient by fitting the time-dependent amplitude to exponential decay. Because we track substantial changes in interface location instead of subtle fluctuations at equilibrium, our method should be relatively insensitive to noise from the order parameter’s structural classifications for individual water molecules. We have computed Γ for two crystal surfaces, namely, basal plane (0001) and primary prism plane (10$1̄$0).³³ The computed kinetic coefficients fulfill two purposes. First, they allow us to use the Sander and Tkachenko model in a predictive, not just interpretive, capacity. Second, they allow us to compare the values of the kinetic coefficient as determined from a curvature-driven velocity to other estimates from a supercooling driven velocity. Accordingly, simulations provide a test of the ansatz that the same kinetic coefficient applies to both situations.

For a single component system, the entropy of fusion (S_L − S_S) and the latent heat Λ are related by⁷ 

where T_m is the melting temperature and subscripts L and S indicate the liquid and solid entropies, respectively. Note that Λ is positive, 333.55 kJ/kg for real water at 273.15 K and ambient pressure. Λ is 299.85 kJ/kg for the TIP4P/Ice model³⁴ at 269.8 K and ambient pressure. Using the definition of free energy of melting

in combination with ΔH_m = Λ, for a small supercooling, the free energy of melting at a temperature of T_m can be written as

Recall that the kinetic coefficient gave the ice growth velocity as

where ΔT = T − T_m represents the undercooling. For a planar interface h$(x⃗)$⁠, the free energy [relative to a system with the ice–water interface at $h(x⃗)=0$] can be written as a functional,

When ΔT < 0 and $h(x⃗)$ moves in the positive direction, the free energy moves in the negative direction as it should for a favorable process. A schematic representation of our sign convention for the interface growth is represented in Fig. 1. This analysis also helps us to choose the correct sign when the melting growth velocity (v) is written in terms of ∂h/∂t,

For a planar interface, the functional derivative $δG/δh(x⃗)=ρiceΛΔT/Tm$⁠. Clearly, Eq. (6) is identical to the expression in Eq. (4), which defines the kinetic coefficient. Now, generalizing Eq. (5) to include effects of surface area changes gives

where γ is the ice–water interfacial free energy. Using the free energy functional that includes surface area changes in Eq. (6) yields the generalized expression for the growth velocity,

Taking the functional derivative³⁵ results in

Equation (9) represents a general starting point for interpreting the results with supercooling and curved solid–liquid interfaces. At equilibrium, Eq. (9) predicts that the radius of curvature for a given supercooling of 1 K is ∼300 Å. Here, we focus on isothermal systems and systems with sufficiently small curvature that $[1+(∇h)2]−12≊1$⁠. In this limit, the interface moves according to the following equation:

To quantify the interface variation with simulation box parameters, we perturb the interface in a known fashion, for example, with a sine wave (or a cosine wave), as shown in Fig. 2.

Assuming the sinusoidal variation for interface height h = A(t) sin(πx/L), the growth velocity of the crystal–melt interface can be written as

On integrating Eq. (11), we have an exponential decay in wave amplitude with time,

Equation (12) predicts the decay time L²ρ_iceΛ/(ΓγT_mπ²) in terms of known parameters L, γ, T_m, ρ_ice, and Λ and unknown parameter Γ. By tracking the decay of amplitude of interface height in molecular dynamics simulations starting from configurations as shown in Fig. 2 and using Eq. (12), we can estimate Γ at the coexistence temperature. The result obtained in Eq. (11) is used to define a Fourier mode component followed by the computation of the correlation function for that Fourier mode. The advantage of using the correlation method does not require special preparation for simulations.

The first step in this analysis is to generate the initial configuration of the ice–water interface exposed to desired crystal planes. Here, we have computed Γ for two principal crystal planes where growth is dominant. We have taken the orthogonal unit cell for hexagonal ice from the work of Hayward and Riemers³⁶ and built an ice crystal by replicating the unit cell in all directions. To generate the ice–water interface, we partition the simulation box into two regions with a sinusoidal boundary, as shown in Fig. 2. We first join the ice slab carved in a sinusoidal fashion with liquid water and then compute the average volume of the simulation box by performing an NpTz simulation (allowing all molecules to move in the simulation). In the computed volume, we place the ice slab on one side and liquid water on another side at a little higher density to account for the vacuum gap between ice and water. Then, we compute the bulk density of water by doing NVT simulation only allowing water molecules to move. Because of the interfacial effects, it may be necessary to add/remove a few water molecules from the bulk to match the equilibrium density corresponding to 269.8 K. The box dimension of the system when the interface is exposed to the primary prism plane is $78.68×44.34×143.40$ ų and to the basal plane is $81.72×44.34×139.10$ ų.

We have used the TIP4P/Ice force field³⁴ to simulate water and ice molecules. Recent studies with the TIP4P/Ice model predicted the coexistence temperature to be 269.8 K at ambient pressure.^29,37,38 The Lennard-Jones interactions were truncated at 8.5 Å, and standard long-range tail corrections were included.³⁴ Similarly, short-range columbic interactions were truncated at 8.5 Å and PPPM³⁹ is used to compute long-range interactions. All simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) molecular dynamics package.⁴⁰ The production runs were carried out in the canonical ensemble at the coexistence temperature. A Bussi thermostat⁴¹ is used to maintain the temperature during the run with a damping factor of 2 ps. The details about the choice of the damping factor are presented in the supplementary material (see Sec. S.I. for details). A time step of 0.5 fs is used for integration, and the shake algorithm⁴² is used to maintain the rigidity of the water molecule. Configurations were saved at 1.0 ps intervals to analyze the interface evolution with time. The total length of simulation for most relaxation trajectories in this study is 15 ns. Simulations were slightly longer in analysis of the system size effects.

Computing the precise location of the interface is a crucial task. Defining an order parameter that can clearly distinguish the ice and liquid-like molecules will aid in the computation of interface location. We use the rotationally invariant, neighbor-averaged bond-order parameter $q̄6$ for this purpose.⁴³ The probability distributions of $q̄6$ for bulk structures of ice and water are shown in Fig. 3. A threshold of 0.35 is used to distinguish ice molecules from water molecules. Two-dimensional histograms of order parameters were made in x and z by averaging over the y-direction. For each x, averaged $q̄6$ values at each z are used to identify the precise interface location. A typical variation of the $q̄6$ vs z at a series of x values is shown in Fig. 4. The interface is not sharp during simulation; the interface is diffusive, indicating interfacial thickness. Once we have quantified the interface location for all bins in x, we compute the magnitude of the first Fourier mode as a measure of relaxation of the ice–water interface over the trajectory.

We computed the kinetic coefficient for two crystal planes of hexagonal ice, namely, the basal plane and primary prism plane. As expected from the theory, the interface relaxation profiles follow an exponential decay with time. Figure 5 represents the exponential decay for the primary prism plane. We fitted the decay to exponential functions, as shown in Eq. (12). Kinetic coefficients were estimated from the fitted exponents. Interfacial free energy values for the ice–water system are needed to compute the kinetic coefficient from the fitted exponent. We have taken γ from the earlier study by Espinosa et al.⁴⁴ The numerical values of γ according to the TIP4P/Ice forcefield are 27.2 mN/m for the basal plane and 31.6 mN/m for the primary prism plane. We computed Γ values using these data and simulation parameters. Computed kinetic coefficients for different hexagonal ice planes are presented in Table I. The error bars are computed from five independent estimates by starting simulation with different initial velocity distributions at 269.8 K. To check for the finite size effects, we have varied the wavelength of the perturbed sine wave by increasing the box length (L). We have tested for finite size effects with three different box sizes in the x-direction, i.e., $L1=62.948$⁠, $L2=78.685$⁠, and $L3=94.422$ Å. The computations of error bars and system size effects are reported in the supplementary material (see Sec. S.II. for details). The computed kinetic coefficients with different system sizes are presented in Table II.

The computed Γ values are consistent within the error bars for different system sizes. Hence, the system sizes considered here are sufficient to eliminate finite-size effects. Qualitatively, the predictions are consistent with both simulations and experimental data. The basal plane will have a slower growth due to a faceted growth at low degree of supercooling. The prism face will have a higher growth because of the rough nature of the surface. The absolute values for Γ in this work are about 30% over-estimated compared to previously reported simulation and experimental values in the literature. A comparison of Γ estimations is presented in Table I. The Γ data from other studies have been computed using the linear relationship between the growth rate and degree of under cooling. The discrepancies may arise from sources of error in this work or earlier studies. In this work, the initial sinusoidal interface preparation may affect the results. In the free solidification studies, there are difficulties maintaining a well-controlled supercooled temperature at the moving interface. In calculations based on equilibrium interface fluctuations, such as attachment frequency estimates for nucleation studies,^49–51 molecular-scale fluctuations in the order parameters used to identify interface locations may be more rapid than actual fluctuations in the interface location. In addition, note that some of the studies mentioned above were carried out with the TIP4P/2005 water model,^10,25 so the results are not expected to precisely agree. The method developed in this work is most appropriate for studying growth kinetics on rough surfaces. The prismatic face is rough, but the basal face becomes faceted at low supercooling. For faceted faces, growth will proceed via spiral growth mechanisms.^52–54 It is not clear whether spiral growth (which depends on factors such as kink density along edges) can be described by a kinetic coefficient, but certainly, the kinetic coefficient cannot describe situations where there is a change in the growth mechanism and, therefore, in the rate at some critical supercooling.^55,56 Crucially, we emphasize that our simulations do not contain any screw dislocations, and therefore, the results cannot correctly describe a spiral growth regime. Accordingly, we cannot be confident in the kinetic coefficient for the basal face even though the numerical errors narrowly bracket its value.

There are several contexts in which the kinetic coefficient will be used. In nucleation theory, the prefactor (which involves a kinetic attachment frequency) should be proportional to Γ. Suppose that we can compute the value of Γ for both the coarse-grained and atomistic models. It should be possible to modify predicted nucleation rates from coarse-grained water models by the ratio of Γ values to obtain an estimate that more accurately reflects interfacial attachment/structuring kinetics in the all-atom model. This scaling trick should help us to get accurate nucleation rates from coarse-grained models, such as mW that get many thermodynamic properties correct,^48,57 but tend to give overly fast dynamics. Similarly, for predictions of growth rates, especially fast growth in the first moments after nucleation, the kinetic coefficient should help predict absolute growth velocities and evaluate criteria related to fingering/dendrite instabilities.⁵⁸ Finally, we plan to use the kinetic coefficient in multistep, multiscale models of growth inhibition by antifreeze proteins.

The kinetic coefficient relates the velocity of a solidification front to the supercooling at the solid–liquid interface. This article developed a general expression, giving the solidification velocity in terms of the functional derivative of free energy with both supercooling and curvature driving forces. We demonstrated a simple method to compute the kinetic coefficient in non-equilibrium simulations at the coexistence temperature. Each simulation is a relaxation trajectory starting from a prepared and annealed two-phase system with a sinusoidal ice–water interface. Theoretical analysis showed that the amplitude of the sinusoidal perturbation should decay exponentially. This was confirmed in simulations, and the results were fit to the theoretical model, providing the kinetic coefficient. Estimated values in Γ for the Ice–Ih–water interface agreed qualitatively with previous data in the literature from simulations and experiments. However, the absolute values were larger by about 30% compared to earlier studies.^10,25,28,29 We discussed potential sources of the discrepancy in our method and in earlier methods. We also noted that our results are from TIP4P/Ice simulations, while others are from different variants of the TIP4P model of water. In comparison to other approaches, the Green–Kubo methods based on equilibrium time correlation functions had the advantage of not requiring a prepared non-equilibrium initial condition. In contrast, the non-equilibrium relaxation should provide a larger signal-to-noise ratio for extracting the kinetic coefficient. We discussed applications, especially those for which curvature effects are important, including nucleation, fingering instabilities, and engulfment by an advancing interface to form inclusions as occurs when antifreeze proteins are engulfed.⁵⁹ 

See the supplementary material for the choice of the thermostat damping factor and the kinetic coefficient data obtained for different ice-Ih crystal planes studied in this work.

We thank Francesco Paesani, Valeria Molinero, Ido Braslavsky, Jon Higdon, and Brian Laird for discussions. This work was supported by the Air Force Office of Scientific Research through MURI Award Grant No. FA9550-20-1-0351.

The authors have no conflicts to disclose.

Ravi Kumar Reddy Addula: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Baron Peters: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Supervision (equal); Writing – original draft (supporting); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.