Article pubs.acs.org/JPCA

Ionization of Acids on the Quasi-Liquid Layer of Ice S. Riikonen,*,† P. Parkkinen,† L. Halonen,† and R. B. Gerber†,‡,§ †

Laboratory of Physical Chemistry, Department of Chemistry, University of Helsinki, P.O. Box 55, FI-00014, Helsinki, Finland Institute of Chemistry and the Fritz Haber Research Center, The Hebrew University, Jerusalem 91904 Israel § Department of Chemistry, University of California Irvine, Irvine, California 92697, United States ‡

S Supporting Information *

ABSTRACT: The ice quasi-liquid layer (QLL) forms on ice surfaces below the bulk ice melting temperature. It is abundant in the atmosphere, and its importance for atmospheric chemistry is recognized. In the present work, we have studied the microscopic mechanisms of acid ionization on the QLL using ab initio molecular dynamics. The model system QLL is established by nanosecond time scale simulations with empirical force fields, while the reactivity of the QLL is studied using ab initio molecular dynamics. Our ab initio simulations reveal that QLL is reactive, exhibiting stable crystalline point defects, which contribute to efficient acid solvation, ionization, and proton transfer. We study in detail deuterated hydrogen iodide (DI) and nitric acid (DNO3). Ionization in both cases benefits from the abundance of weakly bonded hydrogen-bond single-acceptor double-donor water molecular species available on the QLL in high relative concentration. Picosecond time scale ionization is demonstrated for both molecular species. Our results suggest efficient reactivity of acid ionization and proton transfer at temperature ranges appropriate for the upper troposphere and lower stratosphere.

■

INTRODUCTION Premelting of ice surfaces is a phenomenom emerging at temperatures below the bulk melting temperature. This region between vacuum and bulk ice is referred to as the quasi-liquid layer (QLL). Its existence has been a subject of scientific debate since Faraday’s papers on sintering of ice.1 As the QLL exists on ice surfaces under atmospheric conditions, its importance for various physical and chemical environmental processes has to be appreciated.2 The premelting of an ice surface can be divided into three temperature regions:3 (i) At approximately −183 °C, a perfect, crystalline ice surface is preserved, accompanied by large amplitude vibrational motions of the surface molecules, as suggested by low-energy electron diffraction (LEED) studies.4,5 (ii) Approaching −30 °C, these motions become larger, creating crystalline point defects, from which water molecules are expelled to a vacuum. The hydrogen-bond network becomes increasingly distorted as the temperature rises. (iii) After a certain temperature, the topmost layers can be considered to be liquid water, while the thickness of the QLL diverges6 and penetrates deeper into the bulk. Such behavior a “quasi-liquid” layer near the crystalline bulk region and a truly liquid layer on top of the surface having a diverging thickness with increasing temperatureis a general phenomenom of bulk melting processes (see refs 2−7 and references therein). The temperatures cited above should be considered only approximately, as these can depend on the experimental conditions and the orientation of the ice crystal surface. The “onset” of the QLL is usually placed somewhere between i and ii, when molecular orientations and motions become distinguishable from those of ice. References using techniques © XXXX American Chemical Society

with varying surface sensitivities have reported several onset values ranging from −73 °C to a few degrees below the bulk melting temperature8−12 (see particularly ref 8 and references therein). The recent value of ≈ −73 °C refers to increased orientational disorder of the ice surface dangling OH bonds, and has been obtained using the surface sensitive technique of sum frequency generation for the ice/vapor interface.11,12 The most intriguing and complex behavior occurs in the truly “quasi-liquid“ (i.e., neither ice nor liquid) region between ii and iii, which consists of crystalline point defects, broken hydrogen bonds, and large amplitude motion of the molecules. Estimates of the thickness of the QLL as a function of temperature vary with different experimental methods and conditions (see ref 8 and references therein). Some more recent studies6,13−15 give values of the same order of magnitude for the basal plane. Namely, photoelectron studies suggest that the QLL can be thin (∼10 Å, 1 Å = 10−10 m), in accordance with theoretical simulations16 also suggesting that the QLL is far from being completely “liquid-like“ even at ≈ −5 °C.16 In the present work, we study the chemical reactivity and ionization of atmospherically important acidic species in region ii. Such low-temperature ice surfaces are abundant in the atmosphere: the upper troposphere and lower stratosphere contain cirrus clouds, consisting of water ice particles and covering ≈20% of the globe.17−19 Polar stratospheric clouds have long been recognized as scaffolds for heterogeneous chemistry leading to ozone depletion.20,21 Snow can cover up Received: June 6, 2014 Revised: June 13, 2014

A

dx.doi.org/10.1021/jp505627n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

suggest that point defects and ledge sites are reactive for acid ionization: weak water hydrogen bonds, available at the defect sites, can also be exploited by the acid molecule in its effort to solvate and donate the proton.39 Ledge sites and point defects can be abundant in cirrus cloud ice particles, as they exhibit increased surface roughness.45 In the present work, we study acid ionization on a realistic QLL modelto our best knowledgefor the first time. The surface molecular structure we employ is consistent with the extensive empirical potential studies.16,40−44 Details of the lowtemperature QLL model system and simulation techniques are given in the Methods section and in the Supporting Information. In the following, we start by discussing the model system and its main features. All model systems considered are deuterated. When referring to concepts such as “proton transfer”, “proton donor”, “hydrogen bond”, etc., deuterium is assumed. Empirical force field calculations using the SPC/E potential46 were used to equilibrate the systems in the nanosecond time scale, while subsequent reaction ab initio studies were done with density functional theory (DFT). Both simulation stages were done with the CP2K/Quickstep program package.47 More details about the computational scheme are given in the Supporting Information.

