Effect of sodium chloride adsorption on the surface premelting of ice

Abstract

We characterise the structural properties of the quasi-liquid layer (QLL) at two low-index ice surfaces in the presence of sodium chloride (Na⁺/Cl⁻) ions by molecular dynamics simulations. We find that the presence of a high surface density of Na⁺/Cl⁻ pairs changes the surface melting behaviour from step-wise to gradual melting. The ions lead to an overall increase of the thickness and the disorder of the QLL, and to a low-temperature roughening transition of the air–ice interface. The local molecular structure of the QLL is similar to that of liquid water, and the differences between the basal and primary prismatic surface are attenuated by the presence of Na⁺/Cl⁻ pairs. These changes modify the crystal growth rates of different facets and the solvation environment at the surface of sea-water ice with a potential impact on light scattering and environmental chemical reactions.

---

1 Introduction

Ice covers 10% of the Earth and its surface is an important catalyst for atmospheric chemical reactions, such as the adsorption, photodegradation, and release of trace gases in snow-packs,¹ the formation of ozone-depleting species in polar stratospheric clouds,² and chemical reactions in the troposphere.^3–5 From the melting point down to about −70 °C, the surface of ice features a disordered pre-melted layer known as the quasi-liquid layer (QLL). Its thickness depends on the temperature and on the orientation of the ice facets.⁶ The molecular structure of the QLL dictates the morphology of growing ice crystals, it controls the adsorption of molecules at the air–ice interface, and the thermodynamics and kinetics of the chemical reactions involving these molecules.

The thickness, structure and diffusivity of the QLL at pure ice surfaces have been extensively studied by both experiments and computer simulations. There are still substantial discrepancies in the estimate of both the onset temperature of surface melting and the thickness of the QLL, which depend on both the experimental conditions and techniques.^6–13 Molecular dynamics (MD) simulations provide a consistent picture of very few (two or three) molten (bi)layers, corresponding to a thickness of about 1 nm, up to 2 K below the melting point.^14–18 Simulations suggest bilayer-by-bilayer “rounded” melting transitions that correlate with the trends observed in sum-frequency generation (SFG) spectra.^9,19 Furthermore, MD shows that the underlying crystalline layers induce an ordered template in the QLL structure that persists up to 2 K below the melting temperature. This dynamic ordering marks the major differences between the QLL and bulk water,^16,20–22 and it affects the solvation properties of the premelting layer.²³

Whereas the studies of pure-water ice surfaces still leave unanswered questions, even more remain about how the molecular and mesoscale structure of the premelting layer changes in the presence of adsorbed species. In natural environments, ice surfaces are exposed to extrinsic species of various kinds, including ions, inorganic, and organic compounds, which impact the structure of the air–ice interface.^24,25 Sodium and chloride are the most abundant ions in seawater and can be found naturally in the atmosphere. Sea salt aerosols deposit on environmental snow and ice, and frozen seawater has a layer of brine at its surface containing a high surface density of Na⁺ and Cl⁻ ions. Glancing-angle Raman spectroscopy suggested that the surface of seawater ice is markedly different from that of freshwater ice.²⁴ According to these measurements, the liquid brine formed at the air–ice interface of seawater ice provides a chemical environment similar to that of bulk supercooled water, possibly leading to changes in the rates of photolysis reactions.²⁶ However, the molecular mechanisms that lead to such a stark difference between pristine and briny ice surfaces and the details of the surface of seawater ice are still mostly unexplored.

In this work, we investigate the effect of adsorbed sodium and chloride ions on the structure of the premelting layer at low-index surfaces of hexagonal ice I_h by means of classical MD simulations. We analyse the molecular structure of the QLL in the presence of Na⁺ and Cl⁻ as a function of temperature and surface orientation. We consider the two-bilayer low-index basal (0001) and primary prismatic (101̄0) surfaces, and we model two naturally occurring Na⁺/Cl⁻ surface densities of 0.1 NaCl pair nm⁻² and 1 NaCl pair nm⁻² The former is similar to the surface density at the surface of snow atop of first-year sea ice, while the latter is similar to the surface density of brine atop of sea ice.²⁷ Our simulations highlight how the ions affect the premelting behaviour compared to pristine ice surfaces. Specifically, the presence of ion pairs at high surface density turns the discrete bilayer-by-bilayer melting, which was observed for pristine ice basal and primary prismatic surfaces, into gradual surface melting. Furthermore, solvated Na⁺/Cl⁻ pairs alter the surface roughness and its fluctuations thus impacting ice crystal growth, the adsorption of other molecules and the chemical processes at the briny ice surface.

2 Computational methods

