Two-dimensional bilayer ice in coexistence with three-dimensional ice without confinement

Domain: earth_science
Published: 2024-07-09
Editorial outcome: minor_revision

Abstract
Icing plays an important role in various physical-chemical process. Although the formation of two-dimensional ice requires nanoscale confinement, two-dimensional bilayer ice in coexistence with three-dimensional ice without confinement remains poorly understood. Here, a critical value of a surface energy parameter is identified to characterize the liquid-solid interface interaction, above which two-dimensional and three-dimensional coexisting ice can surprisingly form on the surface. The two-dimensional ice growth mechanisms could be revealed by capturing the growth and merged of the metastable edge structures. The phase diagram about temperature and pressure vs energy parameters is predicted to distinguish liquid water, two-dimensional ice and three-dimensional ice. Furthermore, the de-icing characteristics of coexisting ice demonstrate that the ice adhesion strength is linearly related to the ratio of ice-surface interaction energy to ice temperature. In addition, for gas-solid phase transition, the phase diagram about temperature and energy parameters is predicted to distinguish gas, liquid water, two-dimensional ice and three-dimensional ice. This work gives a perspective for studying the singular structure and dynamics of ice in nanoscale and provides a guide for future experimental realization of the coexisting ice.

Introduction
Ice is a common state of water, which is formed by the regular arrangement of water molecules. Ice structure and nucleation play an important role in many fields such as materials science, biology, atmospheric science, cryogenic refrigeration and food engineering. Excessive ice accumulation poses great challenges to infrastructure reliability, traffic and human safety. As early as in the 1920s, Bragg, a famous British physicist and X-ray discoverer, together with other scientists, used X-ray to characterize the ice crystal structure. After one hundred years of research and exploration, 18 kinds of 3-dimensional (3D) ice crystal phases have been discovered so far. Since 1980s, predictions of various ice structures have attracted widespread attention. Over the past two decades, computer simulations have confirmed the existence of a large number of low-dimensional ice phases. The nanoscale confinement is considered as a tunable variable that can affect the overall hydrogen bond network formation of water molecules during the freezing transition. Surface ice nucleation and low-dimensional ice under extreme conditions often appear 1–4, and a large number of reported 2-dimensional (2D) ice exists on metal surfaces 5–11, insulating surfaces 12–16, graphite and graphene surfaces 14,17–22, under strong confinement conditions. Molecular dynamics (MD) simulations and experimental studies have shown that 2D crystalline and amorphous ice 20–35 appear in the supercooled state when water molecules are confined between planar hydrophobic surfaces. Water molecules can form different hydrogen bond networks under nanoscale confinement. Water molecules are confined between two parallel surfaces to form single-layer and double-layer square ice 20,21, twisted double-layer hexagonal ice 22,23, double-layer hexagonal and diamond ice 24,27, double-layer interlocked pentagonal ice 28, single-layer diamond ice 29,30, three-layer ice 31, and three-layer heterogeneous fluid 32.

However, fewer studies have so far shown the possibility of 2D ice formation without nanoscale confinement. The existence of ice-like water was observed at the water-mica interface 15,36,37. The formation mechanism of double-layer hexagonal ice (BHI) on graphene/platinum (111) substrate was revealed by the water vapor deposition experiments at 100–130 K 18. In addition, BHI formation can be observed on Au(111), Pt(111) and Ru(0001) substrates at very low temperatures (<140 K) 9,38. The existence of 2D ice was first confirmed by experiments combined with molecular simulation 39. The liquid-to-BHI transition on different smooth surfaces was revealed based on molecular simulation 40,41. In addition, the ice nucleation characteristics and de-icing mechanism of supercooled water droplets were studied by molecular simulation 42–58.

The knowledge on 2D–3D coexisting ice formation and de-icing characteristics at molecular level remains blank. The influences of the surface wettability, temperature, pressure and phase transition mode on the formation of 2D–3D coexisting ice without confinement has not been studied. Our findings fill the gap about the coexistence of 2D bilayer ice and 3D ice without nanoscale confinement. We explore the possibility of controlling the formation of low- and high-dimensional ices by tuning the interactions of a water droplet on a surface. Here, we show a groundbreaking work about the formation and growth of 2D–3D coexisting ice and de-icing characteristics in which the phase diagrams about temperature and wettability are predicted for liquid-solid and gas-solid phase transitions. Changing wettability leads to different ice growth mechanisms, which helps us to better understand water-ice behavior. This study deepens our understanding of 2D–3D coexisting ice without confinement, and is a step forward towards a predictive model for water-ice phase developments which can be used to guide and interpret future experiments. Indeed, the results presented in our work open the debate about whether it is possible to produce and control the growth of ice phases in experiments and reality.

Results
Correlation between the surface wettability and the contact angle
The simulation system was constructed based on LAMMPS simulation package (Supplementary Fig. 1). The smooth surface wettability is characterized by the contact angle. The correlation between the contact angle θ and the water-surface energy parameter ε is clarified (Supplementary Fig. 2). As the energy parameter ε increases, the contact angle decreases. For θ < 90°, the solid surface is hydrophilic. For θ > 90°, the solid surface is hydrophobic. For θ > 150°, the solid surface is superhydrophobic. For θ < 15°, the solid surface is superhydrophilic. When ε exceeds the critical value, the droplet completely wets the surface with θ = 0°.

Liquid-solid phase transition: the effect of temperature on 2D ice in coexistence with 3D ice without confinement
At a constant temperature, the contact angle of nanodroplets decreases with the enhancement of nanodroplet and surface interaction, which is attributed to the larger surface energy and the stronger wettability of the solid surface (Supplementary Fig. 2). As the temperature increases, the potential energy per water molecule basically increases. The potential energy per water molecule at 250 K is significantly higher than that of per water molecule at other temperatures (Supplementary Fig. 3). In the NVT ensemble, the volume of the simulation system box is constant. As the temperature increases, the pressure of the system increases, so does the potential energy of water molecules. Figure 1a shows the change of the potential energy per water molecule with the surface energy parameters at 205 K. The potential energy per water molecule for the energy parameter ε ≤ 0.43 kcal·mol−1 is significantly higher than that of the energy parameter ε > 0.43 kcal·mol−1. As energy parameter ε increases, the potential energy suddenly drops by ~0.43 kcal·mol−1. The potential energy declines sharply between ε = 0.40 kcal·mol−1 and ε = 0.50 kcal·mol−1. Figure 1b displays the mean-square displacement of water molecules on the solid surface at 205 K. The mean-square displacement of ε = 0.50 kcal·mol−1 corresponds to the vertical axis on the right-hand side, which is one order of magnitude higher than that of ε = 0.43 kcal·mol−1 and 0.30 kcal·mol−1 corresponds to the vertical axis on the left-hand side. The mean square displacement of ε = 0.50 kcal·mol−1 increases sharply at first, showing some water molecules of the droplet diffuse rapidly along the wall to form 2D ice. In this process, the liquid water exists at the same time as the 2D ice solid. Then the inflection point appears and tends to increase slowly, indicating that the liquid water molecules nucleate and grow into 3D ice. At this time, both 2D ice solid phase and 3D ice solid phase are present. The whole freezing process has two phase transitions. The first phase transition produces 2D ice, and the second phase transition produces 3D ice. The appearance of the inflection point means that the liquid water molecules nucleate and crystallize in an instant to form 3D ice, and only the solid state exists. For the mean square displacement of ε = 0.43 kcal·mol−1, there is only liquid phase after equilibrium. The critical value of energy parameter ε for phase transition is estimated to be 0.43 kcal·mol−1.

The dynamic processes of nanodroplets on surfaces with different energy parameters show characteristics of the nanodroplets supercooled to nucleation and crystallization (Supplementary Fig. 4 and Fig. 5). The nucleation process is a phase-transition activation process. Overcoming the free energy barrier to form critical nucleus, the nucleus grows spontaneously (Supplementary Fig. 6). The phase transition and crystallization process of nanodroplets on surfaces with different energy parameters can be summarized as: (1) Droplets supercooled; (2) Critical nucleus formation; (3) Rapid phase transition; (4) Ice crystals growth. Near the gas-liquid, liquid-solid and gas-liquid-solid regions, the appearance of critical nucleus is spatially selective. The simulation results show that the surface nucleation process with different wettabilities after quenching displays different characteristics: For ε = 0.10 kcal·mol−1, the formation of critical nucleus preferentially occurs inside the droplet rather than on the solid surface, which is attributed to the superhydrophobic state of the surface (θ > 150°). The extremely low free energy of the solid surface causes the migration of nanodroplet. For ε = 0.20–0.43 kcal·mol−1, nanodroplet nucleation occurs preferentially on solid surface due to the surface being between hydrophobic and hydrophilic state (70° < θ < 150°). With the increase of wettability, the interaction between nanodroplet and the surface enhances, so the contact angle of nanodroplet decreases and the contact radius increases. Nanodroplet on surface with higher wettability is more likely to nucleate and crystallize due to the very high adsorption energy of the surface. Nanodroplet crystallization is not a single cubic or hexagonal structure, but a mixture of cubic and hexagonal structure stacking on each other. The stacking direction is controlled by the growth direction of ice nucleus. In the process of water freezing, a relatively long-lived hydrogen-bonded network is generated between water molecules. For the liquid-solid phase transition, from the perspective of dynamics, the decrease in temperature leads to slower molecular motion and lower kinetic energy. When the kinetic energy of the molecule is small, the potential energy between the molecules is sufficient to confine the molecules, so that the molecules are usually arranged in a regular crystal structure. The nucleation and crystallization of nanodroplets mainly occur on the solid surface (ε = 0.20–1.0 kcal·mol−1), and the latent heat released by local crystallization can provide enough energy for growth and crystallization 52. The variation of the potential energy with time displays the nucleation and crystallization process of nanodroplets after quenching and cooling (Supplementary Fig. 7). The nucleation event is detected by recording the sudden decrease of the potential energy Epot of water molecules during the formation of critical ice nucleus. For ε = 0.10–0.43 kcal·mol−1, the potential energy change of the nucleation and crystallization process of nanodroplets on the solid surface is divided into three continuous characteristic stages: In the first stage, the potential energy of water molecules changes slightly. In the second stage, the potential energy Epot of water molecules decreases sharply, and the critical nucleus is formed. In the third stage, the potential energy of water molecules is in an equilibrium, and the remaining water molecules slowly crystallize.

For 0.43 kcal·mol−1 < ε ≤ 1.0 kcal·mol−1, the nanodroplets first form 2D ice and then form 3D ice during the supercooled process. Taking ε = 0.50 kcal·mol−1 as an example, the dynamic process of the formation of 2D-3D coexisting ice is shown in Fig. 2a. After quenching and cooling, some water molecules on the surface droplets diffuse rapidly to form 2D ice. The 2D ice is formed by two layers of interlocking 5-, 6-, 7-membered rings water molecules. The double-layers are connected by hydrogen bonds. Each water molecule forms three hydrogen bonds with water molecules in the same layer, and forms a hydrogen bond with water molecules in the upper and lower layers. Therefore, all hydrogen bonds are saturated and their structures are very stable 39. 2D ice formation mechanism without nanoscale confinement is that appropriate water-surface interactions can compensate for the entropy loss in the freezing transition process 41. The supercooled unfrozen droplets are promoted by 2D ice. After overcoming the free energy barrier and forming a critical nucleus, the crystal nucleus grows rapidly and crystallizes to form 2D and 3D coexisting ice. The 3D ice is formed by the mixed accumulation of hexagonal ice and cubic ice. The 2D ice underlying growth mechanisms could be revealed by capturing the growth and merged of the metastable edge structures (Supplementary Fig. 8). Firstly, the edges of the 2D ice on both sides are formed by the non-rotating stacking of double-layer 4-, 5-, 6-, 7-membered ring ice. Then, with the diffusion of water molecules, the 2D ice growth edge is composed of double-layer 5-, 6-, 7-membered ring water molecules. The 2D ice boundary has metastable characteristics. At the moment of merging at the edge of the 2D ice double layer, the 6-membered ring in the white circle in the picture contributes a single water molecule, which is connected with the left 5-membered ring, and then forms two metastable 4-membered rings to connect the water molecules on both sides. After that, the two 4-membered ring water molecules become relatively stable 5-membered ring water molecules, and the water molecules on both sides merge and grow. During the growth process, the zigzag growth mode and the armchair growth mode of the 2D ice appears. When the calculation time reaches 200 ns, the 2D-3D coexisting ice structures are relatively stable. The 2D ice is still composed of double-layer 5-, 6-, 7-membered ring water molecules, and the proportion of 5-, 7-membered ring water molecules is basically unchanged, indicating that the 2D-3D coexisting ice is not sensitive to the calculation time when it is sufficiently long (Supplementary Fig. 9). The intensity of the first peak is slightly larger than that of the second peak as shown in Fig. 2b, revealing that the growth of 2D ice starts from the bottom layer (near the solid wall), which is consistent with the results in ref. 39. Figure 2c shows the growth of the number of ice molecules. After quenching and cooling, 2D ice grows rapidly and its number reaches a stable value. After that, the unfrozen droplets overcome the free energy barrier to form 3D ice quickly, and the number of 2D ice and 3D ice molecules reaches the maximum.

Based on the molecular simulation calculation, the critical value of the surface energy parameter is estimated to be 0.43 kcal·mol−1. The evaluation system of liquid water, 2D ice and 3D ice on the temperature and energy parameter ε on the Pt surface is predicted. As shown in Fig. 2d, the influence of multiple temperature conditions and multiple energy parameters on the liquid-solid phase transition of water molecules is clarified. For 210 K < T ≤ 250 K and ε ≤ 0.43 kcal·mol−1, the supercooled nanodroplet is in liquid phase without phase transition. For 210 K < T ≤ 250 K and ε ≥ 0.43 kcal·mol−1, the supercooled nanodroplet is in the coexistence of 2D ice and liquid phase. For 130 K ≤ T ≤ 210 K and ε ≤ 0.43 kcal·mol−1, the supercooled nanodroplet is in solid phase with phase transition, and the 3D ice is disorderly accumulated by hexagonal ice and cubic ice. For 130 K ≤ T ≤ 210 K and ε ≥ 0.43 kcal·mol−1, the supercooled nanodroplets undergoes phase transition, and the ice is in the coexistence of 2D ice and 3D hexagonal ice and cubic ice. The mechanism for the 2D-3D coexisting ice without nanoscale confinement is attributed to overcoming the free energy barrier and the satiation of the Bernal-Fowler ice rules and the appropriate water-surface interaction, which can compensate for the entropy loss caused by the freezing transition process.

