2.1 Ice nucleation rate

To observe ice nucleation, we equilibrate the water at a supercooled temperature and then cool the water at a constant rate of 0.25 K ns⁻¹ until ice forms. The same cooling rate is used for all simulations in this study. All simulations employ periodic boundary conditions along the three spatial dimensions and use a simulation time step of 5 fs. The starting temperature for the cooling ramp at all pressures is 240 K for ML-mW and 225 K for mW, corresponding to 52 and 48 K of supercooling, respectively, relative to their melting temperatures at 1 atm. Ice is identified using the q₆ order parameter, where clusters of molecules with q₆>0.54 are considered to be ice (Steinhardt et al., 1983; Lupi et al., 2014; Rosky et al., 2022). Using the same method as Rosky et al. (2022), the freezing temperature of each cooling ramp is identified by the sharp increase in the ice-to-liquid ratio. The steepest point on a sigmoidal fit to the ice/liquid ratio curve is used to select the nominal nucleation temperature.

Intensive nucleation rates are found by running the same cooling simulation a minimum of 20 times and recording the freezing temperature of each cooling run. These repeated cooling runs form a statistical distribution of freezing temperatures that can be divided into temperature bins with a corresponding intensive nucleation rate, j_het(T), within each bin. The calculation of intensive nucleation rate values for this process is described in Rosky et al. (2022) using the methods of Zobrist et al. (2007). The widths of the temperature bins are displayed as uncertainty bounds on our data points. As determined from Poisson statistics, we report with 99 % certainty that the intensive nucleation rates shown are contained within the bounds of the temperature bin. We reference Table 2 of Koop et al. (1997) to obtain the 99 % upper and lower confidence bounds for the number of freezing events observed in each temperature bin. The cooling rate and substrate surface area used in our study allow us to sample intensive heterogeneous-nucleation rates on the order of magnitude j_het=10²⁴ m⁻² s⁻¹.

2.2 Ice-nucleating substrate immersed in water

We simulate heterogeneous ice nucleation by inserting an ice-nucleating substrate into a box of water with periodic boundary conditions. This configuration is shown in Fig. 1a. The Fig. 1a simulation box has dimensions of 49.22 × 49.02 × ∼ 58.75 Å, containing 4906 water molecules with a 10 Å thick sheet of ice-nucleating substrate inserted at the center of the z axis. The total surface area of substrate in contact with water (including both the top and bottom of the substrate) is 48.25 nm². The substrate molecules are held fixed with zero velocity. The x and y dimensions of the simulation cell are kept constant, while the z-axis dimension may be adjusted to accommodate changes in pressure. After equilibration of the simulation box, we use cooling ramps to measure j_het(T) at three pressures: 1, −500, and −1000 atm.

Figure 1Heterogeneous ice nucleation simulation configurations with periodic boundary conditions employed along all three dimensions. Water molecules are cyan, and the substrate material is brown. All configurations have the same total surface area of contact between water and substrate to within 6 %. (a) Unconfined water with a barostat. (b) Confined water with a barostat, designed to capture the effect of confining water along the z axis between substrate surfaces. (c) Capillary water bridge, used without a barostat so that negative Laplace pressure may arise naturally within the water capillary bridge due to the curved geometry of the air–water interface (gray shading). (d) Narrow capillary bridge, which simulates confinement of the water along the x axis between the two air–water interfaces (gray shading). In configurations (b), (c), and (d), the periodic boundary conditions cause the bottom surface of the substrate to be in contact with the water molecules at the top of the simulation cell, thus forming the water bridge in (c) and (d).

Download