Studying the premelting of ice surfaces by MD requires the simulation of systems of several thousands of water molecules for time scales of hundreds of nanoseconds. The need for relatively large models was made clear in former works that address finite size effects in the characteristic behaviour of premelting transitions.^15,17 For these reasons, we perform classical MD simulations using the TIP4P/2005 empirical rigid-molecule point-charge model.²⁸ This model, parameterised on a set of experimental properties of water among which the temperature of maximum density, reproduces qualitatively the phase diagram of molecular ice and water,^29,30 and it is among the best in the class of rigid non-polarisable forcefields at reproducing the thermodynamic properties of water and ice.³¹ However, it is important to note that TIP4P/2005 has a melting point of 249.5 ± 0.1 K,³² compared to the experimental value of 273.15 K. To describe the interactions between the water molecules and the Na⁺ and Cl⁻ ions, we employ the Madrid-2019 forcefield,^33,34 which is a point-charge model based on the method of scaled charges. Its parameters are calibrated using hydration numbers, radial distribution functions, and the density of both aqueous solutions and solid NaCl. This forcefield has been used to study the seawater/ice interface, the physical properties of seawater, and the effect of NaCl segregation in ice grain boundaries.^35,36

MD simulations were carried out using the GROMACS 2020.3 program compiled in double-precision.^37,38 The MD equations of motion were integrated with a time step of 1 fs, which guarantees the conservation of total energy in a microcanonical test run with a total energy drift of approximately 5 × 10⁻⁵ eV per atom per ns. Long-range electrostatics were computed using the particle–particle particle–mesh Ewald method.³⁹ An orthorhombic ice I_h slab model made of 6144 H₂O molecules in a supercell with 3D periodic boundary conditions was used for simulations. The slab models were constructed from bulk ice I_h models, in which proton disorder is generated by a Monte Carlo algorithm that minimises the total dipole of the simulation cell.^40,41 Bulk ice was equilibrated in the constant pressure canonical ensemble at a set of temperatures from 190 to 240 K with a stride of 10 K. These well-equilibrated models were then cleaved by introducing a vacuum layer of 100 Å. This thickness of the vacuum layer is sufficient to prevent interactions between the two juxtaposed surfaces through periodic boundary conditions. The slabs were equilibrated for 50 ns in the canonical ensemble, enforced by stochastic velocity rescaling⁴² with a coupling constant of 0.1 ps. The same thermostat was used to perform production runs in the canonical ensemble. For each one of the temperatures considered, we ran production run simulations for 200 ns, and the structural properties reported in the next section (density profiles, surface structure function, and order parameter distributions) are calculated averaging over 2000 frames taken every 1 ps for each trajectory. Uncertainties are computed as the standard deviation over eight 25 ns-long blocks into which each trajectory is split. We have also verified that the differences among the blocks are due to statistical fluctuations rather than a systematic drift, thus proving that our production runs are sufficiently long. Furthermore, we averaged the structural properties over the two surfaces of the slab, after verifying that the computed properties are the same on each of them within the statistical error.

3 Results and discussion

To provide a microscopic insight into the surface structure of ice as a function of temperature and surface density of ion pairs, we computed density profiles, layer-resolved radial distribution functions, crystalline order parameter and surface roughness. Whereas we have conducted MD simulations with both low and high surface densities of ion pairs, hereafter we mostly focus on the comparison between pristine and briny ice surfaces with a surface density of 1 Na⁺/Cl⁻ pair per nm², i.e. a mass concentration of 47 parts per thousand (ppt) (see, for example, Fig. 1). The surface structure of ice with a low surface density of ion pairs (0.1 Na⁺/Cl⁻ pair per nm² or 5 ppt) do not exhibit significant differences from pure ice and the results for this system are reported in the ESI† (Fig. S1, S2, S7, S9 and S11).

Fig. 1 Snapshots of the molecular systems with the three varying surface surface densities of NaCl at the basal surface at 220 K.

Using the experimental cryoscopic constant of water, the freezing point depression of ice with a surface density of 1 NaCl pair nm⁻² of NaCl is 3 K. However, the TIP4P/2005 model slightly overestimates the value of the cryoscopic constant and the system melts at 245 K.³² We verified this prediction simulating the direct melting of ice slabs with sodium chloride ion pairs deposited at their surfaces.⁴³ Hence we consider temperatures up to 240 K, corresponding to −9 K below the melting temperature T_m of pure ice and −5 K below T_m with ions.

3.1 Density profiles