to 50% of the landmass of the northern hemisphere and thus creates a massive and porous media where atmospheric trace gases interact with ice surfaces.22 Important examples are the high concentrations of HONO and other nitrogen contaning species emitted by arctic snow.22−24 The importance of acidic trace gases interacting with ice surfaces25 has inspired many laboratory studies on the ionization of acids in low-temperature ice nanocrystals26,27 and films.28−30 Probably the most studied acid in this context is the hydrochloric acid (HCl). While it is known to ionize on liquid water and on water ice surfaces, the extent of its ionization at cryogenic temperatures (T ≈ −220 °C and below) is being debated:26−28,31−34 Fourier transform infrared (FTIR) spectra of small nanoparticles suggest only ∼15% ionization at 30% coverage, and almost solely due to HCl self-solvation.26 Classical force-field molecular dynamics simulations by the same authors suggest that the nanoparticle ice surface (model particles with ∼1 nm diameter) is “passive” as HCl molecules only rarely find sites that promote ionization.27 Reactive ion scattering studies of Kang et. al.31 on vapor-grown amorphous surfaces give a similar conclusion: no ionization of HCl at T ≈ −220 °C is observed. On the other hand, Parent et. al., based on their RAIRS and NEXAFS spectra, suggest 92% ionization of HCl at 0.16 ML coverage, on crystalline and amorphous ice surfaces at the same temperature. Similar results were obtained by Ayotte et al.28 using RAIRS: HCl ionized at T ≈ −250 °C starting from 0.2 ML HCl concentration on crystalline ice, while similar results, i.e., extensive ionization, have been obtained for nitric acid (HNO3) at T ≈ −230 °C. A possible explanation for the conflicting results could be the different ice morphologies present in the various experiments, as speculated by Ayotte.28 Devlin and Kang33 have suggested that the findings of Parent et. al.32,34 are due to HCl self-solvation and/ or mismeasuring the HCl dose, also calling for a reinterpretation of some of the features in their infrared spectra.33,34 The hydrogen-bonding topologies can be different for amorphous ice, crystalline ice, and rounded nanoparticle surfaces.35 Amorphous, vapor-deposited ice also comes in two different forms (sparse and dense),36 while both of these may, in principle, feature unique surface morphologies. Suggesting an efficient reactivity of crystalline ice, theoretical calculations of Bolton et al.37,38 show that a perfectly crystalline ice basal plane does offer a large number of sites for potential HCl ionization. Stable crystalline point defects and ledge sites result in the increase of such sites.38,39 This situation can be contrasted to the nanoparticle case, where the surface reconstructs considerably, eliminating dangling OH groups.27 The QLL lies in the temperature region between lowtemperature experiments and liquid water. On one hand, one might expect the elimination of dangling OH groups due to the “amorphization” of the strictly crystalline ice surface, while, on the other hand, crystalline point defects might arise, creating similar reactive sites as seen in earlier ab initio studies. The surface topology and reactivity of the QLL may vary in a similar manner as in the corresponding cryogenic systems (i.e., amorphous/crystalline slabs or small nanoparticles at ≈ −220 C), while, in the present case, we are interested in the QLL layer forming on the crystalline ice basal plane. A long-standing effort of empirical water potential studies has been to characterize the QLL on crystalline ice facets at the molecular level.16,40−44 In some such studies, the existence of long-lived ice surface crystalline point defects has been emphasized.40,42,44 As discussed above, theoretical calculations

■

METHODS The starting point for the QLL model system was a slab of hexagonal ice consisting of four bilayers, having lateral dimensions of 17.5 Å × 15.2 Å and a total of 128 water molecules. During the MD simulations, the lowermost bilayer was held fixed. Nanosecond time scale simulations and equilibration were performed using an empirical, nonpolarizable (ENP) water potential, the SPC/E model46 and ab initio density functional theory (DFT) MD simulations were performed on the geometries obtained from ENP calculations. The aim of this scheme is to establish a valid model system complying with some key properties of a low-temperature QLL. A large body of theoretical work using ENPs on the quasi-liquid layer emphasizes the need of nanosecond time scale simulations to achieve convergent QLL properties.16,40−44 Different atomic potentials (ENP vs ab initio) and simulation conditions (unit cell sizes, fixed molecules) exhibit different bulk ice melting temperatures. This makes the estimation of undercooling temperatures in the present scheme difficult, even impossible. In the Supporting Information, this problem is addressed. There, we confirm that our simulations are performed in the QLL region by comparing the self-diffusion coefficient and relative mean square displacement of molecules. The DFT calculations were performed using the BLYP48 exchange and correlation functional, GTH pseudopotentials, and the Grimme D3 dispersion corrections,49 as implemented in the CP2K/Quickstep program package.47 The classical forcefield SPC/E simulations were performed using the CP2K/Fist module. Nanosecond time scales were achieved using SPC/E, while ab initio MD runs consisted of tens of picoseconds. Target temperatures were obtained using the canonical ensemble with massive thermostatting.50 Production runs were microcanonical. In reaction studies, acid molecules were placed in close contact with the equilibrated QLL surface and the system was left to evolve within the microcanonical ensemble. The DZVP-molopt-SR basis51 set was utilized in the DFT calculations. B

dx.doi.org/10.1021/jp505627n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

Figure 1. Normalized histogram of oxygen positions in the direction of the normal to the slab for the QLL model system, from a time period of 8 ps at a temperature of T = 258 K. The average number of water molecules in layers is indicated in the top part of the figure. These are obtained from the integration of the peak areas. Layers are labeled with roman numbers. 32 oxygen atoms correspond to a full bilayer. The distance Z = 0 has been chosen in the middle of the slab.