Simulations containing 4906 molecules (a volume of roughly 125 nm³) are a conventional choice for studies that aim to simulate bulk water properties (e.g., Molinero and Moore, 2009; Lupi et al., 2014; Li et al., 2011). The interaction forces felt between two water molecules in our simulation go to zero when molecules are spaced further than ∼ 4 Å apart, and forces between water molecules and substrate molecules go to zero for separations beyond ∼ 5 Å (Molinero and Moore, 2009). Thus, the chosen thickness of our substrate ensures that water molecules on either side of the substrate will not interact with each other. We refer to the configuration shown in Fig. 1a as being unconfined because the vertical dimension is large enough that much of the water is unaware of the substrate. Structural correlations in the liquid decay to zero beyond 1.5 nm (Cox et al., 2015; Bi et al., 2016; Lupi et al., 2014); the height of the unconfined simulation cell is more than twice this length. To confirm that the nucleation rate on the substrate is not influenced by confinement effects in the unconfined cell, simulations at 1 atm were performed with the z axis dimension being doubled to 100 Å, showing no change in the heterogeneous-nucleation rate.

2.3 Confinement

To simulate the effects of confining the water between the two substrate surfaces, we repeat the simulations with reduced box heights, as shown in Fig. 1b. We test z-axis heights of h=30, 24, and 18 Å, which will later be used as the heights of our water capillary bridges. For both the unconfined and confined configurations, equilibration is carried out with a Berendsen barostat (damping constant of 10 ps) and a Bussi thermostat (damping constant of 5 ps). The cooling runs are carried out in the isothermal–isobaric (NPT) ensemble, implemented by employing an isobaric–isoenthalpic (NPH) ensemble using a Nose–Hoover barostat (damping constant of 10 ps) coupled with the Bussi thermostat (Allen and Tildesley, 2017). Only the ML-mW model is used for the confinement simulations.

2.4 Water capillary bridges

We have taken steps towards addressing whether these magnitudes of negative pressure within water can be found in nature by simulating water capillary bridges. A volume of water placed between two hydrophilic substrate surfaces forms a capillary bridge that is expected to have negative Laplace pressure within the water. As shown in Fig. 1c, we use the same ice-nucleating substrate used in the previously described configurations to construct water capillary bridges with heights h=30, 24, and 18 Å. We apply the same constant-cooling simulation procedure to obtain intensive heterogeneous ice nucleation rates within the capillary bridges. In these simulations, we remove the external barostat so that negative pressure within the water is solely a result of the capillary bridge geometry. Only the ML-mW model is used for the capillary bridge simulations.

Water capillary bridge simulations are equilibrated and carried out in the NVT ensemble (number of molecules, volume, and temperature is conserved during each time step) using a Bussi thermostat with a damping constant of 5 ps. The simulation cell dimensions are all held fixed during the capillary bridge simulations. The widths of these capillary bridges are measured at the water–substrate interface. The total area of contact between the water and substrate remains consistent (∼ 48.25 nm²) to within 6 % for all capillary bridge configurations. The desired contact area is achieved by tuning the y-axis dimension of the simulation cells to adjust the length of the capillaries along the y axis. Intensive nucleation rates are ultimately obtained by dividing by the water–substrate surface area (m⁻² s⁻¹), and maintaining a consistent surface area ensures that all simulations sample the same magnitude of nucleation rate.

We look more closely at the ice nucleation events within the water capillary bridges by identifying the location of all nucleation events. The identification of ice nucleation location is done visually using the clustering capability of the Open Visualization Tool (OVITO), a visualization and analysis tool for molecular-dynamics simulation data (Stukowski, 2010). The center of mass of the initial ice cluster is used as the freezing location. The initial ice clusters contain an average of 25 water molecules, and their positions have an uncertainty of 5 Å along each axis.

Throughout this paper, data points will be presented using circles to indicate unconfined heterogeneous freezing in the configuration shown in Fig. 1a. Squares will indicate confined heterogeneous freezing, as in the example of Fig. 1b. Diamonds will indicate water capillary bridge data, as in Fig. 1c. Narrow diamonds will represent the narrow capillary configuration, as shown in Fig. 1d.

2.5 Interaction potential between water and substrate