Fig. 2 shows the density profile in the direction perpendicular to the plane of the surface. The air–ice interface is on the left side of the plots. Both the basal and the primary prismatic surfaces consist of hydrogen-bonded bilayers. The density profiles show a double peak for each bilayer: a structural fingerprint, the absence of which provides a signature of disordering. As the temperature increases, the outermost bilayers smear and their density peaks decrease, while the valley separating them from the layers underneath increases. As observed in previous works, the outermost layer of pristine ice is already disordered at the lowest temperature considered (190 K). For example, for pristine ice the basal plane surface layer never exhibits the bilayer fingerprint (Fig. 2a). The double peak is preserved in the first prism surface up to 230 K, but the presence of a third low-intensity peak at air–ice interface suggests the partial formation of an adlayer (Fig. 2c). In the presence of ions the outermost bilayer of both the basal and the first prismatic surfaces are further smeared, and in the latter case, the bilayer structure and the adlayer-related peak are reduced to shoulders and eventually completely removed at 240 K (Fig. 2b and d). The difference in the melting behaviour of the second and third bilayer is even starker. The sudden increase of the valley between the first and the second layer of the basal plane of pristine ice confirms the “rounded” bilayer-by-bilayer melting transition identified in previous works.¹⁹ Conversely, in the presence of ions, the valley between the first and the second bilayer increases gradually, suggesting a continuous surface melting transition. The density profiles of the basal and the primary prismatic surfaces with ions (Fig. 2b and d) behave very similarly as a function of the temperature, suggesting that the ion pairs attenuate the structural differences between between these two surfaces.

Fig. 2 Oxygen atoms density profile as a function of temperature for the basal facet (top panels) without (a) and with ions (b) and primary prismatic facet (bottom panels) without (c) and with ions (d), averaged over 200 ns production runs. Statistical uncertainties are of the order of the line widths.

The Na⁺/Cl⁻ density profiles (Fig. 3) shows how ions penetrate in the subsurface layers as the temperature is increased from 220 to 240 K. At all temperatures Cl⁻ is slightly in excess in the outermost part of the QLL, in accordance with the surface propensity of halides observed in previous experiments and simulations of aqueous solutions.^44–49 We have verified that our modeling framework produces consistent results with previous studies on the surface propensity of anions in aqueous solutions, even if the Madrid model is non polarisable. The atomic density profiles computed for a slab of aqueous NaCl solution at 248 K exhibit a slight accumulation of Cl near the surface and a depletion of Na in agreement with former MD simulations with polarisable water models (Fig. S3, ESI†).^47,49 On ice, the accumulation of negative charge at the surface causes a peak of Na⁺ mass density in the first subsurface layer. Chloride is then prevalent in the layers underneath at the interface between QLL and crystalline ice. The density of chloride ions in the third bilayer increases gradually as the temperature rises from 220 to 230 K, and at 240 K Cl⁻ penetrates into the fourth bilayer. This trend is the same regardless of surface orientation. To shed light on the behaviour of ions in the QLL, we computed the charge density profiles and compared them with the particle densities of all four species in the system (O, H, Cl, Na) in Fig. S4 (ESI†). The peaks of the density of Cl⁻ align with the peaks of the density of H which have a positive partial charge, thus partially neutralizing the positive peaks of the charge density profile. Similarly, Na⁺ tends to compensate the negative charge in correspondence to the peaks of the oxygen density. The overall effect is a rounding of the charge density features, which is also the fingerprint of the onset of disorder.

Fig. 3 Cl⁻ and Na⁺ density profiles for the three highest temperatures sampled for both the basal and primary prismatic facet at the NaCl surface density of 1 NaCl pair nm⁻². Respective oxygen atom density profiles are included to show where the density of the ions is in relation to the bilayers. Ion density is scaled by a factor of 5×. The shaded regions corresponds to the standard deviation of the density profiles from block statistical analysis.

3.2 Structure of the quasi liquid layer

Looking further into the structural changes induced by ion pairs, Fig. 4 shows top-down snapshots of the second, third and fourth bilayers of the basal facet at 240 K. Although the QLL is dynamic and diffusive, these representative snapshots allow us to interpret the molecular mechanism of ice premelting. At the highest temperature (240 K), we see that for pristine ice (Fig. 4 top panels) the third and fourth bilayers retain the crystalline structure of hexagonal ice, even though the third bilayer exhibits topological defects.⁵⁰ The second bilayer is much more disordered but it still features areas with crystalline hexagonal structures. For the systems with 1 NaCl pair nm⁻² (Fig. 4 bottom panels), the second bilayer is completely disordered and contains a high surface density of ions – mostly Na⁺ – as also shown by the ion density profiles. The third bilayer is mostly disordered due to the presence of ions, mostly chloride, and, similarly to the second bilayer of pristine ice, it exhibits crystalline patches. The fourth bilayer is mostly crystalline, with defects similar to those observed in the third bilayer of pristine ice. However, it embeds chloride ions that replace water in their crystallographic positions. As chloride is easily embedded into the ice lattice, when the temperature increases Cl⁻ ions gradually penetrate into the deeper bilayers.^51,52 The negative ions cause a reorientation of the hydrogen bonding network around the crystallographic sites where they replace a water molecule. Additionally the preferred embedding of anions at the ice/QLL interfaces engenders a net charge that favours the formation of defects and destabilises the interfacial bilayer.⁵³ Previous work²¹ found that the first two bilayers of pristine ice are more weakly bound and more susceptible to premelting than those deeper in the ice slab. The snapshots in Fig. 4 suggest that the anions embedded in as far as the third and fourth bilayers destabilise the crystalline layer leading to a smoother premelting transition.