Figure 2. QLL model slab viewed from the top. Instaneous atomic configuration and hydrogen bonding of layers v and vi are illustrated in the background. A rectangle shows the dimensions of the unit cell, while arrows indicate layer vi molecules. (a) Histograms of layer v and vi oxygen positions in green and blue, respectively, from the time interval and temperature of Figure 1. A blue dashed line indicates a migration event of a molecule in layer vi. (b) Histogram of free oxygen lone-pair positions in gray color from a time interval half of that of panel a. Dashed polygons (panel a) trace out some hydrogen-bonding motifs, resulting from the crystalline point defects.

■

RESULTS AND DISCUSSION Characteristics of the QLL. Figure 1 shows an oxygen atom z-profile histogram from an ab initio simulation (equilibrated for ≈2 ns with an empirical SPC/E force field, then ≈19 ps with ab initio, and, finally, collecting data for 10 ps). It features well-ordered deeper bilayers, while the topmost (bi)layers v and vi exhibit large amplitude motions. We emphasize that the shape of the profile is similar to empirical potential studies using much larger model slabs16,43,44 where the topmost “quasi-liquid-like” layers v and vi are referred to as ϵ2 and ϵ1. The number of molecules in each (bi)layer is determined by integrating the area of the oxygen atom profile peaks, 32 implying an intact bilayer in our model. According to Figure 1, 10% of v and vi molecules reside in the topmost layer (vi), in an agreement with the 8% value reported in ref 43 and 7.5% reported in ref 44. Figure 1 also suggests an interpretation for the QLL model system: the topmost layers v and vi consist of crystalline point defects and admolecules, respectively: layer v, which still maintains a double peak, is missing approximately three molecules (point defects).

These molecules that were expelled from layer v form layer vi that lacks any ice-like structure. Kroes40 came to a similar interpretation. The presence of “gaps” (i.e., point defects) was also emphasized in ref 44. In the Supporting Information, we have studied the self-diffusion coefficient of water molecules in layers v and vi. The value we obtain using SPC/E simulations is less than 10−6 cm2/s, which is consistent with the experimental results of the QLL in the temperature range −23 → −15 °C.52,53 Considering atmospheric trace gases impinging upon the QLL, empirical potential studies suggest that ice-like bilayers block the penetration of adsorbants deeper into the bulk,54,55 emphasizing further the importance of layers v and vi for atmospheric reactions. To further establish the interpretation of a relatively static ice surface consisting of crystalline point defects and admolecules, we have plotted the xy profile histogram of layers v and vi in Figure 2a. In the 8 ps time scale of the histogram, molecules maintain their positions, while one significant migration event in layer vi takes place. A visual inspection of the instantaneous C

dx.doi.org/10.1021/jp505627n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

Figure 3. Left panel: Probability that a water molecule found along the surface normal is of certain proton-acceptor/donor species. The QLL from the 10 ps time period ab initio simulation. The symbol “A” stands for the hydrogen-bond (proton) acceptor, while “D”, for the donor. The existence of a hydrogen bond is judged from a geometric criterion, namely, from the O−O distance ≤ 3.5 Å and the O−H···O angle ≤ 140°.61 Right panel: Schematic illustration for the acceptor/donor species profile in a perfect crystalline ice at T = 0 K. Ice oxygen atom positions are plotted as black spheres, while hydrogen atoms are omitted. The ideal basal plane has an equal number of ADD and AAD molecules in the topmost layer (eight in our model system), while all other molecules are of the AADD type. For both panels, the number of ADD and AAD molecules in the layers are indicated by N[ADD] and N[AAD], respectively.

Figure 4. Upper part: Simulated infrared spectra of the deuterated crystalline ice (ICE, T = 96 K), the quasi-liquid layer (QLL, T = 258 K), and the vapor−liquid interface (VLI, T = 317 K). The lower part (negative values): Simulated infrared spectra of crystalline ice, decomposed into contributions from two molecular types, namely, from AAD (hydrogen-bond double-acceptor single-donor, i.e., “dangling hydrogen”) and ADD (hydrogen-bond single-acceptor double-donor, i.e., “dangling oxygen”) species.

fluctuate, implying weak hydrogen bonds. Figure 2b shows a histogram of free lone pairs from a 4 ps time period for layers v and vi. A technique for distinguishing free oxygen lone pairs, based on Wannier function shifts, is given in the Supporting Information. The most reactive species are the three “admolecules” of layer vi. The broken hydrogen-bond network of layer v exposes a significant density of free lone pairs and, thus, weak hydrogen bonds. It is worthwile to study the concentration of different molecular species on this “frozen” and broken hydrogenbonded system. Enumerating the hydrogen bonds of molecules results in the “A/D” notation; i.e., the dangling oxygen species is “ADD” (hydrogen-bond proton single-acceptor, doubledonor), while dangling hydrogen “AAD” (hydrogen-bond proton double-acceptor, single-donor). The number of ADD−ADD and AAD−AAD hydrogen bonds is known to mandate the energetics of 3-fold coordinated water molecule systems,56−58 while the relative concentration of these molecular types has been debated in the case of the vapor− liquid interface.59−61