Liquid-solid phase transition: the effect of pressure on 2D ice in coexistence with 3D ice without confinement
In order to investigate the influence of pressure on coexisting ice, five pressure conditions were chosen for molecular dynamic simulation at 205 K. In the NVT ensemble, pressure regulation is achieved by adding different numbers of nitrogen molecules (Supplementary Fig. 10). Due to the limited size of the system box, 20, 50, 100, 200 and 300 nitrogen molecules were added to obtain approximate pressures of 1.0, 2.5, 5.5, 11.0 and 17.5 atmospheres, respectively. Figure 3b shows that the contact angle of nanodroplets decreases with increasing gas pressure, and the droplet curvature decreases first and then increases with increasing gas pressure for energy parameter ε = 0.4 kcal·mol−1. Taking pressure of 11.0 atm as an example, for 0.1 kcal·mol−1 ≤ ε < 0.43 kcal·mol−1, the liquid-solid phase transition only generates 3D ice with disordered stacking of hexagonal ice and cubic ice (Supplementary Fig. 11). For 0.43 kcal·mol−1 ≤ ε ≤ 1.0 kcal·mol−1, the liquid-solid phase transition produces coexisting ice, which is composed of 2D ice and 3D ice stacked disorderly with hexagonal ice and cubic ice (Supplementary Fig. 12). With increasing wetting characteristics, the contact angle of droplets decreases, and droplets gradually spread on the wall. After quenching, 2D ice begins to form instantaneously, and 3D ice is formed after overcoming the nucleation energy barrier. The variation of potential energy with time also reflects the occurrence of nucleation and icing process (Supplementary Fig. 13). When the surface energy parameters is 0.5 kcal·mol−1, 2D ice grows and merges slowly for pressure 11.0 atm (and 17.5 atm) as shown in Fig. 3a, which is attributed to the combined effect of surface wettability and gas pressure. After calculation of 200 ns, the 2D ice is still composed of double-layer 5-, 6-, 7-membered ring water molecules.

The phase diagram about energy parameters and pressure is predicted to distinguish 3D ice and 2D-3D coexisting ice by a large number of molecular dynamics simulations. The critical value of the surface energy parameter is estimated to be 0.43 kcal·mol−1, which is the same as that of the nitrogen-free system, indicating that the formation of 2D-3D coexisting ice on the platinum surface is not very sensitive to the applied gas pressure conditions. The effects of multiple pressure conditions and multiple surface wetting characteristics on the liquid-solid phase transition of water molecules are revealed as shown in Fig. 3c. For 1.0 atm ≤ P ≤ 17.5 atm and ε < 0.43 kcal·mol−1, the supercooled droplet undergoes liquid-solid phase transition to solid phase, and the 3D ice is composed of disordered hexagonal ice and cubic ice. For 1.0 atm ≤ P ≤ 17.5 atm and 0.43 kcal·mol−1 ≤ ε ≤ 1.0 kcal·mol−1, liquid-solid phase transition occurs, and the coexisting ice consists of 2D ice and 3D hexagonal ice and cubic ice. For 0.43 kcal·mol−1 ≤ ε < 0.8 kcal·mol−1, the 2D coexisting ice is composed of double-layer 5-, 6-, 7-membered ring water molecules. For 0.8 kcal·mol−1 ≤ ε ≤ 1.0 kcal·mol−1, the 2D ice in coexisting ice is finally composed of a layer of 5-, 6-, 7-membered ring water molecules and a layer of 4-, 6-membered ring water molecules near solid surface (Supplementary Fig. 14). The appearance of the 4-membered square ring water molecules is attributed to the surface being too hydrophilic and gas pressure.

De-icing mechanism of 2D ice in coexistence with 3D ice
The detachment processes of nanoscale coexisting ice from smooth Pt surface are used to investigate the de-icing characteristics (Supplementary Fig. 15). Ice adhesion strength is the ice adhesion force divided by the ice-solid substrate contact area (σ = F/A). The tension is applied to the ice by applying acceleration to all water molecules in the +z direction (perpendicular to the solid substrate), increasing the vertical acceleration at a constant rate, so that the tension is stably loaded throughout the ice. With the superposition of force fluctuations, the attractive force between the coexisting ice and the solid substrate increases linearly (Supplementary Fig. 16a). When the coexisting ice is completely detached from the solid substrate instantaneously, the attractive force decreases suddenly. The coordinates of the intersection point of fluctuating attractive force fitting curve (the black line) and the expected attractive force (the red line) are the detachment force and detachment time (Supplementary Fig. 16b). When the tension is greater than the maximum adhesion between the ice and the solid substrate, the ice will be completely detached from the solid substrate, and the ice detachment is related to the sudden decrease of the adhesion force. MD simulation is able to capture the ice detachment process at the atomic level 55, and provide nanoscale mechanism of ice detachment dynamics. Figure 4a reveals that the ice adhesion force of 2D ice in coexistence with 3D ice is significantly larger than that of 3D ice. Figure 4b shows that the adhesion strength of 2D-3D coexisting ice is smaller than that of 3D ice, which is attributed to the combined effect of ice adhesion force and ice contact area. Ice adhesion strength depends on the ice adhesion force and the ice contact area. The contact area between 2D-3D coexisting ice and solid substrate is 2.35 times that between 3D ice and solid substrate (A2D-3D / A3D = 2.35). The 2D-3D coexisting ice has a larger contact area with a solid substrate, which could provide a higher adhesion force. However, the adhesion force of 2D-3D coexisting ice is no more than 2 times that of 3D ice (F2D-3D / F3D < 2). Therefore, ice adhesion strength for 2D-3D coexisting ice is lower than that of 3D ice. Figure 4c, d shows the dynamic process of coexisting ice and 3D ice detachment from the solid substrate, respectively. The coexisting ice is first detached from the surface at the left junction of the 2D ice and the 3D ice, then the 2D ice and the 3D ice on the left side are detached from the surface, and finally the 2D ice on the right side is detached from the surface. The 3D ice is first detached from the surface at the right edge, becoming increasingly tilted relative to the substrate and eventually detached from the substrate. This angular ice detachment mechanism is consistent with the experimental observation results 54, which is attributed to the sharp edge or corner of the ice having the maximum tensile stress due to its singularity. Therefore, the thermal fluctuation of ice molecules is most likely to cause ice to break away from the edge corner. The calculated ice adhesion strength is two orders of magnitude larger than the experimental ice adhesion strength (≤1 MPa), which is mainly attributed to the fact that the typical loading rate used in MD simulation is about 6–7 orders of magnitude higher than the experimental study 55, the nano-size of ice and the separation stress of ice increase logarithmically with the loading rate 56. If the same loading rate is used in MD simulation and experiment, the de-icing stress of the two methods is not much different.

The adhesion strength of coexisting ice is significantly lower than that of 3D ice for different temperatures and the same energy parameters (Supplementary Fig. 17). At the same temperature, the adhesion strength of 3D ice and coexisting ice is linearly related to the energy parameter ε between nanodroplets and solid substrate, indicating that the surface wettability energy parameter εw-Pt is directly related to the ice adhesion strength. Based on the concept of depletion layer, the linear correlation between cos θ and ε of water droplets on the solid substrate is derived 57. By assuming the homogeneous solid density ρS and liquid density ρL, and ignoring the electrostatic or interfacial entropy, the work H12 of each surface area required to separate the liquid from the solid is calculated theoretically.

H12 = γSV + γLV − γSL = ∫R Z* dz(z − z*)2 u(z) = −πρLρS ∫R Z* dz(z − z*)2 u(z) (1)

1 + cos θ = (γSV + γLV − γSL) / γLV = H12 / γLV ∼ εw-Pt (2)

σ ∼ 1 + cos θ (3)

1 + cos θ* = φs (1 + cos θ) (4)

where γSV, γLV and γSL are the surface tensions between solid-vapor (SV), liquid-vapor (LV) and solid-liquid (SL), respectively. u(z) represents the solid-liquid interaction energy. The LJ interaction simulation analysis is used. Since the LJ potential is linearly related to ε, H12 is also a linear function of ε. Equation (2) is obtained by combining Young’s equation. The ice adhesion strength is linearly related to the cosine of the receding angle of water droplets on the solid surface 58, that is, σ ∼ 1 + cos θrec. The atomic plane does not exhibit contact angle hysteresis, indicating that the receding angle of the water droplet is the same as the equilibrium contact angle. Therefore, the ice adhesion strength and the linear variation of the contact angle cosine of water droplets on the solid surface is given by Eq. (3). From Eqs. (2) and (3), it is obtained that the ice adhesion strength is a linear function of the water-solid substrate interaction force σ ∼ ε. The slope of linear function of coexisting ice adhesion strength and energy parameter ε on solid surface is smaller than that of linear function of 3D ice adhesion strength and energy parameter ε (Supplementary Fig. 17). As the interaction energy between ice and solid substrate decreases, the ice adhesion strength of solid surface decreases.

The ice temperature and the ice-solid substrate interaction force are important parameters affecting the ice adhesion strength. The increase of temperature enhances the fluidity of water molecules in ice, which reduces the ice adhesion strength. Temperature and ice-solid substrate interaction have opposite effect on ice adhesion strength. For a given ice-solid substrate interaction, higher temperature corresponds to lower ice adhesion strength. For different temperatures, the ice adhesion strength is linearly related to the surface energy parameter. Considering the opposite effects of temperature and water-substrate interaction on ice adhesion strength, Fig. 4e shows the linear relationship between ice adhesion strength and dimensionless number ε / (kBT). The linear correlation slope of ε / (kBT) of coexisting ice on solid surface is smaller than that of 3D ice on solid surface, which is consistent with the results shown in Supplementary Fig. 16. With the decrease of ice-solid substrate interaction, the ice adhesion strength on solid surface decreases. At different temperatures, the ice adhesion strength is linearly related to the ice-solid surface interaction. For the coexisting ice and 3D ice on the solid surface, the MD data of ice adhesion strength are fitted to obtain two straight lines, namely, the ice-solid substrate interaction force and ice temperature ratio function.

In NPxyT ensemble, after 200 ns of calculation, the 2D ice composed of metastable double-layer 5-, 6-, 7-membered ring water molecules is transformed into steady-state 6-membered ring 2D ice. 3D ice structure in coexisting ice is not sensitive to the calculation time (Supplementary Fig. 18). At 30 ns, the 2D ice is composed of double-layer 5-, 6-, 7-membered ring water molecules. At 200 ns, the 2D ice is composed of double-layer 6-membered ring water molecules. With the increase of energy parameter εw-Pt, the ice adhesion force at 200 ns is gradually larger than that at 30 ns, which is attributed to the combined effect of surface wetting characteristics, 2D ice structure and ice freezing time (Supplementary Fig. 19).

Gas-solid phase transition: 2D ice in coexistence with 3D ice without confinement
A gas-liquid-solid three-phase system is established to simulate vapor deposition (Supplementary Fig. 20). The dynamic processes of vapor deposition display the characteristics of nucleation and crystallization (Supplementary Fig. 21 and Fig. 22). For different surface wettabilities, dropwise condensation mode or film condensation mode can appear in supercooled process. For ε = 0.10 − 0.20 kcal·mol−1, water vapor is difficult to deposit on the solid surface because the surface is too hydrophobic. For ε = 0.30 − 0.40 kcal·mol−1, with the progress of vapor deposition, the gas-liquid phase variation occurs, and independent nanodroplets appear on the solid surface. The water vapor condenses in dropwise condensation mode, and then the liquid-solid phase transition occurs. The nucleation of the nanodroplet preferentially occurs on the solid surface. Nanodroplet crystallization is a stacking disordered arrangement of cubic and hexagonal structure. For ε = 0.43 – 0.50 kcal·mol−1, the solid surface is hydrophilic. With the progress of vapor deposition, the gas-solid phase transition occurs, and the gas phase molecules diffuse to form 2D ice. The 2D ice is composed of double-layer 5-, 6-, 7-membered ring ice without rotation and stacking. Then the gas-liquid phase transition occurs, and a large number of gas molecules condense in the dropwise condensation mode. With the progress of vapor deposition, the liquid-solid phase transition is activated. The free energy barrier is overcome to form a critical nucleus. The critical nucleus is close to the 2D ice, and the crystal nucleus grows rapidly and crystallizes to form 3D stacking ice with hexagonal structure and cubic structure. 2D ice promotes the nucleation and growth of 3D ice, thus forming 2D and 3D coexisting ice. The formation mechanism of 2D ice without confinement is attributed to the fact that appropriate water-surface interactions can compensate for the entropy loss 41. For ε = 0.60 kcal·mol−1, the solid surface is hydrophilic. Firstly, the gas-solid phase transition occurs to form 2D ice. Secondly, the gas-liquid phase transition occurs and condenses in dropwise condensation mode. Then the liquid-solid phase transition occurs, the crystal nucleus grows rapidly and crystallizes, and the gas-solid phase transition occurs at the same time. The gas molecule directly becomes a component of the 3D ice, and finally only 3D ice is formed. For ε = 0.70 – 1.0 kcal·mol−1, firstly, the gas-solid phase transition occurs to form 2D ice, and then the gas-liquid phase transition occurs and condenses in dropwise condensation mode. Secondly, 2D ice is covered by a complete liquid film, showing the film condensation mode. Finally, the liquid-solid phase transition occurs, forming 3D hexagonal structure and cubic structure disordered accumulation ice. For ε = 0.10 – 0.20 kcal·mol−1, the potential energy of gas molecules keeps in equilibrium without phase transition. For ε = 0.30 − 1.0 kcal·mol−1, the potential energy of gas molecules decreases twice (Supplementary Fig. 23). The first drop is gas-liquid phase transition or the gas-solid phase transition coexists with gas-liquid phase transition, and the second drop is liquid-solid phase transition, which verifies the different surface condensation crystallization phase transition characteristics mentioned above.

For 0.43 kcal·mol−1 ≤ ε ≤ 1.0 kcal·mol−1, the vapor deposition first forms 2D ice during the supercooled process, and then forms 3D ice. For ε = 0.50 kcal·mol−1, the formation process of 2D and 3D coexisting ice is shown in Fig. 5a. As the vapor deposition progresses, the gas-solid phase transition occurs. The gas molecules first form a single 2D double-layer 5-membered ring ice on the solid surface, then form a dispersed 2-dimensional ice on the solid surface composed of the non-rotating stacking of double-layer 5-, 6-, 7-membered ring ice. The dispersed 2D ice merge and grow. The gas molecules are condensed in dropwise condensation mode. The critical nucleus forms and grows rapidly to form 3D ice. Finally, 2D and 3D coexisting ice is formed. Figure 5b shows the relationship between temperature and potential energy. The first potential energy drop corresponds to the coexistence of 2D ice and nanodroplet, including gas-solid phase transition and gas-liquid phase transition. The second potential energy drop corresponds to the formation of 3D ice critical nucleus, including only liquid-solid phase transition. The second potential energy drop corresponds to the temperature jump, indicating that the heat released by local nucleation crystallization can increase the temperature.

Through a large number of the molecular dynamics simulations, the critical value of surface energy parameter is estimated to be 0.43 kcal·mol−1. The evaluation system of temperature and energy parameter ε of gaseous water, liquid water, 2D ice and 3D ice is predicted as shown in Fig. 5c. For 200 K ≤ T ≤ 250 K and ε ≤ 0.20 kcal·mol−1, there is no phase transition. For 200 K ≤ T ≤ 210 K and 0.20 < ε < 0.43 kcal·mol−1, the gas-liquid-solid phase transition occurs. The system only has dropwise 3D ice. For 200 K ≤ T ≤ 210 K and 0.43 ≤ ε < 0.60 kcal·mol−1, the system has 2D and 3D coexisting ice. For 200 K ≤ T ≤ 210 K and 0.60 ≤ ε ≤ 1.0 kcal·mol−1, the gas-liquid-solid phase transition occurs. The system only has film-like 3D ice. For 210 K < T ≤ 250 K and 0.20 < ε < 0.60 kcal·mol−1, the system exists 2D ice and dropwise droplet. For 210 K < T ≤ 250 K and 0.60 ≤ ε ≤ 1.0 kcal·mol−1, 2D ice and film-like water coexist in the system.