Fig. 4 Snapshots from a top down view of the second, third, and fourth bilayers of the basal facet with and without NaCl at 240 K. Bonds shown are for a bond length cutoff 1.9 Å to display the crystalline structure (or lack there of) in each layer. Na⁺ ions are represented by blue spheres and Cl⁻ ions are represented by green spheres.

Characteristic snapshots of both the basal and the primary prismatic surfaces at different temperatures are reported in Fig. S5 (ESI†) for the second bilayer of pristine ice and in Fig. S6 (ESI†) for third bilayer of ice with 1 NaCl pair nm⁻². The comparison between these two sets of figures suggests that there is a close correspondence in the behaviour of the second bilayer of pristine ice to the third bilayer of briny ice across the temperature range explored. In the latter system the inclusion of chloride marks the onset of the premelting transition, but extended disordering occurs only when sodium penetrates in the bilayer. Merging the information from the density profiles and Fig. 4 and Fig. S5 and S6 (ESI†) indicates that the presence of a high surface density of Na⁺/Cl⁻ ion pairs not only shifts the melting temperature but it also systematically increases the thickness of the QLL by a full bilayer. The structural disordering of subsurface layers corresponds to the incorporation of Na⁺ and Cl⁻ ions, with the latter penetrating in the deeper layers earlier than the former.

The temperature evolution of the layer-resolved RDFs of the surface and subsurface layers (Fig. S7, ESI†) contributes to support the hypothesis that the basal and primary prismatic surfaces, which have distinct premelting behaviour in freshwater ice, melt in a very similar way in the presence of ions and the premelting transition evolves more gradually as a function of temperature.

3.3 Local order parameter distribution

To shed light on the disordering of the premelted ice layer as a function of temperature we calculated the local q₆ Steinhardt order parameter^54–56 as implemented in former nucleation studies.^57,58Fig. 5 presents the density maps of the q₆ order parameter of the first four bilayers of the basal facet at different temperatures. The density maps of the primary prismatic facet and of the basal facet at low surface density of ions are presented in Fig. S8 and S9 in the ESI.† For reference, Fig. S10 (ESI†) reports the q₆ distribution for a water slab at 240 K. The latter is between 0.2 and 0.3 and does not depend on the vicinity to the surface. Fig. 5 shows that for pristine ice at 190 K there is a sharp distinction in the q₆ values among the layers, which at this temperature remain well separated. The surface layer is disordered but its q₆ distribution remains on the right of the liquid range. The subsurface layer presents a broad but mostly crystalline q₆ distribution and the bottom two layers can serve as reference for crystalline ice. With high surface density of Na⁺/Cl⁻ ion pairs the first and the second bilayer are not clearly distinct and the q₆ distribution spans the whole range from the liquid to the ice reference. In particular a substantial part of the second bilayer has q₆ values corresponding to the liquid phase. The third bilayer already shows the onset of disordering with a tail in the q₆ distribution extending below 0.5.

Fig. 5 Local Steinhardt's order parameter q₆ density map of the first four bilayers at each sampled temperature, for the basal facet without and with ions. Averaged over 200 ns. The scatter plots are colored by the calculated Gaussian kernel density estimation for the q₆ values, the scales values are arbitrary units but consistent for every density map.

As the temperature is increased to 210 K, the first two layers of pristine ice form a continuous distribution of q₆ with a marked dependence on the z-coordinate. The q₆ distribution does not change significantly up to 230 K, except for a the shift due to the thermal motion of the water molecules and the onset of dynamic exchange of molecules between the second and the third bilayer. In the system with Na⁺/Cl⁻ ion pairs at 200 K and 210 K the first and second bilayers form a single triangular q₆ distribution centred at ∼0.3 extending to 0.6 for the subsurface layer. At 220 K and above there is a single q₆ distribution that encompasses the top three bilayers and strongly depends on the z-coordinate. The q₆ distribution for the outermost part of the QLL is indistinguishable from liquid water, whereas ordering gradually increases with depth in the subsurface region.

At 240 K the q₆ of the third bilayer of pristine ice shifts to the left of the ice reference, indicating the onset of disordering. In the system with Na⁺/Cl⁻ ion pairs the q₆ distribution encompasses also the fourth bilayer, which retains values of the order parameters close to the ice reference. Similar trends can be observed by comparing the primary prismatic facet with and without 1 NaCl pair nm⁻² (Fig. S6, ESI†).