water molecule positions and hydrogen bonding of Figure 2 reveals a broken hydrogen-bonded network, consisting of squares, pentagons, octagons, etc., instead of the hexagons of the perfect basal plane. This is a consequence of the crystalline point defects at layer v. Layer vi molecules sit on top of this structure, which, during the course of tens of picoseconds, slowly transforms to another hydrogen-bonded motif. The persistence of the molecular structure described above is also evident from the experimental and theoretical calculations of the QLL self-diffusion coefficient, which is orders of magnitude smaller than the self-diffusion coefficient of liquid water (see the Supporting Information). Then, summarizing, the QLL features crystalline point defects in layer v, resulting in a broken, yet stable hydrogen-bonded network. We can obtain a qualitative idea of the reactivity of this broken and stabilized hydrogen-bonding structure by studying how many oxygen atom lone pair electrons it effectively exposes: such “free” lone pair electrons are not involved in a hydrogen bond, and they can accept the incoming proton. They can also be available instantaneously, as the hydrogen bonds D

dx.doi.org/10.1021/jp505627n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

Figure 5. Left panel: Forming the NO3−D+ contact ion pair (CIP) on the QLL: dashed arrows denote a hydrogen bond (proton) donation. The symbol “ADD” stands for a single-acceptor double-donor water molecule species, “DD” for a double-donor species, etc. Right panel: Grothuss migration of the hydronium species. The thickness of the green lines indicates the relative amount of time for which the proton Eigen/Zundel cation stays localized between certain oxygen atoms during 10 ps of an ab initio MD simulation.

contrasted to the more isotropic solvation of a monatomic base, say, to iodide I−. In a recent study, we observed that, when placed and positioned correctly into a crystalline point defect in the ice-Ih basal plane, nitric acid can ionize effectively, at picosecond time scales.39 In the present case, we will study if the broken, yet stable hydrogen-bonded network of the QLL offers possibilities for rapid ionization. As discussed above, the QLL structure features spontaneously formed crystalline point defects, suggesting enhanced reactivity. We placed DNO3 on top of the deuterated QLL model surface at seven different sites of the unit cell depicted in Figure 2a, while the positions and momenta of water molecules were taken from stabilized (NVT, canonical ensemble) and equilibrated (NVE, microcanonical ensemble) ab initio runs (see Table S3 of the Supporting Information for more details). The composite DNO3/QLL system then evolved during ≈20 ps in an ab initio NVE run. More information on the seven simulation runs is given in the Supporting Information (in the following, we will discuss in detail the trajectory labeled as “A” in the Supporting Information). The mechanisms in all trajectories can be summarized by Figure 5 and in a few discrete events (1−4): (1) Once the nitric acid proton creates a hydrogen bond to a water moleculecall this the “target” moleculethat hydrogen bond is persistent throughout the simulation. (2) The type of the target molecule mandates if ionization will happen: an acid−ADD complex is a necessary condition for ionization (see Figure 5). (3) Prior to ionization, nitric acid must be solvated with N ≥ 2 hydrogen bonds. (4) Pathways to achieve an acid−ADD complex are illustrated in the left side of Figure 5: fluctuations of the hydrogen-bond network can lead to a DD species, or alternatively, nitric acid can “steal” a hydrogen bond. Even more complex co-operative behavior in forming the acid−ADD system was observed in some simulation runs, along the lines of Figure 3e,f of ref 39. Of the seven MD runs, two resulted in ionization. In these cases, both (2) an acid−ADD complex formed and (3) nitric

Figure 3 shows the probability of various molecular types along the normal of the QLL surface. The number of ADD and AAD molecular species on the ideal ice basal plane is compared to their numbers in the QLL. In layers v and vi of the QLL, the number of ADD and AAD species is approximately 40% less than in the ideal basal plane. This implies the elimination of dangling OH groups (AAD species) and, at first sight, seems to create an unfavorable solvating environment. However, while on the perfect basal plane the ADD and AAD species are tetrahedrally coordinated and strongly bonded, on the QLL, they form part of a broken hydrogen-bond network in layer v, are weakly bonded, and can take part in solvation and proton transfer, which involves hydrogen-bond breaking (see, for example, Figure 5). From Figure 3, we see that, at the boundary of layers v and vi of the QLL, the most dominant, weakly bonded molecular species is ADD. As will be discussed below, this feature is critical for efficient ionization of atmospheric species on the QLL. Figure 4 shows simulated infrared (IR) adsorption spectra, obtained with the dipole-moment autocorrelation method, similar to ref 62. Systems considered are the QLL, the crystalline ice surface (ICE), and the vapor−liquid interface (VLI). The hydrogen-bonded region of the QLL spectra is blue-shifted to the ICE and red-shifted to the VLI spectra. Blueshifting of the QLL spectra with respect to the ideal ice surface stems from the weaker hydrogen bonds of the QLL. Although not directly comparable to the present simulation, sumfrequency generation (SFG) spectra of the ice basal plane feature a small blue shift with increasing temperature.11,12,63 DNO3 Ionization on the QLL. Nitric acid requires several solvating OH bonds in order to release its proton. Even when placed in liquid watera dynamic and flexible environment where the anion can “sample” a large amount of positions and hydrogen-bonding environments, observing ionization in the picosecond time scale can be challenging.64,65 We may trace this behavior to the fact that nitric acid has a large, polyatomic anion, while its oxygen lone pair electrons are placed in a directional manner, at the tips of the N−O bonds (see Figure 7). It can then be difficult to “fit in” the nitrate ion NO3− into an environment of solvating OH bonds. This situation can be E

dx.doi.org/10.1021/jp505627n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

Figure 6. (a−c) Snapshots from an ab initio MD run for the nitric acid molecule on the QLL model surface. (d) Deuterium iodide after forming an ion pair separated by two water molecules on the QLL.