Methods
Liquid-solid phase transition: MD simulation of nanodroplet nucleation
The nanodroplet-solid surface simulation system was constructed with a platinum surface as the substrate (Supplementary Fig. 1). The research system contains 2430 water molecules and 4000 platinum atoms. The initial nanodroplets are located on the upper part of the rough surface of the metal platinum, and the platinum atoms with a lattice constant of 3.92 Å are arranged in a face-centered cubic (FCC). The x and y directions are periodic boundary conditions, the z direction is non-periodic boundary conditions, and the reflection boundary conditions are set at the top. The size of simulated box is 15.69 × 3.92 × 11.77 nm. To improve the calculation efficiency and reduce the calculation time and cost, the coarse-grained water molecule model (mW) 42 is used to describe the interaction force between water molecules in nanodroplets, and the interaction force between solid wall and water molecules in nanodroplets is described by 12-6 LJ interaction potential. The energy parameter ε ranges from 0.1 to 1.0 kcal·mol−1 by adjusting the surface wettability energy parameters between solid surface and water molecules, covering the solid surface from superhydrophobicity to superhydrophilicity. The interaction between platinum atoms is based on the 12-6 LJ interaction potential model 59. The corresponding LJ interaction potential parameters are 16.004 kcal·mol−1 and 2.47 Å. In the simulation, each wall atom is connected to its original position through a harmonic spring, that is, U = K[r(t) − r(0)]2, where r(0) is the initial position of the platinum atom, r(t) is the position of the platinum atom at t, and K = 13 eV Å−3 is the spring coefficient. The NVT ensemble is used to control the temperature of the simulation system. The relaxation temperature of the nanodroplet is 300 K. The Verlet velocity method is used to integrate the motion equation. The calculation timestep is 10 fs, and the system equilibrium calculation is 1 ns. The quenching method is used to cool down to 130–250 K, while monitoring the changes in the potential energy and structure of the nanodroplets, and calculating until the nanodroplets nucleate and crystallize to achieve the liquid-solid phase transition of the nanodroplets. The mean square displacement (MSD) at time t is calculated by Eq. (5):

MSD(t) = 1/n ∑ni=0 (ri(t) − ri(0))2 (5)

where ri and r0 are the positions of atom i at time t and 0, respectively.

In the NVT ensemble, pressure regulation is achieved by randomly adding different amounts of nitrogen molecules. A coarse-grained model is used for the nitrogen molecules in the system. Due to the limited size of the system, the number of nitrogen molecules is set to 20, 50, 100, 200 and 300 representing approximately 1.0, 2.5, 5.5, 11.0 and 17.5 atmospheres, respectively. The 12-6 LJ interaction potential model is used for the interaction between nitrogen molecules. The relevant interaction potential parameters are 0.0741 kcal·mol−1 and 3.31 Å. The relaxation temperature of the system is 300 K and the quenching temperature is 205 K. The rest of the setup conditions are the same as above without the addition of nitrogen molecules.

Identification of ice
To identify whether there is ice phase in supercooled water, the bond order parameters are used for quantitative discrimination based on the identification method proposed by Steinhardt et al. 60 and developed by Wolde et al. 61.

q lm(i) = 1/Nb(i) ∑Nb(i) j=1 Y lm(rij) (6)

c(i,j) = ∑l=−l l q lm(i)q* lm(j) / [∑l=−l l q lm(i)q* lm(i)]1/2 [∑l=−l l q lm(j)q* lm(j)]1/2 (7)

where Nb(i) is the number of nearest neighbors of particle i, Y lm(rij) is the spherical harmonic function, rij is the unit vector connecting particle i and one of the nearest neighbors of particle j, and * denotes the complex conjugate function. The c(i,j) with a correlation to q3 is used to distinguish crystals from liquid. The number of staggered bonds [c(i,j) ≤ −0.8] and eclipsed bonds [−0.35 ≤ c(i,j) ≤ 0.25] are used to identify cubic ice, hexagonal ice and liquid water 62.

Ice adhesion strength calculation of 2D and 3D coexisting ice
Based on the LAMMPS molecular simulation platform, 3D ice, 2D and 3D coexisting ice obtained from liquid-solid phase transition (temperature influence section) were used as initial configurations to study the molecular mechanism of de-icing. The simulation system adopts NVT ensemble temperature control. The ice relaxation temperatures are 180, 205 and 230 K. The velocity Verlet method is used to integrate the motion equation. The calculation timestep is 2 fs, and the system equilibrium calculation is 0.50 ns. Then, the acceleration field is applied to all ice atoms in the +z direction, and its magnitude increases linearly with time t. The initial value of acceleration is 0 and increases linearly with 5.734 * 10−9 nm·fs−2. A pulling force is applied to the ice, and the base reacts to the ice. The force FA is monitored by collecting the force of the ice on the base.