In molecular-dynamics simulations, the forces between molecules are defined by interaction potentials. Lupi et al. (2014) introduced an interaction potential to model mW water with a carbon substrate, with the contact angle between mW water and carbon tuned to 86^∘. We use this same interaction potential to simulate heterogeneous freezing of the mW water model. Although our focus in this work is not to model any particular substrate, this mW–carbon interaction potential has the benefit of having already been used in studies of ice nucleation (Lupi et al., 2014; Bi et al., 2016). A description of the equations of interaction for this potential and its parameters is included in Appendix A.

We observe that a simulation box containing 4096 mW water molecules at 1 atm with no substrate freezes homogeneously at an average temperature of 202 K when cooled at a constant rate of 0.25 K ns⁻¹. After inserting a layer of substrate into the center of this same volume of water (as in Fig. 1a) and using the same cooling rate, we observe that heterogeneous freezing on the substrate takes place at an average temperature of 217.5 K, a 15.5 K increase over the homogeneous-freezing temperature. This is consistent with the results of Lupi et al. (2014), who reported a 12 ± 3 K difference between heterogeneous and homogeneous-freezing temperatures when cooling at a rate of 1 K ns⁻¹.

Our previous study (Rosky et al., 2022) indicates that the ML-mW water model is more appropriate for studying pressure effects on ice nucleation because it exhibits a molar-volume difference between liquid water and ice that is closer to real water (Chan et al., 2019). Therefore, in order to increase the ability to translate our findings onto real water droplets, we focus our study on the ML-mW model instead of mW. We define a substrate interaction potential for the ML-mW model that is a modified version of the mW–carbon interaction potential defined by Lupi et al. (2014). We adjust the parameters of the interaction potential to form a substrate on which ML-mW water is stable within the range of negative pressures that we study and exhibits a similar temperature difference between heterogeneous and homogeneous freezing as that seen in the mW–carbon interaction. A simulation box containing 4096 ML-mW water molecules at 1 atm freezes homogeneously at an average temperature of 214.5 K when cooled at a constant rate of 0.25 K ns⁻¹. With the substrate inserted, the average heterogeneous-freezing temperature on this substrate is 228.5 K, a change in temperature of 14 K.

The resulting ML-mW–substrate interaction potential has a contact angle between water and substrate that is smaller than that of the mW–carbon potential. This smaller contact angle is consistent with the studies of Bi et al. (2016), Cox et al. (2015), and Lupi and Molinero (2014), in which the interaction potential is adjusted in a similar manner as here to modify the substrate hydrophilicity. Our estimate of the ML-mW–substrate contact angle is 50.5^∘ and will be discussed in Sect. 3.2. The resulting interaction parameters used in this work are summarized in Table A1.

3.1 Heterogeneous-nucleation rate with negative pressure

Freezing of water on an ice-nucleating substrate (Fig. 1a) is simulated at 1, −500, and −1000 atm to identify how the intensive heterogeneous-nucleation rate j_het behaves at negative pressures. The cooling rate and the area of the substrate in contact with water are kept fixed as we change the pressure of the system so that the same magnitudes of j_het are sampled at all pressures. For each pressure setting, we identify the temperature at which the intensive nucleation rate j_het is equal to 10²⁴ m⁻²s⁻¹, thus obtaining contours of constant j_het in pressure–temperature coordinates. Figure 2 shows these results for both the ML-mW model and the mW model. Intensive homogeneous-nucleation rate data j_hom (m⁻³ s⁻¹), as well as equilibrium melting points T_m, are also included in these plots for comparison with Rosky et al. (2022).

Figure 2Melting (open green circles), heterogeneous-freezing (closed red circles) and homogeneous-freezing temperatures (open blue circles) versus pressure. The heterogeneous-freezing results are new, using an ice-nucleating substrate inserted into a filled box of water with periodic boundary conditions (See Fig. 1a). The melting and homogeneous-freezing results are reproduced from Rosky et al. (2022) for comparison. The red and blue shading represents the 99 % confidence intervals (uncertainty) for the simulation data. Dashed lines use the slope predicted by Eq. (1) to obtain a best fit to the intercept. Solid lines are a best linear fit of both slope and intercept. The numeric values of these slopes are given in Table 1. Contours of constant heterogeneous-nucleation rate are, to within sampling uncertainty, linear and nearly parallel with lines of constant homogeneous-nucleation rate for both the ML-mW model (a) and mW model (b). The gray melting point line from Marcolli (2017) is a fit to experimental measurements, which the ML-mW model reproduces more realistically. Note that −100 MPa = −1000 atm.