Figure 7. Left upper panel: number of hydrogen bonds formed by the nitric acid molecule as a function of MD simulation time. The number of hydrogen bonds is calculated using eq 3 of the Supporting Information, with d = 0.325 Å and averaged over 6 fs time windows. The Wannier function centers being analyzed are illustrated in the right panel. Left lower panel: oxygen−proton distances analyzed as a function of MD simulation time, averaged over 6 fs time windows. For the definition of h − l, see the right panel. In detail, h − l > 0 implies intact nitric acid, h − l ≈ 0 the Zundel cation, while h − l < 0 the Eigen cation. Right panel: oxygen, nitrogen, and hydrogen atoms are denoted by red, blue, and white spheres, respectively. Black spheres are the Wannier function centers used in the hydrogen-bond analysis.

acid became highly solvated. Furthermore, once the ionization takes place, there is not a clear distinction between a CIP (contact ion pair) and SSIP (solvent-separated ion pair) as the proton “rattles” nearby the nitrate ion. It fluctuates constantly between the Zundel and Eigen cations,66 while the fluctuating Zundel/Eigen pair can migratevia the Grothuss mechanism67further away from the nitrate anion, as depicted in the right-hand panel of Figure 5. Snapshots along an MD trajectory are illustrated in Figure 6a−c, where the number of hydrogen bonds N solvating the nitric acid (before the proton jump) is seen to fluctuate between N = 2 and N = 3. A more detailed analysis on the number of the solvating hydrogen bonds is given in Figure 7, where we can correlate the ionization and proton migration events to situations where the number of solvating hydrogen bonds to initiate the ionization is N ≥ 2. In Figure 7, the hydrogen-bond criteria is based on the Wannier function shifts (described in the Supporting Information), which makes the process of ionization more evident: shifting the lone-pair electrons away from the oxygen atoms makes the anion more positive, driving away the positive proton. The ionization pathway and events 1−3 discussed above are similar to those described by Bianco, Wang, and Hynes65 for nitric acid in liquid water systems. In contrast to liquid water, in the present case, the QLL structure has a high relative concentration of weakly bonded ADD species, their origin being in the spontaneously formed ice crystalline point defects of the QLL (as discussed above). This facilitates the type of reaction path described in Figure 5.

In the present case, weakly bonded water molecule species are a result of the surface crystalline point defects. Weak and reactive, yet static hydrogen-bond structures in ice-like systems can also be created with surface mismatch: an example of such a system is the water monolayer on a hydroxylated silica surface.68 DI Ionization on the QLL. Ionization of DI was studied in a similar fashion as in the case of DNO3, i.e., placing DI in close contact with the equilibrated QLL structure and letting the composite system evolve in an NVE ensemble MD simulation. In detail, DI was placed in a planar geometry, at ≈4 Å from the layer v molecules, at four different sites. Table 1 gives the time Table 1. Summary of Results for DI in Contact with the QLL: the Total Simulation Time (t) and the Time That It Took for DI to Ionize and Form a Solvent-Separated Ion Pair (tdis) trajectory

t (ps)

tdis (ps)

A B C D

14 14 14 10

4.2 2.8 3.3 1.4

it took for DI to form the contact ion pair: rapid ionization within a few picoseconds took place in all four trajectories considered, while the solvent-separated ion pair succeeded immediately the contact ion pair formation. During each trajectory, iodide steadily breaks an increasing number of the weak hydrogen bonds of the QLL, finally obtaining maximum F

dx.doi.org/10.1021/jp505627n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

solvation. A typical geometry, from trajectory “A”, after ionization is depicted in Figure 6d. In all cases, the iodide becomes embedded into layer v.

ADD species was also observed on the QLL. This feature is desirable for efficient proton hopping and the formation of the Eigen hydronium species. The QLL features slow self-diffusion, leading to a situation where anions are “anchored” on the surface; they can further react with incoming atmospheric species and lead, for example, to the production of chlorine reservoir species.71 The dynamical breaking and forming of water molecule hydrogen bonds is inhibited on the QLL (thus the formation of stable crystalline point defects), which can also facilitate the interaction of hydrophobic species with the surface. These features can be contrasted to completely “ice-like” and “liquidlike” surfaces: a perfect and crystalline ice surface lacks crystalline point defects, and its hydrogen bonds are strong. A liquid surface, on the other hand, while exhibiting weaker hydrogen bonds, lacks long-lived crystalline defects. The present work emphasizes the importance of the ice quasi-liquid layer for atmospheric chemistry. The model developed in the present work facilitates future ab initio studies of chemical reactions on the quasi-liquid layer.

■