The q₆ maps confirm that the thickness of the QLL increases by about one layer when ion pairs are deposited at high surface density. Remarkably, ions modify the local structure of the QLL weakening the regular template induced by the crystalline subsurface layers, so that the surface has the same distribution of the local order parameter as liquid water. This observation agrees with the conclusions of experiments showing that the premelting layer of seawater ice has the same vibrational properties as liquid brine.²⁴ As the ice template is weaker, the differences in the structure and the premelting behaviour between the basal and the primary prismatic surfaces are also weakened for the deeper subsurface layers and vanish for the outermost layers.

3.4 Surface roughness

The roughness of ice surfaces affects the adsorption of reactive species and determines the crystal growth rates of facets with different orientations, thus controlling the shape of crystals formed by condensation from vapour.⁵⁹ The roughness of a surface is characterised by its structure function H(R)² defined as:

H(R)² = 〈(Δh(0) − Δh(R))²〉, (1)

where h(R) is the height of the surface at the planar distance R = |(R_x, R_y)| with respect to the plane origin, Δh(0) = h(0) − 〈h(R)〉, and 〈…〉 indicates an ensemble average that runs over the frames of the MD trajectory and the position of the origin. We define the surface by identifying the oxygen atoms, and sodium and chloride ions at the air/ice interface as those accessible to a spherical probe with a 3.5 Å radius.⁶⁰H(R) in eqn (1) can be rewritten in terms of the normalised height–height correlation function C(R) = 〈Δh(0)h(R)〉/〈Δh(0)²〉 and of the surface roughness as:

H(R)² = 2σ²[1 − C(R)]. (2)

A disordered surface is defined “flat” when it does not have long-range correlations. In this case, which is common for solid surfaces, C(R) decays exponentially with a characteristic length λ and the surface roughness squared σ² is the large-R limit of H(R)².⁶¹ The surface of a liquid, instead, is “rough” and H(R) diverges as log(R). The logarithmic divergence of the interfacial width at the air/liquid interface of water is due to the presence of capillary waves.⁶² Previous studies using the TIP4P/ice water model showed that, even in the presence of a premelting layer, the air–ice interface can be considered flat up to few degrees below the melting point, and it undergoes a sequence of phase transitions.^15,18 Here we revise these findings with the TIP4P/2005 water model and we address the effect of ions. To classify surface structures and transitions, we use the same nomenclature used in previous works addressing surface roughening and premelting.^18,63

Fig. 6 shows H(R)² at each temperature for the basal and primary prismatic facets with and without Na⁺Cl⁻ pairs. The curves show two distinct peaks for R < 7 Å which are characteristic of the short-range hydrogen bonding environment. At larger R the functions differ from one another with a marked dependence on ions surface density, orientation, and temperature. σ, λ, and H(R)² for the two different orientations of pristine ice are rather distinct. For the basal plane (Fig. 6a) from 190 to 230 K, H(R)² reaches its asymptotic roughness very rapidly and the value of σ² remains below 2 Ų, indicating that the surface is flat. At temperatures below 210 K the surface exhibits a very short correlation length <1.5 Å. H(R)² retains periodic fluctuations, which are the signature of both the non-diffusive character of the system and residual hexagonal patterning. Hence, we classify this phase as ordered flat (OF). The periodic features begin to fade for R ≳ 10 Å starting at 210 K and λ increases to ∼3 Å. This indicates a surface phase transition from OF to disordered flat (DOF), which is accompanied by the onset of surface diffusion. At 240 K H(R)² retains similar characteristics as the lower temperatures but the roughness jumps beyond 2 Å, and λ reaches a value of approximately 4.3 Å, indicating a second structural transition to a high temperature reconstructed flat phase (HTRF).

Fig. 6 Structure function as a function of temperature of the basal facet of pristine ice (a) and of ice with sodium chloride ions (b), and of the primary prismatic facet of pristine ice (c) and of ice with sodium chloride ions (d). The shaded regions corresponds to the standard deviation of the structure functions from block statistical analysis.

Conversely, the H(R)² of the primary prismatic surface (Fig. 6c) displays large periodic fluctuations from 190 to 220 K. The convergence toward the asymptotic values is slower, and in two cases (210 and 220 K) the correlation length is comparable or larger than the maximum radius of a sphere that can be enclosed in the simulation cell. The long-range fluctuations indicate that this surface retains the fingerprint of crystalline ordering, although the structure is affected by structural defects that lead to the formation of an adlayer of molecules (see also Fig. 2). Hence, we dubbed this surface ordered flat (OR). The roughness parameter σ² is non-monotonic with the temperature. It reaches a maximum at 210 K, drops significantly at 230 K and increases again at 240 K. At 230 K the periodic features of H(R)² vanish indicating a structural transition to HTRF with λ values of approximately 4.3 Å and roughness similar to that of the basal plane at 240 K.