Download

Comparing the data for intensive heterogeneous-nucleation rates (j_het=10²⁴ m⁻²s⁻¹) with the data for intensive homogeneous-nucleation rates (j_hom=10³² m⁻³s⁻¹), we see that the two follow a similar slope in pressure–temperature coordinates. Most significantly, we observe that the increase in temperature as a function of pressure for j_het can be approximated as linear within the sampling uncertainty, indicating that the use of a linear estimate for $(\mathrm{\Delta }T/\mathrm{\Delta }P{)}_{\mathrm{het}}$ may be applied to heterogeneous ice nucleation. For the mW model in particular, the slope predicted by Eq. (1) fits exceptionally well to both the homogeneous- and heterogeneous-freezing data. Meanwhile, Eq. (1) seems to underestimate the heterogeneous slope of the ML-mW model. The source of excess steepness in the ML-mW heterogeneous slope is not uncovered in this study and merits further investigation. Despite this weakness, these results do indicate that the slope predicted by Eq. (1) may still be applicable to heterogeneous nucleation to within the simulation uncertainty. The similarity between the heterogeneous and homogeneous slopes for each respective water model suggests that l_f and Δv_ls remain key factors in determining the pressure dependence of the heterogeneous-nucleation rate. Indeed, the consistent change in slope between the mW and ML-mW models indicates that the magnitude of density difference between liquid water and ice plays a similar role for heterogeneous freezing as that which was previously found for homogeneous freezing (Rosky et al., 2022).

In classical nucleation theory, j_het has the following form (Lamb and Verlinde, 2011):

where f_het is a heterogeneous compatibility function, typically related to the contact angle of water on the substrate, and Δμ is the chemical potential difference for the phase change. The factor $C=\mathrm{16}\mathit{\pi }{\mathit{\gamma }}_{\mathrm{ls}}^{\mathrm{3}}/\left(\mathrm{3}{k}_{\mathrm{B}}{\mathit{\rho }}^{\mathrm{2}}\right)$ depends on the liquid–solid surface tension (γ_ls) and the density of ice (ρ). The pre-factor A_het is related to the diffusivity of water molecules. The pressure dependence is introduced into this expression by using a formulation for chemical potential difference given by Němec (2013). As with homogeneous freezing, we may approximate that the pre-factor, liquid and ice density, latent heat of fusion, and ice–liquid surface tension remain constant with temperature and pressure on the scales relevant to this work. As long as the heterogeneous compatibility function f_het is not strongly pressure dependent, the approximation of Eq. (1) is valid for heterogeneous, as well as homogeneous, ice nucleation. The heterogeneous compatibility function f_het of the substrate is formulated as a function of the effective contact angle between water and substrate (Lamb and Verlinde, 2011; Zobrist et al., 2007). The roughly linear trend in our data suggests that the contact angle and compatibility function are not strongly dependent on pressure in the negative-pressure regime studied here. This follows from the fact that the derivation of Eq. (1) hinges on approximating many terms as constant along the lines of constant j; if these assumptions are invalid in the pressure regime examined here, we would expect their effects to reveal themselves by producing a highly non-linear trend in pressure–temperature coordinates, which we do not observe.

Table 1Parameters entering Eq. (1) for homogeneous and heterogeneous nucleation for both the mW and ML-mW models. Dashes indicate that a measurement of this value is not known to us.

Sources: ^a Chan et al. (2019). ^b Rosky et al. (2022).

Download Print Version | Download XLSX