CONCLUSIONS Reactivity of ice surfaces for acid ionization plays an important role both in astrochemistry69 and atmospheric25 chemistry, while low temperature (≈ −250 °C) laboratory studies of acid ionization with substrate-grown surfaces and ice nanoparticles give seemingly conflicting results for reactivity for hydrochloric acid.26−28,31−34 While such experiments26−28,31−34 used various model systems for studying ice reactivity, the surface reactivity may well be sensitive to the ice system under study.28 In particular, substrate-grown ice comes in various forms,36 namely, crystalline,28,32 sparse amorphous, and dense amorphous,31,32 while small ice nanoparticles26,27 have surface curvature, not necessarily present in the substrate-grown systems. The extent of acid ionization on low-temperature ice surfaces then still poses open questions, possibly depending on the particular ice system under study (there is no “universal” ice surface) and on the interpretation of the experiments.33,34 Ice under atmospheric conditions comes dominantly in the crystalline form (ice-Ih), exhibiting crystalline facets that can interact with impinging acid molecules: cirrus cloud17,70 ice particles, in particular, inhabit the upper stroposphere at temperatures as cold as ≈ −70 °C, creating then a vast surface area for ice/acid interactions. In the present work, we studied the reactivity of the crystalline ice-Ih basal plane. The so-called “quasi-liquid layer” (QLL)2,3,8 forms on this surface starting at ≈ −70 °C, while this should make the ice-Ih basal plane QLL a prominent system in the atmosphere. In earlier, classical force field molecular dynamics studies,16,40−44 the basal plane QLL was observed to maintain the bulk ice bilayer structure, while “partially melted” top layers form a characteristic density profile normal to the surface. In order to stabilize this structure, nanosecond time scale molecular dynamics simulations are mandatory. We confirmed this structure by a combined classical force field and densityfunctional theory study, while the following interpretation of the QLL structure was emphasized: crystalline point defects are formed spontaneously on the topmost bilayer during the nanosecond time scale molecular dynamics. The stability of this structure in the picosecond time scale was established by oxygen atom position histograms (both normal and parallel to the surface). The structural stability is also manifested by the small self-diffusion coefficient of the QLL when compared to liquid water, while, at the same time, orders of magnitude larger than that of crystalline ice.52,53 Theoretical studies, on the other hand, suggest that crystalline point defects and ledge sites enhance the reactivity of crystalline ice surfaces for acid ionization.37−39 We studied this aspect of the QLL with density functional theory molecular dynamics, by placing deuterated acids, namely, DNO3 and DI, in close contact with the equilibrated QLL model system featuring crystalline point defects. Picosecond time scale ionization for both acids was demonstrated. Ionization mechanisms followed pathways suggested by earlier studies, while the crystalline point defects on the QLL were playing a fundamental role in the ionization process; removing water molecules from crystalline ice results in broken, but at the same time stable, hydrogen-bonding structure. Its hydrogen bonds are weak, while impinging acids can break them in order to solvate properly. A high density of

■

ASSOCIATED CONTENT

* Supporting Information S

Details on simulation times and properties of the simulated vapor−liquid and ice−vacuum interfaces: mass profiles, oxygen−oxygen radial distribution functions, diffusion coefficient calculations, etc. Literature review of quasi-liquid layer calculations. Definition of the hydrogen bond, based on Wannier function centers. Details on acid ionization simulation runs on the quasi-liquid layer. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: sampsa.riikonen@iki.fi. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Finland’s IT Center for Science Ltd. (CSC) for the use of its computational resources. All authors acknowledge support by the Academy of Finland through the FiDiPro and Lastu programs, while R.B.G. acknowledges additionally Israel Science Foundation (grant 172/12) and US National Science Foundation (grant 0909227). S.R. acknowledges useful discussions with Audrey Dell Hammerich, Barak Hirshberg, and Garold Murdachaew. P.P. acknowledges University of Helsinki for a graduate school position.

■

REFERENCES

(1) Faraday, M. Note on regelation. Proc. R. Soc. London 1859, 10, 440−450. (2) Dash, J. G.; Fu, H.; Wettlaufer, J. S. The premelting of ice and its environmental consequences. Rep. Prog. Phys. 1995, 58, 115. (3) Li, Y.; Somorjai, G. A. Surface premelting of ice. J. Phys. Chem. C 2007, 111, 9631−9637. (4) Materer, N.; Starke, U.; Barbieri, A.; van Hove, M.; Somorjai, G. A.; Kroes, G. J.; Minot, C. Molecular surface structure of a lowtemperature ice Ih(0001) crystal. J. Phys. Chem. 1995, 99, 6267−6269. (5) Materer, N.; Starke, U.; Barbieri, A.; Hove, M. V.; Somorjai, G.; Kroes, G.-J.; Minot, C. Molecular surface structure of ice(0001): dynamical low-energy electron diffraction, total-energy calculations and molecular dynamics simulations. Surf. Sci. 1997, 381, 190−210.

G

dx.doi.org/10.1021/jp505627n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