The sequence of preroughening transitions of pristine ice surfaces is similar to those observed in ref. 18 and correlate with changes in the growth behaviour in the Nakaya diagram, which classifies the shape of ice crystals growing from vapour deposition at different temperatures and humidity (Fig. 7).^59,64,65 The transition from OF to DOF of the basal surface, between T_m −49 and −39 K, corresponds to the change from columnar ice crystals to plates at low vapour pressure. The transition to HTRF for the primary prismatic surface at about T_m −25 K corresponds to the preference for plates at any vapour pressure. Finally, the transition to HTRF of the basal surface at T_m −15 K correlates with the occurrence of columnar crystals.

Fig. 7 Schematic of ice crystal growth at low supersaturation in relation to the correlation length values obtained from Fig. 6. Ice crystal shape temperature dependent transitions obtained from ref. 6.

H(R)² for the systems with 1 NaCl pair nm⁻² (Fig. 6b and d) indicates that these surfaces are much rougher than the corresponding pristine ice facets. Additionally, even at the lower temperatures, no features arising from a periodic pattern are present, thus suggesting that solvated Na⁺/Cl⁻ ion pairs disrupt long-range order at the air–ice interface. The primary prismatic surface is generally rougher than the basal surface up to 230 K. At 190 K, the basal surface with ions has a similar structure function and correlation length to that of pristine ice in the DOF phase, which occurs at higher temperature. At temperatures between 200 and 230 K, the basal and primary prismatic surfaces with ions have the same features and correlation lengths, and can be classified as HTRF. Above 230 K both surfaces undergo a roughening transition that leads to a logarithmically divergent H(R). The structure function of this phase, named premelted rough, is indistinguishable from that of a liquid water slab at a similar temperature (Fig. 8). These results suggests that the surface structure of ice with a high density of NaCl is characterized by the occurrence capillary waves. Whereas previous works suggested that anions pin equilibrium capillary waves at the surface of either water or ice,^23,48 we do not observe this effect in our simulations. Fig. 8 shows that, in spite of the surface propensity of Cl⁻, neither the surface of a liquid slab is affected by the dissolved NaCl salt, nor we observe the pinning of capillary waves in the premelting layer.

Fig. 8 Structure function of the surface of water and of aqueous NaCl solution with a concentration of 45 ppt at 240 K, computed for a slab of 6144 water molecules. The two structure functions overlap, and exhibit logarithmic divergence.

Above 200 K (T − T_m ∼ −40 K) the surface roughening transitions occur at the same temperatures for both facets, thus affecting the shape of growing ice crystals. We indeed expect that the presence of ion pairs at typical atmospheric temperatures would level the growth rates of the two low-index surfaces and lead to the formation isotropic crystals.

4 Conclusions

We have characterised the structure and melting behaviour of both the basal and primary prismatic surfaces of ice I_h in the presence of sodium chloride at both high and low surface density. We find that the addition of ions in low surface density does not change the structural properties of the QLL with respect to freshwater ice. Conversely, high surface density of ion pairs corresponding to seawater ice conditions leads to a significant change in QLL thickness, structure and premelting behaviour. Na⁺/Cl⁻ ion pairs disrupt the residual ordering of QLL and increase its thickness by at least one bilayer at all the temperatures considered. Surface premelting becomes more gradual in the presence of ions and the differences in ordering and premelting behaviour between the two surfaces are erased. Close to the melting temperature (T − T_m ∼ −10 K) the QLL with ions exhibits the same structural features as a water slab, thus providing the same solvation environment as a liquid surface.²⁴

The changes induced by Na⁺/Cl⁻ ion pairs at the molecular level lead to major chances in the surface structural transitions as a function of temperature, erasing the differences between the two facets. This effect may impact the growth rate and the shape of ice crystals growing from the vapour phase in clouds containing significant amounts of sodium chloride aerosols. The observed changes in surface roughness upon ion adsorption may also impact the light scattering properties of ice and snow. The observation that adsorbed ions increase surface roughness, knowing that rougher surfaces decrease albedo, would explain why sea ice has a lower albedo than freshwater snow.

In general terms, the key implication of our study is the crucial need to consider the surface density of Na⁺/Cl⁻ present when modelling ice surfaces, given its ability at high surface density to alter the surface structure, thus affecting mechanistic properties such as glacier sliding,⁶⁶ and chemistry at the ice surface.^67–69

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We are grateful to Luis G. MacDowell for useful discussions and the critical reading of the manuscript. This work is supported by the National Science Foundation under Grant No. 1806210 and Grant No. 2053235.

References