The values of T_m, Δν_ls, and l_f used in Eq. (1) to produce the dashed lines in Fig. 2 are listed in Table 1 and correspond to bulk water at 1 atm (no proximity to an interface). As summarized in column four of Table 1, these values predict a slope of $(\mathrm{\Delta }T/\mathrm{\Delta }P{)}_{\mathrm{het}}=-\mathrm{0.069}$ K MPa⁻¹ for the ML-mW model. A line of best fit to the simulated ML-mW heterogeneous-freezing data instead gives a slope of $(\mathrm{\Delta }T/\mathrm{\Delta }P{)}_{\mathrm{het}}=-\mathrm{0.14}$ K MPa⁻¹ (solid red line in Fig.  2a). The predicted slope (dashed red line in Fig. 2a) sits only marginally within the uncertainty bounds of the heterogeneous-nucleation rate data and is 47 % steeper than the best fit of the ML-mW homogeneous-freezing data (see Table 1).

We take some time to consider which factors may contribute to the steeper slope of heterogeneous freezing in the ML-mW model compared to the homogeneous-freezing line. One hypothesis is that Eq. (1) still holds true for heterogeneous ice nucleation, but adjustments need to be made to the values of Δv_ls or l_f used in the equation to achieve quantitative agreement with the heterogeneous-freezing line. While our results support that the trend in $(\mathrm{\Delta }T/\mathrm{\Delta }P{)}_{\mathrm{het}}$ can be approximated as linear in the negative-pressure regime, we are less confident in the quantitative values of Δv_ls and l_f for heterogeneous freezing compared to for homogeneous freezing. In the homogeneous case, the values of T_m, l_f, and Δv_ls are well constrained because they can be calculated from bulk water, which could account for the excellent agreement between Eq. (1) and the homogeneous simulation results of Rosky et al. (2022). For the heterogeneous case, there is more ambiguity around these thermodynamic values because the properties of water near an interface may differ from the bulk properties. We see that, by using bulk thermodynamic values in Eq. (1), we underestimate the slope of $(\mathrm{\Delta }T/\mathrm{\Delta }P{)}_{\mathrm{het}}$ by up to 50 %. Interestingly, this discrepancy between the homogeneous and heterogeneous slopes is seen only in the ML-mW model. Meanwhile the homogeneous- and heterogeneous-freezing lines are nearly parallel for the mW model. This could indicate that the thermodynamic properties of mW water are less influenced near the substrate compared to those of the ML-mW model. This interpretation is supported by evidence from Qiu et al. (2018), who identified that mW water at the mW–carbon interface has thermodynamics similar to those of the bulk liquid.

Another possibility is that the contact angle and the surface free energy (surface tension) between water and substrate exert an influence on the slope of $(\mathrm{\Delta }T/\mathrm{\Delta }P{)}_{\mathrm{het}}$. Evans (1967) showed that a heterogeneous-freezing line for aqueous suspension of ice-nucleating material at positive pressures is actually less steep than the melting-point line and is thus also less steep than the homogeneous-freezing line. Although the mW–substrate and the ML-mW–substrate were designed to have the same freezing-temperature enhancement over the homogeneous-freezing temperatures, they do exhibit different contact angles between the water and substrate. Thus, this could potentially be a factor involved in the discrepancy between ML-mW homogeneous and heterogeneous slopes as compared to the mW case, where the two lines are parallel. Note that, in the previously mentioned study by Qiu et al. (2018), the bulk-like thermodynamics of mW at the substrate interface can be attributed to the nearly 90^∘ contact angle between mW and carbon.

In our simulation results, using bulk water thermodynamic values in Eq. (1) provides a lower bound to the slope of $(\mathrm{\Delta }T/\mathrm{\Delta }P{)}_{\mathrm{het}}$, which can be very useful in estimating the increase in heterogeneous-freezing temperature due to negative pressure in atmospheric and experimental contexts. However, further investigation is needed to identify a robust way to account for the observed steepness of the ML-mW heterogeneous-nucleation slope. Going forward, the linearity and the theoretical prediction of the slope are probably adequate for practical use, given other significant uncertainties in using classical nucleation theory.