(6) Bluhm, H.; Ogletree, D. F.; Fadley, C. S.; Hussain, Z.; Salmeron, M. The premelting of ice studied with photoelectron spectroscopy. J. Phys.: Condens. Matter 2002, 14, L227. (7) Dash, J. G. History of the search for continuous melting. Rev. Mod. Phys. 1999, 71, 1737−1743. (8) Petrenko, V. F.; Whitworth, R. W. Physics of ice; Oxford University Press: New York, 1999. (9) Lied, A.; Dosch, H.; Bilgram, J. H. Surface melting of ice Ih single crystals revealed by glancing angle x-ray scattering. Phys. Rev. Lett. 1994, 72, 3554−3557. (10) Dosch, H.; Lied, A.; Bilgram, J. Glancing-angle x-ray scattering studies of the premelting of ice surfaces. Surf. Sci. 1995, 327, 145−164. (11) Wei, X.; Miranda, P. B.; Shen, Y. R. Surface vibrational spectroscopic study of surface melting of ice. Phys. Rev. Lett. 2001, 86, 1554−1557. (12) Wei, X.; Miranda, P. B.; Zhang, C.; Shen, Y. R. Sum-frequency spectroscopic studies of ice interfaces. Phys. Rev. B 2002, 66, 085401. (13) Dosch, H.; Lied, A.; Bilgram, J. H. Disruption of the hydrogenbonding network at the surface of Ih ice near surface premelting. Surf. Sci. 1996, 366, 43−50. (14) Elbaum, M.; Lipson, S.; Dash, J. Optical study of surface melting on ice. J. Cryst. Growth 1993, 129, 491−505. (15) Furukawa, Y.; Yamamoto, M.; Kuroda, T. Ellipsometric study of the transition layer on the surface of an ice crystal. J. Cryst. Growth 1987, 82, 665−677. (16) Conde, M. M.; Vega, C.; Patrykiejew, A. The thickness of a liquid layer on the free surface of ice as obtained from computer simulation. J. Chem. Phys. 2008, 129, 014702. (17) Liou, K.-N. Influence of cirrus clouds on weather and climate processes: A global perspective. Mon. Weather Rev. 1986, 114, 1167− 1199. (18) Wang, P.-H.; Minnis, P.; McCormick, M. P.; Kent, G. S.; Skeens, K. M. A 6-year climatology of cloud occurrence frequency from stratospheric aerosol and gas experiment II observations (1985− 1990). J. Geophys. Res.: Atmos. 1996, 101, 29407−29429. (19) Wylie, D. P.; Menzel, W. P. Eight years of high cloud statistics using HIRS. J. Clim. 1999, 12, 170−184. (20) Solomon, S.; Garcia, R. R.; Rowland, F. S.; Wuebbles, D. J. On the depletion of antarctic ozone. Nature 1986, 321, 755−758. (21) Solomon, S. Stratospheric ozone depletion: a review of concepts and history. Rev. Geophys. 1999, 37, 275−316. (22) Dominé, F.; Shepson, P. B. Air-snow interactions and atmospheric chemistry. Science 2002, 297, 1506−1510. (23) Zhou, X.; Beine, H. J.; Honrath, R. E.; Fuentes, J. D.; Simpson, W.; Shepson, P. B.; Bottenheim, J. W. Snowpack photochemical production of HONO: A major source of OH in the arctic boundary layer in springtime. Geophys. Res. Lett. 2001, 28, 4087−4090. (24) Dibb, J. E.; Huey, L. G.; Slusher, D. L.; Tanner, D. J. Soluble reactive nitrogen oxides at south pole during ISCAT 2000. Atmos. Environ. 2004, 38, 5399−5409. (25) Huthwelker, T.; Ammann, M.; Peter, T. The uptake of acidic gases on ice. Chem. Rev. 2006, 106, 1375−1444. (26) Devlin, J. P.; Uras, N.; Sadlej, J.; Buch, V. Discrete stages in the solvation and ionization of hydrogen chloride adsorbed on ice particles. Nature 2002, 417, 269−271. (27) Buch, V.; Sadlej, J.; Aytemiz-Uras, N.; Devlin, J. P. Solvation and ionization stages of HCl on ice nanocrystals. J. Phys. Chem. A 2002, 106, 9374−9389. (28) Ayotte, P.; Marchand, P.; Daschbach, J. L.; Smith, R. S.; Kay, B. D. HCl adsorption and ionization on amorphous and crystalline H2O films below 50 K. J. Phys. Chem. A 2011, 115, 6002−6014. (29) Marchand, P.; Marcotte, G.; Ayotte, P. Spectroscopic study of HNO3 dissociation on ice. J. Phys. Chem. A 2012, 116, 12112−12122. (30) Marcotte, G.; Ayotte, P.; Bendounan, A.; Sirotti, F.; Laffon, C.; Parent, P. Dissociative adsorption of nitric acid at the surface of amorphous solid water revealed by x-ray absorption spectroscopy. J. Phys. Chem. Lett. 2013, 4, 2643−2648. (31) Kang, H.; Shin, T.-H.; Park, S.-C.; Kim, I. K.; Han, S.-J. Acidity of hydrogen chloride on ice. J. Am. Chem. Soc. 2000, 122, 9842−9843.

(32) Parent, P.; Lasne, J.; Marcotte, G.; Laffon, C. HCl adsorption on ice at low temperature: a combined x-ray absorption, photoemission and infrared study. Phys. Chem. Chem. Phys. 2011, 13, 7142−7148. (33) Devlin, J. P.; Kang, H. Comment on “HCl adsorption on ice at low temperature: a combined X-ray absorption, photoemission and infrared study” by P. Parent, J. Lasne, G. Marcotte and C. Laffon. Phys. Chem. Chem. Phys. 2012, 14, 1048−1049. (34) Parent, P.; Lasne, J.; Marcotte, G.; Laffon, C. Reply to the ‘Comment on “HCl adsorption on ice at low temperature: a combined X-ray absorption, photoemission and infrared study”’ by J. P. Devlin and H. Kang. Phys. Chem. Chem. Phys. 2012, 14, 1050−1053. (35) Devlin, J. P.; Buch, V. Surface of ice as viewed from combined spectroscopic and computer modeling studies. J. Phys. Chem. 1995, 99, 16534−16548. (36) Stevenson, K. P.; Kimmel, G. A.; Dohnálek, Z.; Smith, R. S.; Kay, B. D. Controlling the morphology of amorphous solid water. Science 1999, 283, 1505−1507. (37) Bolton, K.; Pettersson, J. B. C. Ice-catalyzed ionization of hydrochloric acid. J. Am. Chem. Soc. 2001, 123, 7360−7363. (38) Bolton, K. A QM/MM study of HCl adsorption at ice surface defect sites. J. Mol. Struct.: THEOCHEM 2003, 632, 145−156. (39) Riikonen, S.; Parkkinen, P.; Halonen, L.; Gerber, R. B. Ionization of nitric acid on crystalline ice: the role of defects and collective proton movement. J. Phys. Chem. Lett. 2013, 4, 1850−1855. (40) Kroes, G.-J. Surface melting of the (0001) face of TIP4P ice. Surf. Sci. 1992, 275, 365−382. (41) Picaud, S. Dynamics of TIP5P and TIP4P/ice potentials. J. Chem. Phys. 2006, 125, 174712. (42) Bishop, C. L.; Pan, D.; Liu, L. M.; Tribello, G. A.; Michaelides, A.; Wang, E. G.; Slater, B. On thin ice: surface order and disorder during pre-melting. Faraday Discuss. 2009, 141, 277−292. (43) Neshyba, S.; Nugent, E.; Roeselová, M.; Jungwirth, P. Molecular dynamics study of ice-vapor interactions via the quasi-liquid layer. J. Phys. Chem. C 2009, 113, 4597−4604. (44) Pfalzgraff, W.; Neshyba, S.; Roeselova, M. Comparative molecular dynamics study of vapor-exposed basal, prismatic, and pyramidal surfaces of ice. J. Phys. Chem. A 2011, 115, 6184−6193. (45) Neshyba, S. P.; Lowen, B.; Benning, M.; Lawson, A.; Rowe, P. M. Roughness metrics of prismatic facets of ice. J. Geophys. Res.: Atmos. 2013, 118, 3309−3318. (46) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269−6271. (47) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Comput. Phys. Commun. 2005, 167, 103−128. (48) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38, 3098−3100. (49) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements h-pu. J. Chem. Phys. 2010, 132, 154104. (50) Tobias, D. J.; Martyna, G. J.; Klein, M. L. Molecular dynamics simulations of a protein in the canonical ensemble. J. Phys. Chem. 1993, 97, 12959−12966. (51) VandeVondele, J.; Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. J. Chem. Phys. 2007, 127, 114105. (52) Toubin, C.; Picaud, S.; Hoang, P. N. M.; Girardet, C.; Demirdjian, B.; Ferry, D.; Suzanne, J. Dynamics of ice layers deposited on mgo(001): quasielastic neutron scattering experiments and molecular dynamics simulations. J. Chem. Phys. 2001, 114, 6371−6381. (53) Mizuno, Y.; Hanafusa, N. Studies of surface properties of ice using nuclear magnetic resonance. J. Phys., Colloq. 1987, 48, C1-511− C1-517. (54) Toubin, C.; Hoang, P.; Picaud, S.; Girardet, C. Time study of pollutants at the surface of ice at 200 K. Chem. Phys. Lett. 2000, 329, 331−335. H