Gas-solid phase transition: molecular dynamics simulation of ice nucleation
Supplementary Fig. 19 shows a gas-liquid-solid three-phase system established by the LAMMPS molecular simulation platform. The system contains 5197 mW water molecules and the solid wall contains 8000 Pt atoms. The simulation system mainly includes high temperature wall (heat source), high temperature liquid layer, and low temperature wall (cold source). The sizes of nanochannels in three directions of xyz are 15.69 × 3.92 × 313.76 nm. The Pt-FCC structure is used on the high and low temperature solid walls, and the x, y and z directions are periodic boundary conditions. The interaction force between water molecules in the liquid layer and the wall is described by the 12-6 LJ interaction potential. Changing the interaction strength between the nanodroplet and the solid wall, εw-Pt can realize the regulation of the wall from superhydrophobic to superhydrophilic in the range of 0.1–1.0 kcal·mol−1. Gas-solid phase transition surface properties are the same as liquid-solid phase transition surface properties. The simulation system uses the NVT ensemble to control the temperature, and its relaxation temperature is 500 K. The velocity Verlet method is used to integrate the motion equation. The calculation timestep is 10 fs, and the system equilibrium calculation is 1 ns. Then the NVE ensemble is used for the water molecules. Finally, the solid surface temperature is reduced to the target temperature by quenching. At the same time, the change of potential energy and liquid structure is monitored, and the calculation is carried out until crystallization, so as to obtain the gas-solid phase transition.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors conducted a molecular dynamic simulation study that demonstrated the coexistence of\ntwo-dimensional bilayer ice and three-dimensional ice without confinement. They found that above\na critical value of the surface energy parameter, this coexistence occurs. The study explored the\ncoexistence of ice during liquid-solid and gas-solid phase transitions. The authors discussed the\ndeicing mechanism for 2D-3D coexisting ice, which has lower ice adhesion strength than 3D ice.\nHowever, the mechanism that leads to the coexistence of 2D-3D with the increase in surface\nenergy parameters showing low ice adhesion strength needs further clarification.\nComments and questions:\n1. What is the novelty and motivation of the work? Previous studies have reported the formation of\n2D ice without nanoscale confinement including experimental work with molecular simulations\nshowing practical scenarios (Ref. 36, 42-44). Figure 1 depicts the potential energy variation with\nan increase in energy parameters, which exhibits a trend similar to that reported previously for 2D\nice formation (Ref. 44). What is the research gap? How would the present work provide new\nknowledge to the field?\n2. What are the practical scenarios or conditions in which 2D-3D coexisting ice can form without\nconfinement? How can the authors justify the validity of their study, given that it is purely a\nmolecular simulation study?\n3. What mechanism leads to the coexistence of 2D-3D ice with an increase in the surface energy\nparameter (Fig. 2)? The change in surface wettability from superhydrophobic to superhydrophilic\nhas been well studied with variations in the surface energy parameter (Ref. 44). Why does the\npotential energy increase with temperature (Fig. S3)? The authors should clarify the reasons for\nthe results mentioned in the present work.\n4. Figure 3 shows that 3D ice has higher ice adhesion strength than 2D-3D coexisting ice. What is\nthe reason for the lower ice adhesion strength for 2D-3D coexisting ice? The 2D-3D coexisting ice\nhas a large contact area with a solid substrate which could provide a higher adhesion force.\nMoreover, the authors mention, “Fig. 3(a) and (b) reveal that the ice adhesion force of 2-\ndimensional ice in coexistence with 3-dimensional ice is significantly larger than that of 3-\ndimensional ice”. The explanation contradicts the results shown. The authors should clarify this\ndiscrepancy.\n5. The phase diagram for the gas-solid phase transition (Figure 4c) shows that 3D ice exists for\nenergy parameters from 0.6-1 kcal/mol at low temperatures (< 210 K). However, the phase\ndiagram for the liquid-solid phase transition (Figure 2d) shows the coexistence of 2D-3D ice for\nhigher surface energy parameters and similar temperature conditions. Given that the gas-solid\nphase transition involves gas-liquid and liquid-solid transitions, why does 3D ice occur during the\ngas-solid phase transition for the higher value of the surface energy parameter?", "reviewer_2": "In the manuscript \"Two-dimensional bilayer ice in coexistence with three-dimensional ice without\nconfinement,\" the authors explore the possibility of controlling the formation of low- and high-\ndimensional ices by tuning the interactions of a water droplet on a surface. Understanding the\nconditions under which 2D ices can be produced is an active research topic that still needs to be\nfurther explored. In this regard, the authors offer an alternative perspective that differs from the\nusual ones based on molecular confinement. The fact that changing wettability results in a\ndifferent ice growth mechanism may have important consequences for our understanding of\nwater/ice behavior, as well as for many applications. Indeed, the results presented in the work\nopen a debate about whether it is possible to produce and control the growth of other ice phases\nthrough this method.\nThe manuscript is overall quite convincing. The computational level is well tested and accurate for\nthe results presented. Indeed, the authors have already published an article where a similar\nmethodology was used with reliable results. The way different parameters affect their simulations\nis also well supported in the ESI file. Presenting wettability-temperature phase diagrams offers a\nway to probe their results, which may stimulate the experimental community. I think there is\nexcellent ground for a Nat. Comm. article, but also that some complements and clarifications are\nneeded.\n1. The importance of the work should be addressed in more detail in the introduction. A nice\nbackground on confinement studies on 2D ices is presented as well as some results without\nconfinement. However, the topic of 2D/3D coexistence and the novelty of their work is barely\ndiscussed. An explanation of the potential implications might increase the visibility of their work.\n2. During the result discussion, citations should be added to reinforce some of the conclusions the\nauthors present. For instance, it is said in both the liquid-solid and gas-solid analysis that\n\"thermodynamic analyses show that the tendency to form 2-dimensional ice without nanoscale\nconfinement is that appropriate water-surface interactions can compensate for the entropy loss in\nthe freezing transition process.\" Which thermodynamic analysis? The authors do not present these\nresults in the manuscript.\n3. Regarding the simulations, the systems consist of a drop on a surface. Being a curved surface,\nthere is a pressure difference inside the drop and outside, according to the Young-Laplace\nequation, which affects thermodynamic conditions such as vapor pressure. How is this\nphenomenon included in the model, and how does it affect their stability? Likewise, for a curved\nsurface the Gibbs-Thomson equation predicts a melting point depression. Please provide some\nestimation of the changes in the vapor pressure and melting temperature.\n4. The authors use an NVT ensemble in their simulations, however, nucleation of other competing\nice phases is produced when considering different pressure conditions. Can their studies be\nextended to other pressure conditions?\n5. In their previous work (ref. 44 of the manuscript. Proc. Natl. Acad. Sci. U.S.A. 116, 16723-\n16728 (2019)), the authors comment that using NpxyT simulations, the number of pentagonal,\nheptagonal, and octagonal rings in the interior part of the bilayer ice decreases and even goes to\nzero at longer simulation times. How does this effect affect the coexistence of 2D/3D ice? Since\nthe structure of the 2D ice changes, could the adhesion energy of the 2D ice (de-icing mechanism)\nbe modified? This point should be addressed in the manuscript as it might have implications in the\nconclusions of the work.\n6. Unequivocally distinguishing ice types and water local structure is crucial in describing\nnucleation phenomena. The performance of a given order parameter has been widely discussed in\nthe literature, since it may yield different results. (J. Chem. Phys. 144, 034501 (2016)) How do\ndifferent order parameters affect the formation of 2D and 3D ice in their simulations? Have the\nauthors used a minimum mislabeling procedure to determine the values used for the q6 and q3\nparameters? The influence of these parameters on the results presented in the manuscript should\nbe stated, if any.\n7. I believe that the manuscript would benefit from a discussion on how the simulations can be\ntranslated into experiments. For instance, mentioning that changes in the ε parameter could\npotentially be achieved through modifications of lattice planes or surface characteristics would\nmake the paper more accessible to experimentalists and broaden its audience.\nMinor comments:\n1. There are some sentences that are quite similar in the solid-gas and solid-liquid discussion. It\nwould be beneficial for the authors to unify the text.\n2. In Figure 1 (left), the red axis is not explained in either the text or the caption. Please include\nan explanation.\n3. The citation format needs to be revised and unified.", "reviewer_3": "The authors are using computer simulations to observe the solidification process of water on solid\nsurfaces.\nWhile I think their point of view is interesting, they do not seem to have investigated it in sufficient\ndepth. They give various descriptions of the crystallization process, but they are not accompanied\nby data to corroborate them, so we can only guess them from the context. By the time I got to\npage 5, I was confronted with various points of uncertainty. I think this study (or this paper) is\nincomplete.\nLine 94: “The potential energy per water molecule at 250 K is significantly higher than that of per\nwater molecule at other temperatures. “: This sentence is unclear; it mentions 250 K, but without\nany data, I have no idea what they are referring to.\nLine 109 Figure 1(b) The scale on the right has no explanation and makes no sense.\nRegarding the plot of 0.5 kcal/mol in red dots, the authors say it shows that liquid and solid\ncoexist, but if they coexist, wouldn't the graph reflect only the liquid character with a large\ndiffusion coefficient?\nAround line 124: various explanations are given, but no evidence seems to be presented. For\nexample,\nLine 122: They state that nucleation occurs in the droplet and not on the solid surface because the\nsurface is superhydrophobic, but they do not show a causal relationship between the two, which\nreads as just an unfounded prediction.\nLine 128: “With increasing wettability, the contact angle of nanodroplet decreases and the contact\nradius increases, which enhances the interaction between droplet and the surface.” I think the\ncausal relationship is the opposite.\nLine 131: “Nanodroplet on surface with higher wettability is more likely to nucleate and crystallize\ndue to the very high adsorption energy of the surface” It is just a guess and no evidence is\nprovided.\nLine 134: I don't know what a “large-scale hydrogen bond” is.\nLine 135: The phenomenon of molecules aligning to form crystals is not limited to water. I don't\nthink water is particularly prone to ordering, but if you claim that, you need evidence.\nLine 136: It is an obvious statement over 10 lines.\nLine 143: This is not true. The most stable phase is the phase with the lowest free energy, not the\nlowest energy.\nLine 145 is correct, but it is an obvious thermodynamic fact and I am not sure what they are trying\nto emphasize.\nLine 147: “The nucleation and crystallization of nanodroplets mainly occur on the solid surface”\nSince no conditions are set and no evidence is given, it is not clear what kind of situation is\nenvisioned in the statement."}, "author_response": "Response Sheet\nTitle: Two-dimensional bilayer ice in coexistence with three-\ndimensional ice without confinement\nPaper number: NCOMMS-23-05786\nWe are glad to see that all reviewers have recommended our manuscript for publication.\nWe would like to thank the reviewers for their constructive suggestions. In the following,\nwe have responded point-by-point to all their comments and questions and indicated\nthe corresponding changes made in the revised manuscript. All the changes are marked\nin blue in the revised manuscript.\nReviewer #1 (Remarks to the Author):\nThe authors conducted a molecular dynamic simulation study that demonstrated the\ncoexistence of two-dimensional bilayer ice and three-dimensional ice without\nconfinement. They found that above a critical value of the surface energy parameter,\nthis coexistence occurs. The study explored the coexistence of ice during liquid-solid\nand gas-solid phase transitions. The authors discussed the deicing mechanism for 2D-\n3D coexisting ice, which has lower ice adhesion strength than 3D ice. However, the\nmechanism that leads to the coexistence of 2D-3D with the increase in surface\nenergy parameters showing low ice adhesion strength needs further clarification.\nReply: We appreciate these positive comments about our work. In the following,\nwe address the reviewer’s questions and indicate revisions made in the revised\nmanuscript and Supplementary Information (SI) accordingly.\nComments and questions:\n1. What is the novelty and motivation of the work? Previous studies have reported the\nformation of 2D ice without nanoscale confinement including experimental work with\n1\nmolecular simulations showing practical scenarios (Ref. 36, 42-44). Figure 1 depicts\nthe potential energy variation with an increase in energy parameters, which exhibits a\ntrend similar to that reported previously for 2D ice formation (Ref. 44). What is the\nresearch gap? How would the present work provide new knowledge to the field?\nReply: Thank you for your constructive question.\na) The knowledge on 2D-3D coexisting ice formation and de-icing characteristics\nat molecular level remains blank. The influences of the surface wettability,\ntemperature, pressure and phase transition mode on the formation of 2D-3D\ncoexisting ice without confinement has not been studied. Our findings fill the\ngap about the coexistence of 2D bilayer ice and 3D ice without nanoscale\nconfinement. We explore the possibility of controlling the formation of low-\nand high-dimensional ices by tuning the interactions of a water droplet on a\nsurface. Here, we show a groundbreaking work about the formation and growth\nof 2D-3D coexisting ice and de-icing characteristics in which the phase\ndiagrams about temperature and wettability are predicted for liquid-solid and\ngas-solid phase transitions. Changing wettability leads to different ice growth\nmechanisms, which helps us to better understand water-ice behavior. This study\ndeepens our understanding of 2D-3D coexisting ice without confinement, and\nis a step forward towards a predictive model for water-ice phase developments\nwhich can be used to guide and interpret future experiments. Indeed, the results\npresented in our work open the debate about whether it is possible to produce\nand control the growth of ice phases in experiments and reality.\nb) We used the same method, i.e. by tuning the interactions of a water droplet on\na surface, to determine a critical value of the surface energy parameter\ncharacterizing the liquid-solid interface interaction. Above this critical value,\n2D and 3D coexisting ice can be formed on the surface. This finding is new and\nsignificant for exploring the possibility of controlling the formation of low- and\nhigh-dimensional ices by tuning the interactions of a water droplet on a surface.\nc) Our study has elucidated the physical mechanism behind the formation of 2D-\n3D coexisting ice, which offers further insight into how the formation of 2D-3D\n2\ncoexisting ice is related to the water-surface interactions. Our findings provide\na molecular-level explanation of 2D-3D coexisting ice on surfaces without\nconfinement, thereby offering a guide for growing 2D-3D coexisting ice on\ndifferent solid surfaces. Our work will hopefully motivate future experimental\nstudies of the formation of low- and high-dimensional ice on a variety of\nrealistic surfaces. 2D ice on the surface can promote or inhibit the formation of\n3D ice, which has potential application value for the design and development of\nanti-icing materials.\nThe updated content above have been incorporated to the revised manuscript on\nlines 68-83, 468-475.\n2. What are the practical scenarios or conditions in which 2D-3D coexisting ice can\nform without confinement? How can the authors justify the validity of their study, given\nthat it is purely a molecular simulation study?\nReply: Thank you for your question.\na) The growth of 2D ice has been widely observed experimentally on graphene [1],\nAu(111) surface [2], Ru(0001) substrate and Pt(111) surface [3] at low\ntemperature without confinement. For graphene surface, thin films of water\nwere deposited on the substrate at temperatures ranging from 20-152K. The 2D\nice was only observed for 100-135 K, while 3D ice was observed at T>140K.\nFor Au(111) surface, the 2D ice was observed for T=120K. For comparison, the\n2D ice can be observed on the Pt(111) surface and Ru(0001) substrate for\nT＜140K. Hence, the experimental results for Pt(111) surface are consistent\nwith our predictions for 130K≤T<140K, that is, there is 2D ice. By adjusting the\ndroplet size and surface wettability, the practical working conditions of 2D-3D\ncoexisting ice can be obtained by experiment. The practical conditions for 2D-\n3D coexisting ice without confinement are estimated as surface energy\nparameter 0.43 kcal/mol<ε ≤1.0 kcal/mol and temperature range\nw-Pt\n130K≤T<140K. The contact angle of the water droplet on the platinum metal\nsurface is θ<75°, which is in the hydrophilic and superhydrophilic state. It is\n3\nimportant to obtain the hydrophilic and superhydrophilic surface by\nconstructing rough structure and adjusting wettability. There are two main ways\nto prepare hydrophilic and superhydrophilic surface. One is light-induced\nmethod, which uses the photocatalytic process of photocatalytic materials such\nas TiO films to produce hydrophilicity and superhydrophilicity. The second is\n2\nto construct microscopic rough structures on the hydrophilic surface to make it\nmore hydrophilic, by increasing the roughness of the substrate surface.\nb) For Pt(111) surface, 2D ice can be observed in the experimental temperature\nrange of T<140K [3]. For 130K≤T<140K, the experimental results are already\nconsistent with our predictions for liquid-solid phase transition by MD\nsimulation. In addition, the intensity of the right low peak is significantly greater\nthan that of the left low peak as shown in Fig. 2b, revealing that the growth of\nthe 2D double ice begins at the bottom, which is consistent with the results in\nReference [40]. The 3D ice is first detached from the surface at the right edge,\nbecoming increasingly tilted relative to the substrate and eventually detached\nfrom the substrate. This angular ice detachment mechanism is consistent with\nthe experimental observation results [55]. Actually, we have validated the\nmolecular simulation study in our previous work [4] (3.1. Validation of the\nsimulation method. To further verify the water model and the simulation setup,\na homogeneous ice nucleation simulation is performed and its nucleation rate is\nestimated as 1.92±0.12×1032 m−3s−1, which is comparable to the order of\nmagnitude of the simulation results in Ref. [26]). In addition, we can find similar\nmethods from previous studies and five cited references [6-10] below are listed\nas examples,\n[1] Kimmel, G.A. et al. No confinement needed: Observation of a metastable hydrophobic\nwetting two-layer ice on graphene. J. Am. Chem. Soc. 131, 12838-12844 (2009).\n[2] Ma, R. et al. Atomic imaging of edge structure and growth of a two-dimensional\nhexagonal ice. Nature 577, 60-63 (2020).\n[3] Maier, S., Lechner, B.A., Somorjai, G.A. & Salmeron, M. Growth and structure of the\nfirst layers of ice on Ru(0001) and Pt(111). J. Am. Chem. Soc. 138, 3145-3151 (2016).\n[4] Jiang J., Li G.X., Sheng Q. & Tang G.H. Microscopic mechanism of ice nucleation:\nThe effects of surface rough structure and wettability. Appl. Surf. Sci. 510, 145520 (2020).\n4\n[5] Matsumoto, M., Saito, S. & Ohmine, I. Molecular dynamics simulation of the ice\nnucleation and growth process leading to water freezing. Nature 416, 409-413 (2002).\n[6] Moore, E.B. & Molinero, V. Structural transformation in supercooled water controls\nthe crystallization rate of ice. Nature 479, 506-509 (2011).\n[7] Lupi, L. et al. Role of stacking disorder in ice nucleation. Nature 551, 218-222 (2017).\n[8] Zhu, C.Q. et al. Direct observation of 2-dimensional ices on different surfaces near room\ntemperatature without confinement. Proc. Natl. Acad. Sci. USA 116, 16723-16728 (2019).\n[9] Bi, Y.F., Cao, B.X. & Li, T.S. Enhanced heterogeneous ice nucleation by special surface\ngeometry. Nat. Commun. 8, 15372 (2017).\nThe updated content have been incorporated to the revised manuscript on lines\n199-202 and 308-310.\n3. What mechanism leads to the coexistence of 2D-3D ice with an increase in the\nsurface energy parameter (Fig. 2)? The change in surface wettability from\nsuperhydrophobic to superhydrophilic has been well studied with variations in the\nsurface energy parameter (Ref. 44). Why does the potential energy increase with\ntemperature (Fig. S3)? The authors should clarify the reasons for the results mentioned\nin the present work.\nReply: Thank you for your constructive question.\na) The coexistence of 2D-3D ice start to appear as the surface energy parameter\nincreases over a critical value, and the specific critical value depends on the\nsurface wettability and temperature. The mechanism for the 2D-3D coexisting\nice without nanoscale confinement is attributed to overcoming the free energy\nbarrier and the satiation of the Bernal-Fowler ice rules and the appropriate\nwater-surface interaction, which can compensate for the entropy loss caused by\nthe freezing transition process. For the liquid-solid phase transition part, the\npurpose of our investigation is to predict the phase diagram of temperature and\nenergy parameters to distinguish liquid water, 2D double-layer ice and 3D ice.\nIn addition, 2D-3D coexisting ice is also obtained. However, Ref. 42 studied 2D\nice formation on various surfaces and the dependence of the 2D crystalline\nstructure on the hydrophobicity and morphology of the underlying surface.\n5\nb) In the NVT ensemble, the volume of the simulation system box is constant. As\nthe temperature increases, the pressure of the system increases, so does the\npotential energy of water molecules.\nThe updated content above have been incorporated to the revised manuscript on\nlines 224-227 and 101-104.\n4. Figure 3 shows that 3D ice has higher ice adhesion strength than 2D-3D coexisting\nice. What is the reason for the lower ice adhesion strength for 2D-3D coexisting ice?\nThe 2D-3D coexisting ice has a large contact area with a solid substrate which could\nprovide a higher adhesion force. Moreover, the authors mention, “Fig. 3(a) and (b)\nreveal that the ice adhesion force of 2-dimensional ice in coexistence with 3-\ndimensional ice is significantly larger than that of 3-dimensional ice”. The explanation\ncontradicts the results shown. The authors should clarify this discrepancy.\nReply: Thank you for your constructive question.\na) Ice adhesion strength is the ice adhesion force divided by the ice-solid substrate\ncontact area (σ=F/A). Ice adhesion strength depends on the ice adhesion force\nand the ice contact area. The contact area of 2D-3D coexisting ice and solid\nsubstrate is 2.35 times that of 3D ice and solid substrate (A /A =2.35). The\n2D-3D 3D\n2D-3D coexisting ice has a large contact area with a solid substrate, which could\nprovide a higher adhesion force. However, the adhesion force of 2D-3D\ncoexisting ice is no more than 2 times that of 3D ice (F /F <2). Therefore,\n2D-3D 3D\nice adhesion strength for 2D-3D coexisting ice is lower than that of 3D ice.\n  F A  F A  2.35F F 1\n3D 2D3D 3D 3D 2D3D 2D3D 3D 2D3D\nb) Apologies for the writing mistake. We have revised the sentence to “Fig. 4a\nreveals that the ice adhesion force of 2D ice in coexistence with 3D ice is\nsignificantly larger than that of 3D ice. Fig. 4b shows the adhesion strength of\n2D-3D coexisting ice is smaller than that of 3D ice, which is attributed to the\ncombined effect of ice adhesion force and ice contact area.”\nlines 292-302.\n5. The phase diagram for the gas-solid phase transition (Figure 4c) shows that 3D ice\nexists for energy parameters from 0.6-1 kcal/mol at low temperatures (< 210 K).\nHowever, the phase diagram for the liquid-solid phase transition (Figure 2d) shows the\ncoexistence of 2D-3D ice for higher surface energy parameters and similar temperature\nconditions. Given that the gas-solid phase transition involves gas-liquid and liquid-solid\ntransitions, why does 3D ice occur during the gas-solid phase transition for the higher\nvalue of the surface energy parameter?\nReply: Thank you for your constructive question.\nFor 200 K≤T≤210 K and 0.6≤ε≤1.0 kcal/mol, gas-solid phase transition shows\nthat 3D ice exists, and the liquid-solid phase transition shows the presence of\n2D-3D coexisting ice. The difference is mainly attributed to the phase change\nform, supercooled degree and droplet shape before freezing. For the gas-solid\nphase transition, the solid surface temperature is reduced from 500 K to the\ntarget temperature by quenching. The gas-solid phase transition involves gas-\nliquid and liquid-solid transitions. The water vapour first condenses into film\nwater and then freezes in the form of 3D ice. For the liquid-solid phase transition,\nthe nanodroplets are quenched from 300K to the target temperature. The water\nmolecules close to the wall in the nanodroplets grow rapidly to form 2D ice, and\nthe remaining water molecules exist in the form of dropwise-water. When the\ncritical nucleus is formed by overcoming the free energy barrier, the 3D ice\nnucleus grows spontaneously and crystallizes. Therefore, 2D ice and 3D ice\ncoexist.\nIn the manuscript \"Two-dimensional bilayer ice in coexistence with three-dimensional\nice without confinement,\" the authors explore the possibility of controlling the\nformation of low- and high-dimensional ices by tuning the interactions of a water\ndroplet on a surface. Understanding the conditions under which 2D ices can be\nproduced is an active research topic that still needs to be further explored. In this regard,\nthe authors offer an alternative perspective that differs from the usual ones based on\nmolecular confinement. The fact that changing wettability results in a different ice\ngrowth mechanism may have important consequences for our understanding of\nwater/ice behavior, as well as for many applications. Indeed, the results presented in\nthe work open a debate about whether it is possible to produce and control the growth\nof other ice phases through this method.\nThe manuscript is overall quite convincing. The computational level is well tested and\naccurate for the results presented. Indeed, the authors have already published an article\nwhere a similar methodology was used with reliable results. The way different\nparameters affect their simulations is also well supported in the ESI file. Presenting\nwettability-temperature phase diagrams offers a way to probe their results, which may\nstimulate the experimental community. I think there is excellent ground for a Nat.\nComm. article, but also that some complements and clarifications are needed.\nReply: We would like to thank the reviewer for their very positive and constructive\ncomments. In the following, we address the reviewer’s questions and indicate\nrevisions made in the revised manuscript and Supplementary Information (SI)\naccordingly.\n1. The importance of the work should be addressed in more detail in the introduction.\nA nice background on confinement studies on 2D ices is presented as well as some\nresults without confinement. However, the topic of 2D/3D coexistence and the novelty\nof their work is barely discussed. An explanation of the potential implications might\nincrease the visibility of their work.\nReply: Thank you for your constructive question.\nice formation without confinement, while the 2D-3D coexisting ice formation\nand de-icing characteristics at molecular level remains blank. The influence of\nthe surface wettability, temperature, pressure and phase transition mode on the\nformation of 2D-3D coexisting ice without confinement has not been\nsystematically studied. Our novel findings fill the gap where 2D bilayer ice\ncoexists with 3D ice without nanoscale confinement. We explore the possibility\nof controlling the formation of low- and high-dimensional ices by tuning the\ninteractions of a water droplet on a surface. Here, we report a groundbreaking\nwork about the formation and growth of 2D-3D coexisting ice and de-icing\ncharacteristics in which the phase diagrams about temperature and wettability\nare predicted for liquid-solid and gas-solid phase transitions. Changing\nwettability leads to different ice growth mechanisms, which helps us to better\nunderstand water/ice behavior. This study deepens our understanding of 2D-\n3D coexisting ice without confinement and is a step forward towards a\npredictive model for the formation and growth of other ice phases. Such a\nmodel can help to guide and interpret future experiments. Indeed, the results\npresented in our work open a debate about whether it is possible to produce and\ncontrol the growth of other ice phases through our method. Physical mechanism\nrelated to the formation of 2D-3D coexisting ice is elucidated, which offers\ninsight into how the 2D-3D coexisting ice formation is related to the water-\nsurface interactions. Our findings provide a molecular-level understanding of\n2D ice formation on surfaces without confinement, thereby offering a guide for\ngrowing 2D ices on different solid surfaces. The formation of 2D-3D coexisting\nice without confinement revealed in this work will hopefully motivate future\nexperimental studies on a variety of realistic surfaces. 2D ice has an important\ninfluence on the growth of 3D ice. If there is a 2D ice, the 3D ice will attach\nand grow on the surface and be very stable. However, if there is no 2D ice, the\nformed 3D ice has a small contact surface with the surface and is easily blown\n9\naway by the wind. Therefore, we can design and develop anti-icing materials\nmore specifically according to the structure of 2D and 3D coexisting ice.\nThe updated content above have been incorporated to the revised manuscript on\nlines 68-83, 468-475.\n2. During the result discussion, citations should be added to reinforce some of the\nconclusions the authors present. For instance, it is said in both the liquid-solid and gas-\nsolid analysis that \"thermodynamic analyses show that the tendency to form 2-\ndimensional ice without nanoscale confinement is that appropriate water-surface\ninteractions can compensate for the entropy loss in the freezing transition process.\"\nWhich thermodynamic analysis? The authors do not present these results in the\nmanuscript.\nReply: Thank you for your kind comments. We have cited the relevant reference [42]\nto reinforce the conclusions.\nThe updated content have been incorporated to the revised manuscript on lines\n177-179 and 401-403.\n3. Regarding the simulations, the systems consist of a drop on a surface. Being a curved\nsurface, there is a pressure difference inside the drop and outside, according to the\nYoung-Laplace equation, which affects thermodynamic conditions such as vapor\npressure. How is this phenomenon included in the model, and how does it affect their\nstability? Likewise, for a curved surface the Gibbs-Thomson equation predicts a\nmelting point depression. Please provide some estimation of the changes in the vapor\npressure and melting temperature.\nReply: Thank you for your kind comments.\na) Pressure difference inside the drop and outside is indeed an existing\nphenomenon especially for a nanodroplet based on the Young-Laplace equation.\nFirstly, the vapor pressure mentioned above is almost independent of the droplet\nradius based on previous research result [1]. If the effect of droplet size is\nconsidered, it should have a higher vapor pressure than the flat liquid surface.\nsystem setup in MD simulation. Considering the mean free path of vapor\nmolecules, the rarefaction effect will appear in such a simulation system. It\nshould be noted that for a droplet system in MD simulation the rarefaction effect\nis almost unavoidable due to the limited size of the simulation system. Therefore,\ncombining the Young-Laplace equation and rarefaction effect, the vapour-\npressure independence of droplet size may appear. Of course, current research\nhas not yet clearly verified this speculation, which should be an interesting topic\nfor future research. Here, we just provide a reasonable speculation based on\nprevious research results.\nb) The melting point depression for a particle with a certain radius of curvature\ncould be described using the Gibbs-Thomson equation. However, the original\nGibbs-Thomson equation may not be suitable for this simulation size, according\nto a previous study [2]. At least one correction for surface tension related to the\ndroplet size is required to accurately describe the melting point depression.\nFurthermore, in our simulation, the decrease in the melting point does not affect\nthe ice progress. For example, we set a lower icing temperature in which the ice\nstructure does not show any differences from the results at the current\ntemperature. However, as mentioned by the reviewer, the Gibbs-Thomson\nequation is quite important for the prediction of melting point. But, this is not\nthe focus of this study in which a lower temperature ensuring the occurrence of\nthe freezing process is sufficient.\n[1] Langroudi S M M, Ghassemi M, Shahabi A, et al. A molecular dynamics study of\neffective parameters on nano-droplet surface tension. Journal of Molecular Liquids, 2011,\n161(2): 85-90.).\n[2] Wu N, Lu X, An R, et al. Thermodynamic analysis and modification of Gibbs–Thomson\nequation for melting point depression of metal nanoparticles. Chinese Journal of Chemical\nEngineering, 2021, 31: 198-205.\n4. The authors use an NVT ensemble in their simulations, however, nucleation of other\ncompeting ice phases is produced when considering different pressure conditions. Can\ntheir studies be extended to other pressure conditions?\nmechanism of pressure on coexisting ice in second section, as follows,\nLiquid-solid phase transition: the effect of pressure on 2D ice in\ncoexistence with 3D ice without confinement\nIn order to investigate the influence of pressure on coexisting ice, five pressure\nconditions were chosen for molecular dynamic simulation at 205K. In the NVT\nensemble, pressure regulation is achieved by adding different numbers of\nnitrogen molecules (Supplementary Fig. 9). Due to the limited size of the\nsystem box, 20, 50, 100, 200 and 300 nitrogen molecules were added to obtain\napproximate pressures of 1.0, 2.5, 5.5, 11.0 and 17.5 atmospheres, respectively.\nFig. 3a shows that the contact angle of nanodroplets decreases with increasing\ngas pressure, and the droplet curvature decreases first and then increases with\nincreasing gas pressure for energy parameter ε=0.4kcal/mol. Taking pressure\nof 11.0 atm as an example, for 0.1 kcal/mol≤ε<0.43kcal/mol, the liquid-solid\nphase transition only generates 3D ice with disordered stacking of hexagonal\nice and cubic ice (Supplementary Fig. 10). For 0.43 kcal/mol≤ε≤1.0 kcal/mol,\nthe liquid-solid phase transition produces coexisting ice, which is composed of\n2D ice and 3D ice stacked disorderly with hexagonal ice and cubic ice\n(Supplementary Fig. 11). With increasing wetting characteristics, the contact\nangle of droplets decreases, and droplets gradually spread on the wall. After\nquenching, 2D ice begins to form instantaneously, and 3D ice is formed after\novercoming the nucleation energy barrier. The variation of potential energy\nwith time also reflects the occurrence of nucleation and icing process\n(Supplementary Fig. 12). When the surface energy parameters is 0.5 kcal/mol,\n2D ice grows and merges slowly for pressure 11.0 atm (and 17.5 atm) as shown\nin Fig. 3b, which is attributed to the combined effect of surface wettability and\ngas pressure. After calculation of 200ns, the 2D ice is still composed of double-\nlayer 5-, 6-, 7-membered ring water molecules.\ncontact angle and curvature versus pressure. b A sequence of snapshots of the coexisting\nice growth process. c Phase diagram of 3D ice and 2D-3D coexisting ice with respect to\ngas pressure and energy parameter ε.\nThe phase diagram about energy parameters and pressure is predicted to\ndistinguish 2D ice, 3D and 2D-3D coexisting ice by a large number of\nmolecular dynamic simulations. The critical value of the surface energy\nparameter is estimated to be 0.43 kcal/mol, which is the same as that of the\nnitrogen-free system, indicating that the formation of 2D-3D coexisting ice on\nthe platinum surface is not very sensitive to the applied gas pressure conditions.\nThe effects of multiple pressure conditions and multiple surface wetting\ncharacteristics on the liquid-solid phase transition of water molecules are\nrevealed as shown in Fig. 3c. For 1.0 atm≤P≤17.5 atm and ε＜0.43 kcal/mol,\nthe supercooled droplet undergoes liquid-solid phase transition to solid phase,\nand the 3D ice is composed of disordered hexagonal ice and cubic ice. For 1.0\natm≤P≤17.5 atm and 0.43kcal/mol≤ε≤1.0 kcal/mol, liquid-solid phase\ntransition occurs, and the coexisting ice consists of 2D ice and 3D hexagonal\nice and cubic ice. For 0.43 kcal/mol≤ε＜0.8 kcal/mol, the 2D ice in coexisting\nice is composed of double-layer 5-, 6-, 7-membered ring water molecules. For\n0.8 kcal/mol≤ε≤1.0 kcal/mol, the 2D ice in coexisting ice is finally composed\nof a layer of 5-, 6-, 7-membered ring water molecules and a layer of 4-, 6-\n13\nmembered ring water molecules near solid surface (Supplementary Fig. 13).\nThe appearance of the 4-membered square ring water molecules is attributed to\nthe surface being too hydrophilic and gas pressure.\nThe updated figure and the explanation above have been incorporated to the\nrevised manuscript on lines 228-274, 501-508 and Supplementary Information\non lines 85-117.\n5. In their previous work (ref. 44 of the manuscript. Proc. Natl. Acad. Sci. U.S.A. 116,\n16723-16728 (2019)), the authors comment that using NpxyT simulations, the number\nof pentagonal, heptagonal, and octagonal rings in the interior part of the bilayer ice\ndecreases and even goes to zero at longer simulation times. How does this effect affect\nthe coexistence of 2D/3D ice? Since the structure of the 2D ice changes, could the\nadhesion energy of the 2D ice (de-icing mechanism) be modified? This point should be\naddressed in the manuscript as it might have implications in the conclusions of the work.\nReply: Thank you for your constructive question. In the NVT ensemble, when the\ncalculation time reaches 200 ns, the 2D and 3D coexisting ice structures are\nrelatively stable. The 2D ice is still composed of double-layer 5-, 6-, 7-\nmembered ring water molecules, and the proportion of 5-, 7-membered ring\nwater molecules is basically unchanged, indicating that the 2D and 3D\ncoexisting ice is not sensitive to the calculation time when it is sufficiently long\n(Supplementary Fig. 8). The ice adhesion force of the 2D and 3D coexisting ice\nis not affected by the calculation time.\nIn NPxyT ensemble, the pressure is set to Pxy=1atm. With the increase of\ncalculation time, the 5-, 7-membered rings in the double-layer ice are reduced.\nAfter 200ns of calculation, the 2D ice composed of metastable double-layer 5-,\n6-, 7-membered ring water molecules is transformed into steady-state 6-\nmembered ring 2D ice. 3D ice structure is not sensitive to the calculation time.\nSupplementary Fig. 17 2D and 3D coexisting ice structure variation at 205K\nmolecules. At 200ns, the 2D ice is composed of double-layer 6-membered ring\nwater molecules. With the increase of energy parameter ε , the ice adhesion\nw-Pt\nforce at 200ns becomes gradually larger than that at 30 ns, which is attributed\nto the combined effect of surface wetting characteristics, 2D ice structure and\nice freezing time.\nSupplementary Fig. 18 Relationship between the ice adhesion force and wettability\nenergy parameter ε for the coexisting ice at 205K\nw-Pt\nThe updated figure and the explanation above have been incorporated to the\nrevised manuscript on lines 195-199, 364-372 and Supplementary Information\non lines 84-89, 134-148.\n6. Unequivocally distinguishing ice types and water local structure is crucial in\ndescribing nucleation phenomena. The performance of a given order parameter has\nbeen widely discussed in the literature, since it may yield different results. (J. Chem.\nPhys. 144, 034501 (2016)) How do different order parameters affect the formation of\n2D and 3D ice in their simulations? Have the authors used a minimum mislabeling\nprocedure to determine the values used for the q6 and q3 parameters? The influence of\nthese parameters on the results presented in the manuscript should be stated, if any.\nReply: Thank you for your comments.\n16\na) To identify whether there is ice phase in supercooled water, the bond order\nparameters are used for quantitative discrimination based on the identification\nmethod proposed by Steinhardt et al. [62] and developed by Wolde et al. [63].\n1\nq i N b i Y  r  (5)\nlm N i j1 lm ij\nb\nl q iq*  j\nci, j\n ml lm lm (6)\n l q iq* i 12 l q  jq*  j 12\nml lm lm ml lm lm\n \nwhere N (i) is the number of nearest neighbors of particle i, Y r is a\nb lm ij\nspherical harmonic function, r is the unit vector connecting particle i and one\nij\nof the nearest neighbors of particle j, and * denotes the complex conjugate\nfunction. Adopting the CHILL+ algorithm, the c(i, j) with a correlation to q\n3\n(l=3) is used to distinguish crystals (cubic and hexagonal ice) from liquid. The\nnumber of staggered bonds [c(i, j)≤−0.8] and eclipsed bonds [−0.35≤c(i, j)≤\n0.25] are used to identify cubic ice, hexagonal ice and liquid water [64].\nAccording to the reference (J. Chem. Phys. 144, 034501 (2016)), we calculate\nand analyze the effects of local order parameters q (l=3), q (l=4), q (l=6) on\n3 4 6\nliquid water, 2D ice and 3D ice on Pt surface. Fig.1 shows the distribution of\nthe local order parameters q , q , q for 3D ice, 2D ice and liquid water. The\n3 4 6\ndistributions of q (l=3) for liquid water and 2D ice show a large overlap in Fig.\n3\n1a. The 3D ice can be distinguished from 2D ice and liquid water. The\ndistributions of q (l=4) for liquid water, 2D ice and 3D ice show a fairly large\n4\noverlap in Fig. 1b. It is difficult to identify ice molecular from liquid water\nmolecular. In contrast, the distributions of q (l=6) for liquid water and 2D ice\n6\nshow a fairly large overlap. There is a small overlap in the distribution of local\nq parameters between 2D ice and 3D ice. It is difficult to identify 2D ice from\n6\nliquid water. Therefore, we choose CHILL+ algorithm (about q ) to identify\n3\ncubic ice, hexagonal ice and liquid water. In our simulation system, there are\n2430 water molecular. After quenching and cooling, 2D ice grows close to the\n17\nwall, 3D ice growth is based on overcoming the free energy barrier. The\nproportion of 2D ice in coexisting ice is very small. The number of 2D ice\nmolecules is obtained by direct observation based on reference [42]. Previous\nmethods for identifying ice phase from liquid water are mostly based on bulk\nwater and homogeneous nucleation. The identification of 2D ice depended on\nsolid surface growth is a challenging subject, and we need to further explore the\nrelevant local order parameters. It will be our next research direction.\nFig. 1 Distribution of the local order parameters a q for 3D ice, 2D ice and liquid water.\n3\nb q for 3D ice, 2D ice and liquid water. c q for 3D ice, 2D ice and liquid water.\n4 6\nb) We are very sorry for not using a minimum mislabeling procedure to determine\nthe bond order parameter values because we adopt the CHILL+ algorithm [64]\nrelated to bond order parameter.\nThe updated content above have been incorporated to the revised manuscript on\nlines 510-519.\n7. I believe that the manuscript would benefit from a discussion on how the simulations\ncan be translated into experiments. For instance, mentioning that changes in the ε\nparameter could potentially be achieved through modifications of lattice planes or\nsurface characteristics would make the paper more accessible to experimentalists and\nbroaden its audience.\nReply: Thank you for your kind comments. Our work will hopefully motivate future\nexperimental studies of the formation of low- and high-dimensional ice on a\nvariety of realistic surfaces. The change of the ε parameter could potentially be\nachieved through modifications of lattice planes or surface characteristics. The\nadjustment of surface characteristics is achieved by adjusting the wettability or\n18\nconstructing rough structures. 2D ice on the surface can promote or inhibit the\nformation of 3D ice, which has potential application value for the design and\ndevelopment of anti-icing materials. The conditions for 2D-3D coexisting ice\nwithout confinement are surface energy parameter 0.43 kcal/mol＜ε ≤1.0\nw-Pt\nkcal/mol. The contact angle of the water droplet on the platinum metal surface\nis θ<75°, which is in the hydrophilic and superhydrophilic state. It is\nimportant to obtain the hydrophilic and superhydrophilic surface by\nconstructing rough structure and adjusting wettability. There are two main ways\nto prepare hydrophilic and superhydrophilic surface. One is light-induced\nmethod, which uses the photocatalytic process of photocatalytic materials such\nas TiO films to produce hydrophilicity and superhydrophilicity. The second is\n2\nto construct microscopic rough structures on the hydrophilic surface to make it\nmore hydrophilic, by increasing the roughness of the substrate surface.\nThe updated content above have been incorporated to the revised manuscript on\nlines 468-475.\nMinor comments:\n1. There are some sentences that are quite similar in the solid-gas and solid-liquid\ndiscussion. It would be beneficial for the authors to unify the text.\nReply: Thank you for your question. We have revised the section on solid-gas\ndiscussion accordingly.\n2. In Figure 1 (left), the red axis is not explained in either the text or the caption. Please\ninclude an explanation.\nReply: Thank you for your question. We have added the sentence “The mean-square\ndisplacement of ε=0.50 kcal/mol corresponds to the vertical axis on the right-\nhand side” in the revised manuscript.\n3. The citation format needs to be revised and unified.\nReply: Thank you for your question. We have revised and unified the citation format in\nthe manuscript.\nThe authors are using computer simulations to observe the solidification process of\nwater on solid surfaces.\nWhile I think their point of view is interesting, they do not seem to have investigated it\nin sufficient depth. They give various descriptions of the crystallization process, but\nthey are not accompanied by data to corroborate them, so we can only guess them from\nthe context. By the time I got to page 5, I was confronted with various points of\nuncertainty. I think this study (or this paper) is incomplete.\nReply: We sincerely thank the reviewer’s positive comments and useful suggestions. In\nthe following, we further clarify the reviewer’s concerns and indicate revisions\nmade in the revised manuscript and Supplementary Information (SI)\naccordingly.\n1. Line 94: “The potential energy per water molecule at 250 K is significantly higher\nthan that of per water molecule at other temperatures. “: This sentence is unclear; it\nmentions 250 K, but without any data, I have no idea what they are referring to.\nReply: Thank you for your question. The potential energy per water molecule at 250 K\nis significantly higher than that of per water molecule at other temperatures as\nshown in Supplementary Fig. 3.\nSupplementary Fig. 3 Potential energy per water molecule vs temperature\n20\n2. Line 109 Figure 1(b): The scale on the right has no explanation and makes no sense.\nRegarding the plot of 0.5 kcal/mol in red dots, the authors say it shows that liquid and\nsolid coexist, but if they coexist, wouldn't the graph reflect only the liquid character\nwith a large diffusion coefficient?\nReply: Thank you for your constructive question, which help us improve the Fig.1b to\nthe version below. In Fig.1b, the vertical axis on the right-hand side corresponds\nto the mean-square displacement of ε=0.50 kcal/mol. The mean-square\ndisplacement of ε=0.50 kcal/mol is one order of magnitude higher than that of\nε=0.43 kcal/mol and 0.30 kcal/mol. The mean square displacement of ε=0.50\nkcal/mol increases sharply at first, showing there are both liquid phase and 2D\nice solid phase, then the inflection point appears and tends to increase slowly,\nindicating there are both 2D and 3D coexisting ice solid phase.\nFig.1 b Mean-square displacement of water molecules at T=205 K\nThe updated content above have been incorporated to the revised manuscript on\nlines 110-116.\n3. Around line 124: various explanations are given, but no evidence seems to be\npresented. For example,\nReply: Thank you for your question. Around line 124: For ε=0.10 kcal/mol, the\nformation of critical nucleus preferentially occurs inside the droplet rather than\non the solid surface, which is attributed to the superhydrophobic state of the\nsurface (θ>150°). The relevant figure and data are shown in Supplementary Fig.\n21\n2 and Fig. 4. For ε=0.10 kcal/mol, nucleation and growth occurs inside the\nnanodroplet and then diffuses to the outside, which belongs to homogeneous\nnucleation. The contact angle of nanodroplet is 163º, which is larger than 150º,\nso the surface property is superhydrophobic.\nSupplementary Fig. 2 and Fig. 4 CA and Ice dynamic formation process\n4. Line 122: They state that nucleation occurs in the droplet and not on the solid surface\nbecause the surface is superhydrophobic, but they do not show a causal relationship\nbetween the two, which reads as just an unfounded prediction.\nReply: Thank you for your question. For ε=0.10 kcal/mol, the contact angle of\nnanodroplet is 163º (θ>150°) in Supplementary Fig. 2, so the surface property\nis superhydrophobic. After quenching, nucleation occurs in the droplet and not\non the solid surface as shown in Supplementary Fig. 4, which is attributed to\nthe superhydrophobicity of the surface.\n5. Line 128: “With increasing wettability, the contact angle of nanodroplet decreases\nand the contact radius increases, which enhances the interaction between droplet and\nthe surface.” I think the causal relationship is the opposite.\nReply: Thank you for pointing out this issue. We have revised the sentence to “With\nthe increase of wettability, the interaction between nanodroplet and the surface\nenhances, so the contact angle of nanodroplet decreases and the contact radius\nincreases.”\nlines 139-141.\n6. Line 131: “Nanodroplet on surface with higher wettability is more likely to nucleate\nand crystallize due to the very high adsorption energy of the surface” It is just a guess\nand no evidence is provided.\nReply: In our previous work [1], the relationship between the nucleation rate of\nnanodroplets and the surface wettability was investigated. The surface energy\nis the nanodroplet-surface binding energy, which governs the surface\nwettability. The range of ε from 0.115 to 0.346 kcal/mol is able to cover the\nw-Pt\nsurface wettability from hydrophobic to hydrophilic. The hydrophilic surface\nattracts water molecules more easily than the hydrophobic surface owing to the\nstrong solid-water characteristic energy. It is seen that the ice nucleation rate is\nslightly enhanced by increasing the energy parameter ε . The enhancement of\nw-Pt\nice nucleation rate is attributed to the increased ice nucleation on hydrophilic\nsurfaces with very high adsorption energy. The heterogeneous ice nucleation\ncan be promoted due to the nucleation on solid surfaces with very high\ninteraction strength [2].\n3.0\n2.5\n2.0\n1.5\n1.0\n0.5\n0.0\n0.1 0.2 0.3 0.4\nFig.2 The nucleation rate as a function of surface wettability\n[1] Jiang J., Li G.X., Sheng Q., Tang G.H. Microscopic mechanism of ice nucleation:\nThe effects of surface rough structure and wettability. Appl. Surf. Sci. 510, 145520 (2020).\n[2]M. Fitzner, G.C. Sosso, S.J. Cox, A. Michaelides. The many faces of heterogeneous ice\nnucleation: interplay between surface morphology and hydrophobicity. J. Am. Chem. Soc.\n137, 13658–13669 (2015).\n7. Line 134: I don't know what a “large-scale hydrogen bond” is.\n3882.0\nSurface wettability (kcal/mol)\nReply: Apologies for the unclear expression. We have changed “large-scale hydrogen\nbond” to “a large number of hydrogen bonds”. Ice nucleation occurs once a\nsufficient number of relatively long-lived hydrogen bonds develop\nspontaneously at the same location to form a fairly compact initial nucleus. The\ninitial nucleus then slowly changes shape and size until it reaches a stage that\nallows rapid expansion, resulting in crystallization of the entire system [1].\nFig. 3 The hydrogen bond network structure of water [1]\nIt is certainly interesting to discuss the role of hydrogen bond, but we are\nunfortunately not able to do it using mW model. The choice of mW model here\nshould be attributed to its high calculation efficiency, especially for ice-\nnucleation, a very slow process for molecular dynamics simulations. Moreover,\nmW model has been verified for its application in ice-nucleation.\n[1] Matsumoto, M., Saito, S. & Ohmine, I. Molecular dynamics simulation of the ice\nnucleation and growth process leading to water freezing. Nature 416, 409-413 (2002).\n8. Line 135: The phenomenon of molecules aligning to form crystals is not limited to\nwater. I don't think water is particularly prone to ordering, but if you claim that, you\nneed evidence.\nstacking disordered hexagonal ice and cubic ice. We have deleted it.\n9. Line 136: It is an obvious statement over 10 lines.\nReply: Thank you for your question. We have modified it to make it less than 5 lines.\nThe updated content above have been incorporated to the revised manuscript on\nlines 147-151.\n10. Line 143: This is not true. The most stable phase is the phase with the lowest free\nenergy, not the lowest energy.\nReply: Apologies for the unclear expression. Thank you for your question, which help\nus improve the readability of the manuscript. We have rewritten the sentence to\n“From the thermodynamic point of view, the most stable phase is the phase with\nthe lowest free energy, depending on the size of the solid phase free energy and\nliquid phase free energy.”\nThe updated content above have been incorporated to the revised manuscript on\nlines 151-153.\n11. Line 145 is correct, but it is an obvious thermodynamic fact and I am not sure what\nthey are trying to emphasize.\nReply: Thank you for your question. We have deleted the sentence.\n12. Line 147: “The nucleation and crystallization of nanodroplets mainly occur on the\nsolid surface” Since no conditions are set and no evidence is given, it is not clear what\nkind of situation is envisioned in the statement.\nReply: Thank you for your question. The relevant evidence is shown in Supplementary\nFig. 4 and Fig. 5. For ε=0.10 kcal/mol, the nucleation and growth of\nnanodroplets occur inside the nanodroplet and then diffuses to the outside. For\nε=0.20-1.0 kcal/mol, the nucleation and crystallization of nanodroplets occur on\n25\nthe solid surface. Therefore, the nucleation and crystallization of nanodroplets\nmainly occur on the solid surface.\nSupplementary Fig. 4 and Fig. 5 Ice dynamic formation process on surfaces\nReviewer #2 (Remarks to the Author):\nAfter the revisions, the manuscript has undergone significant improvements. I would like to extend\nmy appreciation to the authors for their thorough review and thoughtful incorporation of the\nsuggested changes. Their responsiveness to the required revisions is commendable. Considering\nthe enhancements made, I recommend the publication of this manuscript in Nature\nCommunications.\nReviewer #3 (Remarks to the Author):\nIn general, there are many parts that are written as if they are facts, and it is difficult for the\nreader to discern them. It is necessary to clearly state whether it is a logically derived general fact,\nan experimental fact, or a fact discovered in this study. For example, in line 149, there is a\nstatement that small kinetic energy causes stacking disorder, but this is a causal law that is not at\nall obvious to the reader, and a rationale should be provided.\nThere are mentions of free energy in various places, but they are just rhetoric, since no free\nenergy calculations were actually done in this paper. Low free energy and thermodynamically\nstable are synonymous, and it makes no sense to substitute the latter for the former without a\nnumerical argument. If the authors are going to discuss the free energy barriers, they should have\nan estimate of their height, etc. Otherwise, it is mere speculation. Such expressions are present\nthroughout the paper. (e.g. lines 127, 151, 180, 204)\nFigure 1b Two snapshots are indicated by arrows from the red line, but the meaning is unclear.\nFigure 1b It is difficult to understand why the graph of mean-square displacement, which is the\naverage of many samples, bends at a particular position. The caption reads \"Mean-square\ndisplacement of water molecules to the solid surface at T=205 K.\" It is not clear what is meant by\n\"to the solid surface”. Since the main text says \"on the solid surface,\" I assumed that it refers to\nwater molecules in contact with the surface at time 0, but I still do not understand the mechanism\nby which the diffusion velocities of all the molecules in question slow down at 0.4 ns after the start\nof the measurement. A more detailed explanation of this situation is needed.\nLine 127 In nucleation theory, it is accepted that nucleation proceeds by crossing a free energy\nbarrier, but it is unclear whether nucleation theory can be applied in a finite size system and in a\nspatially anisotropic system. What does free energy mean in such a case? If the authors are going\nto refer to free energy, they need to discuss based on the component values quantitatively.\nLine 146 This expression is unclear and should be presented quantitatively. As far as I understand,\nthe number of hydrogen bonds in the liquid and solid phases should not be significantly different.\nLine 149 The causal relationship between the small kinetic energy and the occurrence of stacking\ndisorders is unknown.\nLine 151 Since no free energy calculations were performed in this paper, what is said here is\nmerely a reiteration of the obvious (as described in textbooks) that when a phase transition\noccurs, the phase that occurs has a lower free energy. I consider it almost meaningless and\nmisleading to mention free energy in a paper that does not perform free energy calculations or\ncomponent analysis of thermodynamic quantities based on thermodynamic equations.\nLine 154 Again, it is unclear whether this is a general fact or a result found in this study. It is also\nnot clear to what degree \"mainly\" represents.\nLine 155 Numerical data is not provided. Without a discussion of how high the energy barrier is\nand, in contrast, how much kinetic energy is required to cross the barrier (in which case, not all\nkinetic energy may be directed toward crossing the barrier), it is impossible to determine whether\nthe statement is even correct. If such an argument has been made in previous studies, it should\nbe cited. If it is a fact found in this study, it should be written as such and accompanied by\nnumerical evidence.\nLine 160 The definition of induction time varies from paper to paper. If the author defines\ninduction time as the time it takes for the potential energy to begin to decrease, then it should be\nwritten as such. However, there is no mention of induction time in the text hereafter, and if only\none simulation was performed for each condition, it is not possible to discuss the statistical trend\nof induction time, so I think this line is needless.\nLine 172 “Non-rotation stacking” is not a commonly used term. What is referring to rotating in\neither direction? What rotates? A water molecule or the whole two-dimensional ice? If it is already\ndefined in another paper, cite it.\nLine 199 What is written in this paragraph cannot be read from Figure 2b. If \"intensity\" refers to\nthe area of the peak, then adding an integrated area plot would make the claim more credible.\nLine 201 It needs to be explained how one can read from Figure 2b that ice growth begins at the\nbottom.\nFigure 2 Water that was droplet in the liquid state crystallizes as bilayer ice covering the ice\nsurface, and after the bilayer ice can no longer grow, the droplet begins to crystallize. In this case,\nthe ratio of the amount of water in the bilayer ice to that in the droplet should depend on the\nsubstrate surface area, and it is unlikely that what happens in the paper will always happen. I feel\nthat the paper describes a phenomenon that occurs in one particular situation as if it can be\ngeneralized.\nFig. 2a No explanation of color coding.\nFig. 2d There is no legend for the marks.\nFigure 2 shows that the droplet itself moves around by thermal fluctuations. This motion\ncorresponds to the translation of the entire system in terms of the bulk system and should be\nexcluded from the calculation of the diffusion coefficient. A detailed explanation of how to calculate\nthe mean-square displacement is needed.\nResponse Sheet\nTitle: Two-dimensional bilayer ice in coexistence with three-\ndimensional ice without confinement\nPaper number: NCOMMS-23-05786A\nWe are glad to see that reviewers have recommended our manuscript for publication.\nWe would like to thank the reviewers for their constructive suggestions. In the following,\nwe have responded point-by-point to all their comments and questions and indicated\nthe corresponding changes made in the revised manuscript. All the changes are marked\nin blue in the revised manuscript.\nReviewer #2 (Remarks to the Author):\nAfter the revisions, the manuscript has undergone significant improvements. I would\nlike to extend my appreciation to the authors for their thorough review and thoughtful\nincorporation of the suggested changes. Their responsiveness to the required revisions\nis commendable. Considering the enhancements made, I recommend the publication of\nthis manuscript in Nature Communications.\nReply: We are grateful for the reviewer’s efforts to improve the quality of our paper.\nIn general, there are many parts that are written as if they are facts, and it is difficult for\nthe reader to discern them. It is necessary to clearly state whether it is a logically derived\ngeneral fact, an experimental fact, or a fact discovered in this study.\nReply: We thank the reviewer’s positive comments. In the following, we have\naddressed the reviewer’s concerns and indicated our revisions in the revised\nmanuscript and Supplementary Information (SI) accordingly.\n1. For example, in line 149, there is a statement that small kinetic energy causes stacking\ndisorder, but this is a causal law that is not at all obvious to the reader, and a rationale\nshould be provided.\nReply: When the kinetic energy of the molecule is small, the potential energy between\nthe molecules is sufficient to confine the molecules, so that the molecules are\nusually arranged in a regular crystal structure. Simulations with various water\nmodels find that ice nucleated and grown is at all sizes stacking-disordered,\nconsisting of random sequences of cubic and hexagonal ice layers [1-5]. Our\nsimulation results reveal that the three-dimensional ice crystal structure consists\nof stacking-disordered hexagonal ice and cubic ice. We have revised the\nsentence to “When the kinetic energy of the molecule is small, the potential\nenergy between the molecules is sufficient to confine the molecules, so that the\nmolecules are usually arranged in a regular crystal structure.”\n[1] Haji-Akbari, A. & Debenedetti, P. G. Direct calculation of ice homogeneous nucleation\nrate for a molecular model of water. Proc. Natl Acad. Sci. USA 112, 10582-10588 (2015).\n[2] Li, T., Donadio, D., Russo, G. & Galli, G. Homogeneous ice nucleation from\nsupercooled water. Phys. Chem. Chem. Phys. 13, 19807-19813 (2011).\n[3] Moore, E. B. & Molinero, V. Is it cubic? Ice crystallization from deeply supercooled\nwater. Phys. Chem. Chem. Phys. 13, 20008-20016 (2011).\n[4] Malkin, T. L. et al. Stacking disorder in ice I. Phys. Chem. Chem. Phys. 17, 60-76 (2015).\n[5] Reinhardt, A. & Doye, J. P. K. Free energy landscapes for homogeneous nucleation of\nice for a monatomic water model. J. Chem. Phys. 136, 054501 (2012).\n155-157.\n2. There are mentions of free energy in various places, but they are just rhetoric, since\nno free energy calculations were actually done in this paper. Low free energy and\nthermodynamically stable are synonymous, and it makes no sense to substitute the latter\nfor the former without a numerical argument. If the authors are going to discuss the free\nenergy barriers, they should have an estimate of their height, etc. Otherwise, it is mere\nspeculation. Such expressions are present throughout the paper. (e.g. lines 127, 151,\n180, 204)\nReply: Thank you for your question. Nucleation is the process with which the formation\nof new phases begins and is thus a widely spread phenomenon in both nature\nand technology. Classical nucleation theory (CNT) was formulated 90 years ago\nthrough the contributions of Volmer and Weber, Farkas, Becker and Dö ring,\nand Zeldovich, on the basis of the pioneering ideas of none other than Gibbs\nhimself. CNT was created to describe the condensation of supersaturated vapors\ninto the liquid phase, but most of the concepts can also be applied to the\ncrystallization of supercooled liquids and supersaturated solutions. Clusters of\ncrystalline atoms occur within the supercooled liquid by spontaneous,\ninfrequent fluctuations, which eventually lead the system to overcome the free\nenergy barrier for nucleation triggering the actual crystal growth [1] as shown\nin Fig.1.\nFig.1 Sketch of the free energy difference ΔG , as a function of the crystalline nucleus size\nN\nn [1]. A free energy barrier for nucleation, ΔG* , must be overcome to proceed from the\nN\n3\n(metastable) supercooled liquid state to the thermodynamically stable crystalline phase\nthrough homogeneous nucleation (purple). Heterogeneous nucleation (green) can be\ncharacterized by a lower free energy barrier, ΔG* ,het, and a smaller critical nucleus size,\nN\nn *.\nhet\nThe existence of critical ice nucleus was first confirmed by experimental study.\nThe experimental results are in good agreement with the predictions of the\nclassical nucleation theory (CNT) on the critical nucleus and free energy. It is\nthe first time to experimentally confirm the existence of the critical ice nucleus\nand the dependence of its size on the supercooled temperature during water\nfreezing [2]. This study effectively clarifies the general doubts in recent decades\nabout the effectiveness of the classical nucleation theory (CNT) in describing\nthe critical nuclear characteristics at the atomic scale, and deepens our\nunderstanding of the microscopic mechanism for the important phase transition\nduring water freezing.\nA basic physical understanding of homogeneous and hetegeneous ice nucleation\nis provided by the classical nucleation theory (CNT). The CNT suggests that\nwater molecules must form an ice nucleus of critical size (critical ice nucleus)\nbefore a water crystallization process occurs. The formation of the critical ice\nnucleus is associated with a free energy barrier, which needs to be overcome to\ntrigger the further growth and crystallization of ice.\nThe research objectives of this paper are liquid-solid, gas-solid phase diagram\nprediction and de-icing mechanism analysis. The free energy barrier is\ncalculated using the kinetic reconstruction method [3] during the ice nucleation\nprocess. This method requires two quantities: the steady-state probability\ndistribution P (n), which represents the probability of finding the system in a\nst\nstate with the largest ice nucleus having n ice molecules, and the mean first-\npassage time (MFPT), τ(n), which is defined as the average time needed for the\nlargest nucleus to contain n molecules for the first time [4]. This method has\nbeen used to study the free energy barrier for ice nucleation [5, 6].\ninitial conditions were performed by MD simulations. For each simulation, the\nsize of the largest cluster in the system is continuously monitored at regular\nintervals (every 1000 time steps), and the corresponding first appearance time\nfor each size n, t(n), is recorded. The mean first-passage time τ(n) for each size\ni\nn is simply obtained by averaging t(n) over multiple repetitions of simulation\ni\nwith different initial configurations. Additionally, we obtain the steady-state\nprobability P (n): at every 1000 time steps during each repetition, we sample\nst\nthe size n of the largest cluster in the system. Subsequently, a histogram is\nconstructed by counting how frequently a given size n occurs as the largest\ncluster. By normalizing this count with respect to the total number of sampled\ncluster sizes, we derive the steady-state probability distribution P (n).\nst\nSupplementary Fig. 6c Free energy as a function of the cluster size n (β=1/k T).\nB\nSupplementary Fig. 6c depicts βΔG as a function of the number of ice molecules\nn in the largest nucleus for this case. It can be observed that βΔG increases with\nincreasing nucleus size after the inception of a small ice nucleus. When the\nnucleus reaches its critical size, n=23, βΔG approaches its maximum value.\nSubsequently, as the ice nucleus continues to grow, the free energy decreases\nand triggers further ice formation. The value ΔG* = 4.67k T represents the\nB\nenergy difference between the maximum free energy and the minimum free\nenergy after inception of an ice nucleus. ΔG* = 4.67k T is associated with the\nB\ndevelopment of the ice nucleus on the solid wall.\n5\n[1] Gabriele C., et al. Crystal Nucleation in Liquids: Open Questions and Future Challenges\nin Molecular Dynamics Simulations. Chem. Rev. 116, 7078-7116 (2016).\n[2] Bai G.Y., Gao D., Liu Z., Zhou X., Wang J.J. Probing the critical nucleus size for ice\nformation with graphene oxide nanosheets. Nature, 576, 437-441(2019).\n[3] Wedekind J., Reguera D. Kinetic reconstruction of the free-energy landscape. J. Phys.\nChem. B 112, 11060-11063 (2008).\n[4] Wedekind J., Strey R., Reguera D. New method to analyse simulations of activated\nprocesses. J. Chem.Phys. 126, 134103 (2007).\n[5] Luo S., Wang J., Li Z.G. Homogeneous ice nucleation under shear. J. Phys. Chem. B\n124, 3701-3708 (2020).\n[6] Li C., et al. Enhancing and impeding heterogeneous ice nucleation through nanogrooves.\nJ. Phys. Chem. C 122, 25992-25998 (2018).\nThe updated figure and the explanation above have been incorporated to the\nrevised manuscript on line 134 and Supplementary Information on lines 68-96.\n3. Figure 1b Two snapshots are indicated by arrows from the red line, but the meaning\nis unclear.\nReply: Thank you for your kind question. In Fig.1b, the vertical axis on the right-hand\nside corresponds to the mean-square displacement of ε=0.50 kcal/mol. The\nmean-square displacement of ε=0.50 kcal/mol is one order of magnitude higher\nthan that of ε=0.43 kcal/mol and 0.30 kcal/mol. The mean square displacement\nof ε=0.50 kcal/mol increases sharply at first, showing some water molecules of\nthe droplet diffuse rapidly along the wall to form 2D ice. In this process, the\nliquid water exists at the same time as the 2D ice solid. Then the inflection point\nappears and tends to increase slowly, indicating that the liquid water molecules\nnucleate and grow into 3D ice. At this time, both 2D ice solid phase and 3D ice\nsolid phase are present. The whole freezing process has two phase transitions.\nThe first phase transition produces 2D ice, and the second phase transition\nproduces 3D ice. The appearance of the inflection point means that the liquid\nwater molecules nucleate and crystallize in an instant to form 3D ice, and only\nthe solid state exists.\nlines 113-121.\n4. Figure 1b It is difficult to understand why the graph of mean-square displacement,\nwhich is the average of many samples, bends at a particular position. The caption reads\n\"Mean-square displacement of water molecules to the solid surface at T=205 K.\" It is\nnot clear what is meant by \"to the solid surface”. Since the main text says \"on the solid\nsurface,\" I assumed that it refers to water molecules in contact with the surface at time\n0, but I still do not understand the mechanism by which the diffusion velocities of all\nthe molecules in question slow down at 0.4 ns after the start of the measurement. A\nmore detailed explanation of this situation is needed.\nReply: Apologies for the unclear expression. We have revised the sentence to \"Mean-\nsquare displacement of water molecules on the solid surface at T=205 K.\" The\nmean square displacement of ε=0.50 kcal/mol increases sharply at first, showing\nsome water molecules of the droplet diffuse rapidly along the wall to form 2D\nice. In this process, the liquid water exists at the same time as the 2D ice solid.\nThen the inflection point appears and tends to increase slowly, indicating that\nthe liquid water molecules nucleate and grow into 3D ice. At this time, there are\nboth 2D ice solid phase and 3D ice solid phase. The whole freezing process has\ntwo phase transitions. The first phase transition produces 2D ice, and the second\nphase transition produces 3D ice. The appearance of the inflection point means\nthat the liquid water molecules nucleate and crystallize in an instant to form 3D\nice, and only the solid state exists.\nThe updated content above have been incorporated to the revised manuscript on\nlines 113-121 and 128-129.\n5. Line 127 In nucleation theory, it is accepted that nucleation proceeds by crossing a\nfree energy barrier, but it is unclear whether nucleation theory can be applied in a finite\nsize system and in a spatially anisotropic system. What does free energy mean in such\n7\na case? If the authors are going to refer to free energy, they need to discuss based on\nthe component values quantitatively.\nReply: Thank you for your question. According to classical nucleation theory, the free\nenergy barrier is used to evaluate nucleation crystallization, which is the most\nconventional method. Some literatures have calculated and analyzed the free\nenergy in a finite size system and in a spatially anisotropic system [1-4].\nDue to the exceedingly small time and length scales involved in early stages of\nnucleation, atomistic computer simulations can provide unique insights into\nmicroscopic aspects of nucleation and crystallization. The dynamics over the\nfree-energy landscape dictate time evolution of systems out of equilibrium. In\nmany systems such as chemical reactions, phase transitions, conformational\nchanges in biomolecules etc., this dynamics is controlled by a free-energy\nbarrier that must be surmounted; these processes are called activated and play a\ncrucial role in many areas of science and technology.\nThe evaluation of activation barriers and free-energy landscapes from\nsimulations can be accomplished using a variety of equilibrium techniques, such\nas umbrella sampling or transition path sampling. However, it should be noted\nthat these equilibrium techniques may not always suffice when studying\nprocesses that occur in an out-of-equilibrium fashion. There is no guarantee\nbeforehand that the landscape reconstructed in an equilibrium simulation truly\nreflects the governing landscape of the real nonequilibrium dynamics involved\nin the process. The introduction of additional constraints in an equilibrium\nsimulation or the application of bias or forces in a nonequilibrium approach can\nsignificantly alter the conditions under which the process takes place and thus\ndistort the true energy landscape. Therefore, it is highly desirable to have a\nmethod for directly evaluating this landscape based on the actual kinetic\nevolution of the system without imposing any artificial external forces or\nconstraints.\nenergy landscape from out-of-equilibrium kinetics observed during molecular\ndynamics (MD) simulations, Brownian dynamics (BD), stochastic simulations,\nand even direct experiments. This technique only requires two key components:\nsteady-state probability distribution and mean first-passage time. Both can be\ncalculated through simulations as well as experimental measurements.\nThe free energy barrier ΔG(n) for ice nucleation is obtained adopting the method\ndeveloped by Wedekind and Reguera [5]. First, the number of ice molecules is\nobtained based on ice molecule identification using the average bond order\nparameter. Here, q3 is applied to distinguish the molecular structure of cubic\nand hexagonal ice from liquid water. The number of ice molecules in the largest\nice cluster, n, is then used to denote the state of the system. To obtain the free\nenergy, two quantities are required in this method. One is the steady probability\nP (n) of the system assuming a state characterized by n. The other is the mean\nst\nfirst passage time (MFPT), τ(n), which is the mean time the system takes to\nreach the state for the first time. Both P (n) and τ(n) can be obtained in MD\nst\nsimulations and the free energy is then calculated through\nn\ndn\nGnlnBn C (1)\n \nwhere β = 1/(k T) with k is the Boltzmann constant, C is a constant, and B(n)\nB B\nis given by\nBn 1   n P ndn n  (2)\nP st n  1 st n 1  \nwhere n is the maximum value of n, which is set as 200. The dependence of\nn  0  1erf  nn* c  (3)\n2\n9\nwhere τ is a characteristic time, c is the local curvature around the top of the\n0\nbarrier, and n* is the critical nucleus size.\nTo obtain these quantities, more than 60 independent simulations with different\ninitial conditions were performed by MD simulations. For each simulation, the\nsize of the largest cluster in the system is continuously monitored at regular\nintervals (every 1000 time steps), and the corresponding first appearance time\nfor each size n, t(n), is recorded. The mean first-passage time τ(n) for each size\ni\nn is simply obtained by averaging t(n) over multiple repetitions of simulation\ni\nwith different initial configurations. Additionally, we obtain the steady-state\nprobability P (n): at every 1000 time steps during each repetition, we sample\nst\nthe size n of the largest cluster in the system. Subsequently, a histogram is\nconstructed by counting how frequently a given size n occurs as the largest\ncluster. By normalizing this count with respect to the total number of sampled\ncluster sizes, we derive the steady-state probability distribution P (n).\nst\nSupplementary Fig. 6 The calculation of ice nucleation free energy barrier at T=205K\n(ε=0.4kcal/mol). a Mean first-passage time (MFPT) τ(n) as a function of the cluster size n\nobtained from the MD simulations. The red line represents a fit of Eq.3. b Growth curves\nof the largest cluster from the MD simulations. c Kinetic reconstruction of the free energy\nof cluster formation obtained from the MD simulations (β=1/k T).\nB\nFig.6a shows that the MFPT curve has a sigmoidal shape that reaches a well-\ndefined plateau at larger values of n, indicating that overcoming critical size is\nrate-limiting step in cluster formation and subsequent growth of cluster takes\nnegligible time compared to activation time. Fig.6b illustrates the growth curves\n10\nof the largest ice cluster from the MD simulations. The free-energy barrier ΔG\nare calculated through Eqs.1 and 2. Fig. 6c depicts βΔG as a function of the\nnumber of ice molecules n in the largest nucleus for this case (ε=0.4kcal/mol).\nIt can be obversed that βΔG increases with increasing nucleus size after the\ninception of a small ice nucleus. It approaches its maximum when the nucleus\nreaches the critical size, n=23. The free energy then decreases with increasing\nice nucleus and triggers further ice formation. The value ΔG* = 4.67k T\nB\nrepresents the energy difference between the maximum free energy and the\nminimum free energy after inception of an ice nucleus. ΔG* = 4.67k T is\nB\nassociated with the development of the ice nucleus on the solid wall.\n[1] Qiu Y.Q., et al. Ice nucleation efficiency of hydroxylated organic surfaces is controlled\nby their structural fluctuations and mismatch to ice. J. Am. Chem. Soc. 139, 3052-3064\n(2017).\n[2] Bi, Y.F., Cao, B.X. & Li, T.S. Enhanced heterogeneous ice nucleation by special surface\ngeometry. Nat. Commun. 8, 15372 (2017).\n[3] Luo S., Wang J., Li Z.G. Homogeneous ice nucleation under shear. J. Phys. Chem. B\n124, 3701-3708 (2020).\n[4] Li C., et al. Enhancing and impeding heterogeneous ice nucleation through nanogrooves.\nJ. Phys. Chem. C 122, 25992-25998 (2018).\n[5] Wedekind J., Reguera D. Kinetic reconstruction of the free-energy landscape. J. Phys.\nChem. B 112, 11060-11063 (2008).\n[6] Wedekind J., Strey R., Reguera D. New method to analyse simulations of activated\nprocesses. J. Chem.Phys. 126, 134103 (2007).\nThe updated figure and the explanation above have been incorporated to the\nrevised manuscript on line 134 and Supplementary Information on lines 68-96.\n6. Line 146 This expression is unclear and should be presented quantitatively. As far as\nI understand, the number of hydrogen bonds in the liquid and solid phases should not\nbe significantly different.\nReply: Apologies for the unclear expression. We have changed “a large number of\nhydrogen bonds” to “a relatively long-lived hydrogen-bonded network”. Ice\nnucleation occurs once a sufficient number of relatively long-lived hydrogen\n11\nbonds develop spontaneously at the same location to form a fairly compact\ninitial nucleus. The initial nucleus then slowly changes shape and size until it\nreaches a stage that allows rapid expansion, resulting in crystallization of the\nentire system [1]. The ice crystal has a long-range ordered structure, while water\ndoes not have a long-range ordered structure. Only some water molecules\nmaintain short-range order with neighboring water molecules. This orderly\n“water cluster” is constantly generating, destroying and moving, and its degree\nis related to temperature. The fluidity of water molecules in liquid water is too\nstrong to form a large number of stable hydrogen bonds like ice. We note that\nthe enthalpy of fusion of ice at atmospheric pressure is 6.01 kJ/mol, whilst the\nenthalpy of vaporisation at the triple point is 51.059 kJ/mol; in the latter process,\nhydrogen bonds must clearly be broken, but the low enthalpy offusion suggests\nthat the majority of hydrogen bonds survive in the liquid phase, although\nperhaps on a transient timescale [2].\nFig. 3 Number of water molecules having long-lasting hydrogen bonds and the\nsize of the largest cluster in inherent structures [1]\nIt is certainly interesting to discuss the role of hydrogen bond, but we are\nunfortunately not able to do it using mW model. The choice of mW model here\nshould be attributed to its high calculation efficiency, especially for ice-\nnucleation, a very slow process for molecular dynamics simulations.\n[1] Matsumoto, M., Saito, S. & Ohmine, I. Molecular dynamics simulation of the ice\nnucleation and growth process leading to water freezing. Nature 416, 409-413 (2002).\n[2] Debenedetti, P.G. ‘Supercooled and glassy water’. J Phys: Condens Mat 15, R1669-\nR1726 (2003).\nThe updated content have been incorporated to the revised manuscript on lines\n152-153.\n7. Line 149 The causal relationship between the small kinetic energy and the occurrence\nof stacking disorders is unknown.\nReply: Apologies for the unclear expression. We have rewritten the sentence to “When\nthe kinetic energy of the molecule is small, the potential energy between the\nmolecules is sufficient to confine the molecules, so that the molecules are\nusually arranged in a regular crystal structure.”\nThe updated content above have been incorporated to the revised manuscript on\nlines 155-157.\n8. Line 151 Since no free energy calculations were performed in this paper, what is said\nhere is merely a reiteration of the obvious (as described in textbooks) that when a phase\ntransition occurs, the phase that occurs has a lower free energy. I consider it almost\nmeaningless and misleading to mention free energy in a paper that does not perform\nfree energy calculations or component analysis of thermodynamic quantities based on\nthermodynamic equations.\nReply: Thank you for your kind question. We have deleted the sentence.\n9. Line 154 Again, it is unclear whether this is a general fact or a result found in this\nstudy. It is also not clear to what degree \"mainly\" represents.\nReply: Thank you for your question. This is a result found in our simulation study. The\nrelevant evidence is shown in Supplementary Fig. 4 and Fig. 5. For ε=0.10\nkcal/mol, the nucleation and growth of nanodroplets occur inside the\nnanodroplet and then diffuses to the outside. For ε=0.20-1.0 kcal/mol, the\nnucleation and crystallization of nanodroplets occur on the solid surface.\nTherefore, the nucleation and crystallization of nanodroplets mainly occur on\n13\nthe solid surface. \"mainly\" represents that the nucleation and crystallization of\nnanodroplets occur on the solid surface of ε=0.20-1.0 kcal/mol.\nSupplementary Fig. 4 and Fig. 5 Ice dynamic formation process on surfaces\nThe updated content above have been incorporated to the revised manuscript on line\n158.\n10. Line 155 Numerical data is not provided. Without a discussion of how high the\nenergy barrier is and, in contrast, how much kinetic energy is required to cross the\nbarrier (in which case, not all kinetic energy may be directed toward crossing the\nbarrier), it is impossible to determine whether the statement is even correct. If such an\nargument has been made in previous studies, it should be cited. If it is a fact found in\nthis study, it should be written as such and accompanied by numerical evidence.\nReply: We have revised the sentence to “the latent heat released by local crystallization\ncan provide enough energy for growth and crystallization [53].” The relevant\nreference [53] have been cited to reinforce the conclusions.\nThe updated content above have been incorporated to the revised manuscript on\nlines 159-160.\n11. Line 160 The definition of induction time varies from paper to paper. If the author\ndefines induction time as the time it takes for the potential energy to begin to decrease,\nthen it should be written as such. However, there is no mention of induction time in the\ntext hereafter, and if only one simulation was performed for each condition, it is not\npossible to discuss the statistical trend of induction time, so I think this line is needless.\nReply: Thank you for your kind question. We have deleted the sentence.\n14\n12. Line 172 “Non-rotation stacking” is not a commonly used term. What is referring\nto rotating in either direction? What rotates? A water molecule or the whole two-\ndimensional ice? If it is already defined in another paper, cite it.\nReply: Apologies for the unclear expression. Thank you for your question, which help\nus improve the readability of the manuscript. “Non-rotation stacking” means\nthat the positions of water molecules in the upper and lower layers of 2D ice are\ncompletely coincident. We have revised the sentence to “The 2D ice is formed\nby two layers of interlocking 5-, 6-, 7- membered rings water molecules.”\nThe updated content above have been incorporated to the revised manuscript on\nlines 175-176.\n13. Line 199 What is written in this paragraph cannot be read from Figure 2b. If\n\"intensity\" refers to the area of the peak, then adding an integrated area plot would make\nthe claim more credible.\nReply: Apologies for the unclear expression. In Fig.2b, “intensity” represents the size\nof density. “Height” is the size of z axis. The intensity of the first peak is slightly\nlarger than that of the second peak, indicating that the growth of 2D ice starts\nfrom the bottom layer (near the solid wall).\nThe updated content above have been incorporated to the revised manuscript on\nlines 202-205 and 214.\n14. Line 201 It needs to be explained how one can read from Figure 2b that ice growth\nbegins at the bottom.\nReply: Apologies for the unclear expression. We have revised the sentence to “The\nintensity of the first peak is slightly larger than that of the second peak as shown\nin Fig. 2b, showing that the growth of 2D ice starts from the bottom layer (near\nthe solid wall), which is consistent with the results in Reference [40]”.\nThe updated content above have been incorporated to the revised manuscript on\nlines 202-205.\n15\n15. Figure 2 Water that was droplet in the liquid state crystallizes as bilayer ice covering\nthe ice surface, and after the bilayer ice can no longer grow, the droplet begins to\ncrystallize. In this case, the ratio of the amount of water in the bilayer ice to that in the\ndroplet should depend on the substrate surface area, and it is unlikely that what happens\nin the paper will always happen. I feel that the paper describes a phenomenon that\noccurs in one particular situation as if it can be generalized.\nReply: Thank you for your question. We have studied many working conditions, found\nthis new phenomenon occurring under all the conditions. At present, the\nknowledge on 2D-3D coexisting ice formation and de-icing characteristics at\nmolecular level remains blank. The influences of the surface wettability,\ntemperature, pressure and phase transition mode on the formation of 2D-3D\ncoexisting ice without confinement has not been studied. Our findings fill the\ngap about the coexistence of 2D bilayer ice and 3D ice without nanoscale\nconfinement. We explore the possibility of controlling the formation of low-\nand high-dimensional ices by tuning the interactions of a water droplet on a\nsurface. Our work will hopefully motivate future experimental studies of the\nformation of low- and high-dimensional ice on realistic surfaces.\n16. Fig. 2a No explanation of color coding.\nReply: Thank you for your kind question. We have added the sentence “Dark yellow,\nlight blue, light gray (3.08ns) and light red balls represent hexagonal ice, cubic\nice, 2D ice and solid substrate, respectively”.\nThe updated content above have been incorporated to the revised manuscript on\nlines 212-213, 261-263 and 446-447.\n17. Fig. 2d There is no legend for the marks.\nReply: Thank you for your kind question. We have added the sentence “Phase diagram\nconsists of liquid water, 2D ice in coexistence with liquid water, 3D ice and\n2D-3D coexisting ice.”\nlines 216-217, 264 and 449-451.\n18. Figure 2 shows that the droplet itself moves around by thermal fluctuations. This\nmotion corresponds to the translation of the entire system in terms of the bulk system\nand should be excluded from the calculation of the diffusion coefficient. A detailed\nexplanation of how to calculate the mean-square displacement is needed.\nReply: Thank you for your question. Our research ignores the influence of thermal\nfluctuation on the mean square displacement, and we refer to the relevant\nliterature [1, 2] for calculation and analysis. The mean square displacement\n(MSD) at time t is calculated by Eq. (5):\n1 n 2\nMSDt  rt r0 (5)\nn i i\ni0\nwhere rt and r0 are the positions of atom i at time t and 0, respectively.\ni i\n[1] Zhu, C.Q. et al. Direct observation of 2-dimensional ices on different surfaces near room\ntemperatature without confinement. Proc. Natl. Acad. Sci. USA 116, 16723-16728 (2019).\n[2] Luo S., Wang J., Li Z.G. Homogeneous ice nucleation under shear. J. Phys. Chem. B\n124, 3701-3708 (2020).\nThe updated content above have been incorporated to the revised manuscript on\nlines 511-514.\nReviewer #2 (Remarks to the Author):\nThe manuscript has seen considerable enhancements as a result of the revisions made by the\nauthors in response to the reviewer's feedback. The authors have diligently reviewed every\ncomment and suggestion put forth by the reviewer. These modifications have greatly enhanced the\noriginal work. With these improvements in mind, I enthusiastically endorse the publication of this\nmanuscript in Nature Communications."}]

Decision letter