3.2 Heterogeneous ice nucleation in water capillary bridges

One way for water to exist stably under negative pressure is through geometric configurations which produce high degrees of negative surface curvature at the air–water interface, such as inside a water capillary bridge. The negative pressure experienced by the water in these cases is a result of Laplace pressure, $\mathrm{\Delta }P={\mathit{\sigma }}_{\mathrm{lv}}\left(\frac{\mathrm{1}}{r_{\mathrm{1}}}+\frac{\mathrm{1}}{r_{\mathrm{2}}}\right)$, where σ_lv is the surface tension between liquid and vapor, and r₁ and r₂ are the radii of curvature of the air–water interface. Negatively valued radii of curvature cause the equilibrium pressure within the water to be smaller than the external environmental pressure, allowing for negative pressure to exist within water that is otherwise in a 1 atm environment. The Laplace pressure associated with different capillary geometries is summarized by Elliott (2021). For the capillary bridge configuration used in this study, shown in Fig. 1c, the expected Laplace pressure within the water is

where θ is the contact angle between water and the substrate. By substituting the above equation into Eq. (1), we obtain an expression to predict the temperature increase for a given nucleation rate j_het as a function of the inverse capillary bridge height:

Given the previous inference that the terms σ_lv and θ do not change significantly with pressure, we expect a linear relationship between freezing temperature and inverse capillary height $\mathrm{1}/h$. Figure 3 shows the freezing temperatures corresponding to an intensive heterogeneous-nucleation rate j_het=10²⁴ m⁻² s⁻¹ inside water capillary bridges with heights h=30, 24, and 18 Å. We find that the data can be described by a linear trend as anticipated. As will be discussed in Sect. 3.4, we have excluded the 18 Å capillary bridge from the current analysis because this scale of confinement of the water between the substrate surfaces causes an increase in ice nucleation rate that cannot be attributed to negative pressure alone.

Figure 3Increase in heterogeneous-freezing temperature ΔT versus $\mathrm{1}/h$ for capillary bridge heights of h=30 (filled green diamond), 24 (filled blue diamond), and 18 Å (open red diamond), as well as for bulk water (filled black circle) in contact with identical ice-nucleating substrates. Excluding the 18 Å capillary (red diamond), which is influenced by confinement effects, the line of constant nucleation rate as a function of inverse capillary height follows a linear trend, as shown by the linear best fit to the data (dashed black line). The 18 Å capillary data are not included in the linear fit to obtain the indicated slope ΔT⋅h.

Download

We can now use the linear slope $-\mathrm{2}{\mathit{\sigma }}_{\mathrm{lv}}\cos \left(\mathit{\theta }\right)\left[\frac{\mathrm{\Delta }{\mathit{\nu }}_{\mathrm{ls}}{T}_{\mathrm{m}}}{l_{\mathrm{f}}}\right]$ from Eq. (4) to analyze the results in Fig. 3. A linear fit to the data, shown by the dashed black line in Fig. 3, produces a slope of $(\mathrm{\Delta }T\cdot h)=\mathrm{9.3}$ K m. We also know from our analysis of $(\mathrm{\Delta }T/\mathrm{\Delta }P{)}_{\mathrm{het}}$ in Sect. 3.1 and Table 1 that the best-fit value of $\left[\frac{\mathrm{\Delta }{\mathit{\nu }}_{\mathrm{ls}}{T}_{\mathrm{m}}}{l_{\mathrm{f}}}\right]$ for the ML-mW model is −0.14 K MPa⁻¹. We substitute these two slopes into Eq. (4) to solve for the value σ_lvcos (θ)=0.042 J m⁻². The surface tension σ_lv of the mW model was reported by Molinero and Moore (2009) to be 0.066 J m⁻², and by following the methods of Li et al. (2009, their Supplement), we found the same value for the ML-mW model. Using this value of σ_lv, we solve for the contact angle θ=50.5^∘. This value of θ is consistent with estimates we obtain by measuring the radius of curvature of the capillary bridge air–water interfaces using a method similar to Giovambattista et al. (2007). With these estimates of σ_lv and θ, we can also use Eq. (3) to calculate the magnitude of Laplace pressure that may be present within the water capillary bridges. The 24 Å capillary bridge has a pressure of −345 atm, and the 30 Å bridge has a pressure of −275 atm.