1. M. A. Tolbert, M. J. Rossi, R. Malhotra and D. M. Golden, Science, 1987, 238, 1258–1260.
2. L. Chu and C. Anastasio, J. Phys. Chem. A, 2003, 107, 9594–9602.
3. A. R. Ravishankara, Science, 1997, 276, 1058–1065.
4. A. S. McFall, K. C. Edwards and C. Anastasio, Environ. Sci. Technol., 2018, 52, 5710–5717.
5. T. Hullar, F. C. Bononi, Z. Chen, D. Magadia, O. Palmer, T. Tran, D. Rocca, O. Andreussi, D. Donadio and C. Anastasio, Environ. Sci.: Processes Impacts, 2020, 22, 1666–1677.
6. B. Slater and A. Michaelides, Nat. Rev. Chem., 2019, 3, 172–188.
7. M. J. Shultz, Annu. Rev. Phys. Chem., 2017, 68, 285–304.
8. Y. Nagata, T. Hama, E. H. G. Backus, M. Mezger, D. Bonn, M. Bonn and G. Sazaki, Acc. Chem. Res., 2019, 52, 1006–1015.
9. M. A. Sánchez, T. Kling, T. Ishiyama, M.-J. van Zadel, P. J. Bisson, M. Mezger, M. N. Jochum, J. D. Cyran, W. J. Smit, H. J. Bakker, M. J. Shultz, A. Morita, D. Donadio, Y. Nagata, M. Bonn and E. H. G. Backus, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 227–232.
10. H. Bluhm, D. F. Ogletree, C. S. Fadley, Z. Hussain and M. Salmeron, J. Phys.: Condens. Matter, 2002, 14, L227–L233.
11. H. Dosch, A. Lied and J. Bilgram, Surf. Sci., 1995, 327, 145–164.
12. G. Sazaki, S. Zepeda, S. Nakatsubo, M. Yokomine and Y. Furukawa, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 1052–1055.
13. K.-i Murata, H. Asakawa, K. Nagashima, Y. Furukawa and G. Sazaki, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, E6741–E6748.
14. M. M. Conde, C. Vega and A. Patrykiejew, J. Chem. Phys., 2008, 129, 014702.
15. J. Benet, P. Llombart, E. Sanz and L. G. MacDowell, Phys. Rev. Lett., 2016, 117, 096101.
16. T. Kling, F. Kling and D. Donadio, J. Phys. Chem. C, 2018, 122, 24780–24787.
17. Y. Qiu and V. Molinero, J. Phys. Chem. Lett., 2018, 9, 5179–5182.
18. P. Llombart, E. G. Noya and L. G. MacDowell, Sci. Adv., 2020, 6, eaay9322.
19. P. Llombart, E. G. Noya, D. N. Sibley, A. J. Archer and L. G. MacDowell, Phys. Rev. Lett., 2020, 124, 065702.
20. N. Kastelowitz, J. C. Johnston and V. Molinero, J. Chem. Phys., 2010, 132, 124511.
21. M. Watkins, D. Pan, E. G. Wang, A. Michaelides, J. VandeVondele and B. Slater, Nat. Mater., 2011, 10, 794–798.
22. X. Wei, P. B. Miranda, C. Zhang and Y. R. Shen, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 66, 085401.
23. S. P. Niblett and D. T. Limmer, J. Phys. Chem. B, 2021, 125, 2174–2181.
24. T. F. Kahan, S. N. Wren and D. J. Donaldson, Acc. Chem. Res., 2014, 47, 1587–1594.
25. S. Zimmermann, M. Kippenberger, G. Schuster and J. N. Crowley, Phys. Chem. Chem. Phys., 2016, 18, 13799–13810.
26. T. F. Kahan and D. J. Donaldson, Environ. Sci. Technol., 2010, 44, 3819–3824.
27. L. Krnavek, W. Simpson, D. Carlson, F. Domine, T. Douglas and M. Sturm, Atmos. Environ., 2012, 50, 349–359.
28. J. L. F. Abascal and C. Vega, J. Chem. Phys., 2005, 123, 234505.
29. E. Sanz, C. Vega, J. L. F. Abascal and L. G. MacDowell, Phys. Rev. Lett., 2004, 92, 255701.
30. J. L. F. Abascal and C. Vega, Phys. Rev. Lett., 2007, 98, 237801.
31. C. Vega and J. L. F. Abascal, Phys. Chem. Chem. Phys., 2011, 13, 19663.
32. M. M. Conde, M. Rovere and P. Gallo, J. Chem. Phys., 2017, 147, 244506.
33. R. García Fernández, J. L. F. Abascal and C. Vega, J. Chem. Phys., 2006, 124, 144506.
34. A. L. Benavides, M. A. Portillo, V. C. Chamorro, J. R. Espinosa, J. L. F. Abascal and C. Vega, J. Chem. Phys., 2017, 147, 104501.
35. M. M. Conde, M. Rovere and P. Gallo, J. Mol. Liq., 2018, 261, 513–519.
36. I. de Almeida Ribeiro, R. Gomes de Aguiar Veiga and M. de Koning, J. Phys. Chem. C, 2021, 125, 18526–18535.
37. H. J. C. Berendsen, D. van der Spoel and R. van Drunen, Comput. Phys. Commun., 1995, 91, 43–56.
38. D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. C. Berendsen, J. Comput. Chem., 2005, 26, 1701–1718.
39. J. W. Eastwood, R. W. Hockney and D. N. Lawrence, Comput. Phys. Commun., 1984, 35, C-618.
40. N. Grishina and V. Buch, J. Chem. Phys., 2004, 120, 5217–5225.
41. M. Matsumoto, T. Yagasaki and H. Tanaka, J. Comput. Chem., 2018, 39, 61–64.
42. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
43. C. Vega, M. Martin-Conde and A. Patrykiejew, Mol. Phys., 2006, 104, 3583–3592.
44. L. Piatkowski, Z. Zhang, E. H. G. Backus, H. J. Bakker and M. Bonn, Nat. Commun., 2014, 5, 4083.
45. E. M. Knipping, M. J. Lakin, K. L. Foster, P. Jungwirth, D. J. Tobias, R. B. Gerber, D. Dabdub and B. J. Finlayson-Pitts, Science, 2000, 288, 301–306.
46. S. Ghosal, J. C. Hemminger, H. Bluhm, B. S. Mun, E. L. D. Hebenstreit, G. Ketteler, D. F. Ogletree, F. G. Requejo and M. Salmeron, Science, 2005, 307, 563–566.
47. P. Jungwirth and D. J. Tobias, Chem. Rev., 2006, 106, 1259–1281.
48. D. E. Otten, P. R. Shaffer, P. L. Geissler and R. J. Saykally, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 701–705.
49. D. J. Tobias, A. C. Stern, M. D. Baer, Y. Levin and C. J. Mundy, Annu. Rev. Phys. Chem., 2013, 64, 339–359.
50. D. Donadio, P. Raiteri and M. Parrinello, J. Phys. Chem. B, 2005, 109, 5421–5424.
51. E. J. Smith, T. Bryk and A. D. J. Haymet, J. Chem. Phys., 2005, 123, 034706.
52. L. Vrbka and P. Jungwirth, Phys. Rev. Lett., 2005, 95, 148501.
53. J. S. Kim and A. Yethiraj, J. Chem. Phys., 2008, 129, 124504.
54. P. J. Steinhardt, D. R. Nelson and M. Ronchetti, Phys. Rev. B: Condens. Matter Mater. Phys., 1983, 28, 784–805.
55. P. R. ten Wolde, M. J. Ruiz-Montero and D. Frenkel, J. Chem. Phys., 1996, 104, 9932–9947.
56. A. Reinhardt, J. P. K. Doye, E. G. Noya and C. Vega, J. Chem. Phys., 2012, 137, 194504.
57. D. Moroni, P. R. ten Wolde and P. G. Bolhuis, Phys. Rev. Lett., 2005, 94, 235703.
58. G. C. Sosso, J. Chen, S. J. Cox, M. Fitzner, P. Pedevilla, A. Zen and A. Michaelides, Chem. Rev., 2016, 116, 7078–7116.
59. K. G. Libbrecht, Annu. Rev. Mater. Res., 2017, 47, 271–295.
60. M. Sega, G. Hantal, B. Fábián and P. Jedlovszky, J. Comput. Chem., 2018, 39, 2118–2125.
61. K. Rommelse and M. den Nijs, Phys. Rev. Lett., 1987, 59, 2578–2581.
62. A. E. Ismail, G. S. Grest and M. J. Stevens, J. Chem. Phys., 2006, 125, 014702.
63. E. A. Jagla, S. Prestipino and E. Tosatti, Phys. Rev. Lett., 1999, 83, 2753–2756.
64. U. Nakaya, Snow Crystals, Cambridge, MA, Harvard University Press, 1954.
65. K. G. Libbrecht, Rep. Prog. Phys., 2005, 68, 855–895.
66. K. M. Cuffey, H. Conway, A. Gades, B. Hallet, C. F. Raymond and S. Whitlow, J. Geophys. Res., 2000, 105, 27895–27915.
67. T. D. Shepherd, M. A. Koc and V. Molinero, J. Phys. Chem. C, 2012, 116, 12172–12180.
68. T. Bartels-Rausch, H.-W. Jacobi, T. F. Kahan, J. L. Thomas, E. S. Thomson, J. P. D. Abbatt, M. Ammann, J. R. Blackford, H. Bluhm, C. Boxe, F. Domine, M. M. Frey, I. Gladich, M. I. Guzmán, D. Heger, T. Huthwelker, P. Klán, W. F. Kuhs, M. H. Kuo, S. Maus, S. G. Moussa, V. F. McNeill, J. T. Newberg, J. B. C. Pettersson, M. Roeselová and J. R. Sodeau, Atmos. Chem. Phys., 2014, 14, 1587–1633.
69. J. P. D. Abbatt, Chem. Rev., 2003, 103, 4783–4800.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2cp02277j

---