dx.doi.org/10.1021/jp505627n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

(55) Girardet, C.; Toubin, C. Molecular atmospheric pollutant adsorption on ice: a theoretical survey. Surf. Sci. Rep. 2001, 44, 159− 238. (56) Parkkinen, P.; Riikonen, S.; Halonen, L. Global minima of protonated water clusters (H2O)20H+ revisited. J. Phys. Chem. A 2012, 116, 10826−10835. (57) Parkkinen, P.; Riikonen, S.; Halonen, L. (H2O)20 water clusters at finite temperatures. J. Phys. Chem. A 2013, 117, 9985−9998. (58) Parkkinen, P.; Riikonen, S.; Halonen, L. Configurational entropy in ice nanosystems: tools for structure generation and screening. J. Chem. Theory Comput. 2014, 10, 1256−1264. (59) Kuo, I.-F. W.; Mundy, C. J. An ab initio molecular dynamics study of the aqueous liquid-vapor interface. Science 2004, 303, 658− 660. (60) Kuühne, T. D.; Pascal, T. A.; Kaxiras, E.; Jung, Y. New insights into the structure of the vapor/water interface from large-scale firstprinciples simulations. J. Phys. Chem. Lett. 2011, 2, 105−113. (61) Baer, M. D.; Mundy, C. J.; McGrath, M. J.; Kuo, I.-F. W.; Siepmann, J. I.; Tobias, D. J. Re-examining the properties of the aqueous vapor−liquid interface using dispersion corrected density functional theory. J. Chem. Phys. 2011, 135, 124712. (62) Pham, T. A.; Huang, P.; Schwegler, E.; Galli, G. First-principles study of the infrared spectra of the ice Ih(0001) surface. J. Phys. Chem. A 2012, 116, 9255−9260. (63) Groenzin, H.; Li, I.; Buch, V.; Shultz, M. J. The single-crystal, basal face of ice Ih investigated with sum frequency generation. J. Chem. Phys. 2007, 127, 214502. (64) Shamay, E. S.; Buch, V.; Parrinello, M.; Richmond, G. L. At the water’s edge: nitric acid as a weak acid. J. Am. Chem. Soc. 2007, 129, 12910−12911 (PMID: 17915872). (65) Wang, S.; Bianco, R.; Hynes, J. T. An atmospherically relevant acid: HNO3. Comput. Theor. Chem. 2011, 965, 340−345. (66) Miller, Y.; Gerber, R. B. Dynamics of proton recombination with NO3- anion in water clusters. Phys. Chem. Chem. Phys. 2008, 10, 1091−1093. (67) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. The nature of the hydrated excess proton in water. Nature 1999, 397, 601− 604. (68) Murdachaew, G.; Gaigeot, M.-P.; Halonen, L.; Gerber, R. B. Dissociation of HCl into ions on wet hydroxylated (0001) α-quartz. J. Phys. Chem. Lett. 2013, 4, 3500−3507. (69) Hama, T.; Watanabe, N. Surface processes on interstellar amorphous solid water: adsorption, diffusion, tunneling reactions, and nuclear-spin conversion. Chem. Rev. 2013, 113, 8783−8839. (70) Baran, A. J. On the scattering and absorption properties of cirrus cloud. J. Quant. Spectrosc. Radiat. Transfer 2004, 89, 17−36. (71) Hammerich, A. D.; Finlayson-Pitts, B. J.; Gerber, R. B. NOX reactions on aqueous surfaces with gaseous hcl: formation of a potential precursor to atmospheric cl atoms. J. Phys. Chem. Lett. 2012, 3, 3405−3410.

I

dx.doi.org/10.1021/jp505627n | J. Phys. Chem. A XXXX, XXX, XXX−XXX