Our main findings are that negative Laplace pressure created within water capillary bridges increases the temperature of j_het in a manner that is consistent with the linear slope of $(\mathrm{\Delta }T/\mathrm{\Delta }P{)}_{\mathrm{het}}$ combined with the expected Laplace pressure for this capillary geometry. A 24 Å capillary bridge can exhibit a ≈ 3 K increase in heterogeneous-freezing temperature and a ≈ 2 K increase within a 30 Å water capillary bridge.

3.3 Effect of confinement between substrate layers

Capillary theory is expected to remain valid at the nano-scale (Elliott, 2021); however, the ice nucleation rate in water can be affected by confined geometries on these scales (Cao et al., 2019; Roudsari et al., 2022; Hussain and Haji-Akbari, 2021). When confined between two flat surfaces (sometimes referred to as a slit pore), density oscillations in the water induced by the flat interfaces can interfere constructively and influence the ice-forming probability (Cox et al., 2015; Bi et al., 2016; Lupi et al., 2014). When analyzing the increase in freezing temperature within the water capillary bridges, we need to disentangle the effects of confinement from the effects of Laplace pressure. To do so, we ran simulations to observe how confinement alone affects ice nucleation rate on the substrate. We use the configuration shown in Fig. 1b, where the spacing between the substrate surfaces along the z-axis dimension has been reduced to 18, 24, and 30 Å to match the levels of confinement present in our water capillary bridges. In these simulations, water molecules are confined between the substrate surfaces but have no air–water interface that gives rise to Laplace pressure. The pressure in the boxes is set to 1, −500, and −1000 atm for direct comparison against the unconfined configuration across the full range of pressures. Results are plotted in Fig. 4a, showing that 24 and 30 Å separations (blue and green data points) have identical freezing temperatures at all pressures compared to those of the unconfined configuration (black data points). This allows us to conclude that the freezing-temperature enhancement reported previously, e.g., for 24 and 30 Å capillary bridges in Fig. 3, is due to pressure alone. The 18 Å setup (red data points in Fig. 4a) exhibits a significant increase in freezing temperature as a result of the confined geometry. This behavior is indeed evident in the departure of the 18 Å capillary bridge freezing data point from the linear trend in Fig. 3. This is why we have excluded the 18 Å capillary bridge data from the slope analysis in the previous section.

Figure 4(a) Heterogeneous-freezing temperature versus pressure with varying separation between substrates, illustrating the effect of confinement along the z axis. We see that confinement influences the ice nucleation rate for only the 18 Å configuration (red squares). The 24 Å (blue squares) and 30 Å (green squares) confinement configurations exhibit the same freezing temperatures at all pressures as the unconfined reference simulations (black circles). (b) Heterogeneous-freezing temperature versus inverse capillary bridge height for varying widths of the capillary bridge, showing the effect of confinement along the x axis between air–water interfaces. Confining the water within a narrow capillary bridge (orange diamonds) suppresses the increase in freezing temperature that is seen in the capillary bridges that are twice as wide (purple diamonds).

Download

3.4 Freezing locations relative to the air–water interface

We now consider whether there is an analogous confinement effect along the x axis between the two air–water interfaces of our capillary bridge simulations. A narrow capillary bridge configuration can be seen in Fig. 1d. These narrow capillary bridges are 30 Å wide, while the capillary bridges used in our previous results (henceforth referred to as wide capillary bridges) are 60 Å wide. The total surface area of substrate–water contact is kept constant between the narrow and wide capillary bridges by doubling the y dimension of the narrow capillary simulation box. Figure 4b shows heterogeneous-freezing rate temperature as a function of $\mathrm{1}/h$ for the wide capillary bridges (purple diamonds) compared with narrow capillary bridges (orange diamonds). We observe that confining the water within a narrow capillary bridge eliminates the temperature enhancement from Laplace pressure that is observed in the wide capillaries. This result may be connected to the apparent suppression of ice nucleation near the air–water interface, as discussed below.

Figure 5Spatial locations of ice nucleation within water capillary bridges. (a) Locations in a 30 Å tall capillary bridge in the x–z plane. (b) Locations of freezing events in the x–y plane for both the 24 Å and 30 Å capillary bridges combined, viewed from above. (c) Probability density of ice nucleation initiating a distance away from the substrate using freezing locations inside the 24 and 30 Å capillary bridges. The dashed line is a gamma distribution fit. (d) Probability density of ice nucleation initiating a distance away from the air–water interface (red shading in panels a and b) using freezing locations inside the 24 and 30 Å capillary bridges.

Download

Figure 5a shows the x–z plane locations of ice nucleation events within the wide 30 Å tall capillary bridge viewed from the side. The air–water interface is shown in red shading. Figure 5b shows the combined locations of ice nucleation events in the wide 24 and 30 Å tall capillary bridges, as viewed from above in the x–y plane. In this case, the red shading indicates the position of the air–water interface at the substrate. We see that nucleation is preferred near the substrate, as expected, but also that it is suppressed near the air–water interfaces. In Fig. 5c, we plot freezing probability density as a function of distance from the nearest substrate surface using all freezing events in the wide 24 and 30 Å tall capillary bridges. Because this distribution is asymmetric (being affected by the substrate on one side) we fit a gamma distribution to the data to find that 99 % of nucleation events occur at distances between 3.2 and 6.9 Å from the substrate, along the z axis. Figure 5d shows the probability density of ice nucleation as a function of distance from the nearest air–water interface using all freezing events in the wide 24 and 30 Å capillary bridges. Analyzing this spatial distribution of ice nucleation events, we observe that ice nucleation never occurs within approximately 10 Å of the air–water interface. This lack of heterogeneous ice nucleation in the immediate vicinity of the air–water interface could be related to premelting at the ice–vapor interface, as described by Qiu and Molinero (2018).

The narrow capillary bridges being 30 Å wide implies that most of the water molecules are positioned within 10 Å of an air–water interface where no ice nucleation events are observed. In all of the narrow capillary bridges (orange diamonds in Fig. 4b) the heterogeneous-freezing temperature is the same as the 1 atm reference temperature (black data point in Fig. 4b). Furthermore, the enhancement in freezing temperature previously seen in the wide 18 Å tall capillary due to confinement between the substrate layers is eliminated inside the narrow 18 Å tall capillary bridge.

Suppression of ice nucleation in the vicinity of flat air–water interfaces has been noted by Haji-Akbari et al. (2014) using the mW model. They explain this effect through the observation that ice embryos near the interface tend to be less spherical compared to in the bulk, thus imposing a larger ice–liquid surface energy term that inhibits nucleation. For the more detailed TIP4P/Ice model, Haji-Akbari and Debenedetti (2017) saw that ice nucleation was not only suppressed within 10 Å of the air–water interface but also showed evidence of higher ice nucleation rates in the sub-surface region of the interface compared to in the bulk. In Fig. 5d, we see that the freezing locations in the wide capillary bridge simulations hint at a slight preference for ice nucleation to occur between ∼ 20 and 25 Å from the air–water interface. However, when we repeat these simulations with an even wider (120 Å) capillary bridge, we find that this propensity for nucleation in the sub-surface region near the air–water interface is not broadly observed for all geometries (see Appendix B). Whether or not the ML-mW exhibits any surface freezing propensity is inconclusive from our results.

In summary, confinement between the ice-nucleating substrates at scales smaller than 20 Å tends to enhance nucleation in the ML-mW model, whereas confinement between air–water interfaces on scales smaller than 30 Å tends to inhibit heterogeneous nucleation.