A Molecular Dynamics Study of the Long-Time Ice Ih Surface Dynamics

the water molecules in the top two to three bilayers is substantially faster than ... where molecules must move to the upper bilayers before diffusion...
0 downloads 0 Views 69KB Size
Download PDF
1590

J. Phys. Chem. B 2000, 104, 1590-1595

A Molecular Dynamics Study of the Long-Time Ice Ih Surface Dynamics Kim Bolton*,†,‡ and Jan B. C. Pettersson†,§ School of Engineering, UniVersity of Borås, SE-501 90, Borås, Sweden, Department of Chemistry, Physical Chemistry, Go¨ teborg UniVersity, SE-412 96, Go¨ teborg, Sweden, and School of EnVironmental Sciences, Go¨ teborg UniVersity, SE-412 96, Go¨ teborg, Sweden ReceiVed: September 29, 1999; In Final Form: December 8, 1999

Molecular dynamics simulations of the ice Ih surface between 180 and 210 K showed that the dynamics of the water molecules in the top two to three bilayers is substantially faster than that of the molecules in the lower (bulk) bilayers. Within the simulation time of tens of nanoseconds, there is rapid exchange of molecules between these upper bilayers, but not between these and the lower bilayers. In-plane translation of the molecules in the top surface bilayer leads to rapid surface reconstruction, and a structure consisting primarily of water heptagons, hexagons, and pentagons replaces the crystal geometry. An Arrhenius analysis of the diffusion of molecules in the top bilayer yields a barrier for in-plane diffusion E0 ) 23.2 ( 2.9 kJ mol-1 and a preexponential factor D0 ) 0.77 cm2 s-1 (ln(D0) ) -0.26 ( 1.99). The activation barrier, which is the same as the experimentally measured value, is similar to the energy required to break a hydrogen bond and supports a diffusion mechanism where water molecules move by repeated breaking and formation of hydrogen bonds. The similarity between the simulated and experimental barrier heights supports an in-plane diffusion mechanism where molecules must move to the upper bilayers before diffusion can occur.

1. Introduction Ice surfaces play an important role in many areas of environmental science.1 An example is the production of chlorine on Type II polar stratospheric clouds (PSCs).2,3 Although the detailed mechanism is not known, the ice surface plays a catalytic role in producing Cl2 from chlorine-containing molecules such as HCl, HOCl, and ClONO2. These chlorine molecules photolyze to chlorine radicals which react with and deplete stratospheric ozone. Ice-catalyzed production of chlorine comprises many steps.2 These include adsorption of chlorine-containing compounds onto the surface, ionization, molecular and ionic diffusion on and into the ice surface, reaction between surface molecules and ions, and desorption. The nature of the PSC ice surface is expected to play an important role in the mechanism and kinetics of the individual steps, and hence also in the overall reaction. Under stratospheric conditions the ice surface is dynamic, with condensation and evaporation rates of 1015 molecules cm-2 s-1.4 This feature led Gertner and Hynes to propose a mechanism for ionization of stratospheric HCl, where the adsorbed molecule is trapped in the crystalline ice surface by the rapid deposition of water molecules.5,6 They estimate that ionization will occur in ≈ 1-10 ms,6 which is the time required to incorporate the HCl into the crystal lattice. Their mechanism assumes that the ice surface is crystalline and that the deposited HCl and H2O molecules adsorb at, or diffuse rapidly to, the crystal lattice sites. It has, however, been suggested that the ice surface is liquidlike at 190 K (which is relevant to stratospheric ice chemistry) in that there is substantial rotational mobility of the surface water molecules.7 This may enhance HCl ionization at the ice surface. * Author to whom correspondence should be addressed. Tel: +46-33174602. E-mail: [email protected]. † School of Engineering, University of Bora ˚ s. ‡ Department of Chemistry, Physical Chemistry, Go ¨ teborg University. § School of Environmental Sciences, Go ¨ teborg University.

There have been numerous theoretical and experimental studies of ice surface dynamics. In a classical molecular dynamics study of the TIP4P ice surface, Karim and Haymet8 observed surface melting at 270 K and a stable ice surface, within their simulation time of 160 ps, at 240 K. In another study Kroes observed TIP4P surface melting at temperatures larger than 230 K.9 At temperatures below 210 K no diffusional motion was observed within the simulation time, although there was substantial rotational motion and significant surface reconstruction. Water molecules in the uppermost surface bilayer formed heptagon and pentagon structures at temperatures as low as 190 K. For surface temperatures below 210 K an upper limit of 10-6 cm2 s-1 was determined for both in-plane, Dxy, and outof-plane, Dz, diffusion. Early calculations based on a rigid ice surface with an effective central force H2O-H2O potential of Stillinger and Rahman gave a barrier height for surface diffusion of 10.6 kJ mol-1,10 and more recent calculations based on a flexible TIP4P surface yielded barrier heights between 15.4 and 27.0 kJ mol-1.11 Experimental studies of ice surface dynamics include determination of the rates of evaporation, condensation, surface diffusion, and bulk diffusion.4,12,13 George and co-workers measured bulk diffusion coefficients, Dz, that ranged from 2.2 × 10-16 to 3.9 × 10-14 cm2 s-1 for temperatures between 153 and 170 K.14-16 An Arrhenius analysis yielded a bulk diffusion barrier of 71.4 kJ mol-1. They also determined an upper limit of 10-9 cm2 s-1 for Dxy at 140 K.15,16 However, they note that since there is fairly rapid diffusion of water molecules into the ice surface their Dxy is not solely for molecules in the uppermost bilayers, but also for molecules that penetrate into the lattice. Since surface penetration may retard in-plane diffusion, molecules in the top bilayers may have a larger diffusion coefficient than their experimental value.13 Nuclear magnetic resonance (NMR) studies for ice temperatures between 253 and 273 K yielded in-plane diffusion coefficients between 1.35 × 10-9 and

10.1021/jp9934883 CCC: $19.00 © 2000 American Chemical Society Published on Web 02/01/2000

Long-Time Ice Ih Surface Dynamics 3.08 × 10-9 cm2 s-1 and a barrier of 23.2 kJ mol-1.17 This barrier is in the range determined for TIP4P ice,11 and is similar to the barrier for self-diffusion in liquid water, which is 19.3 kJ mol-1.18 Comparison between simulation and experimental data for ice nanocrystals indicate that there is significant disordering in several bilayers at temperatures as low as 2583 K, and that adsorption of molecules such as H2S increases the ordering in these bilayers.19 In this work molecular dynamics simulations were used to investigate the dynamics of water molecules at the ice Ih surface over 48 ns. In and out-of-plane motion were investigated for ice temperatures of 180-210 K. In the next section the potential energy surface (PES) and trajectory methods are described. In Section 3, the simulation results are discussed and compared to existing experimental and theoretical data. Reasons for some of the apparent discrepancies9,13 between experimental and theoretical results are suggested, and the relevance of the surface dynamics to stratospheric chemistry is discussed. The conclusions are presented in Section 4. 2. Methods The PES and trajectory methods have been detailed in a previous paper.20 The TIP4P model,21 which assumes rigid water molecules, was used to describe the H2O-H2O intermolecular interactions. This model has been used to simulate, for example, ice surfaces,22 ice-liquid interfaces,23 bulk water,24 and bulk crystalline25 and amorphous ice.26 A switching function was used to smoothly adjust the intermolecular forces to zero for H2OH2O separations larger than 10 Å.20 An infinitely large ice Ih surface was simulated using a 27 × 31 × 29 Å3 slab with periodic boundary conditions in the surface xy plane. The slab had eight bilayers with 96 molecules in each bilayer. The initial (unthermalized) crystal structure was based on the unit cell geometry determined by Hayward and Reimers, which accounts for proton disorder.20,27 The positions of the water molecules in the bottom two bilayers were held fixed to preserve the bulk crystalline structure, and the molecules in the top six bilayers were thermalized by coupling to a temperature bath. This was done by integrating the molecules in the top five bilayers using the Verlet method28 and using the Langevin equation29 to integrate the molecules in the sixth bilayer. This incorporates frictional and stochastic energy exchange with the “constant temperature bulk” into the dynamics. For each of the temperatures investigated the thermalized surfaces were generated from the crystalline ice Ih surface by propagating the molecular coordinates at the desired temperature for 0.5-1.0 ns. This was sufficient time for thermalization of the translational and rotational velocities and to allow for reconstruction of the ice surface. The effect of slab thickness on the ice surface dynamics was tested by simulating a twelvebilayer slab at 210 K. The similarity of these results and those obtained for the eight-bilayer surface confirmed that the thinner surface could be used to simulate ice surface properties. Surfaces at 180, 190, 200, and 210 K were simulated by propagating a single trajectory for each temperature. Each trajectory was propagated for 48 ns which required ≈ 3 weeks on a CRAY ORIGIN 2000 computer. As in the previous study,20 a 6 fs integration time step was used. The dynamics of the top three bilayers was analyzed. Molecules in the lower bilayers did not have large amplitude motion within 48 ns and are not included in the analysis presented below. The z (out-of-plane) center-of-mass coordinates were monitored every 10 ps and plotted in histograms. The histograms show the probability of finding water molecules at various depths below the ice surface.

J. Phys. Chem. B, Vol. 104, No. 7, 2000 1591 The in-plane diffusion coefficient, Dxy,30 was obtained for each molecule from

Dxy )

MSD 4t

(1)

The mean square displacement, MSD, over the time interval t is tf-t

MSD )

[(x(ti + t) - x(ti))2 + (y(ti + t) - y(ti))2] ∑ t )0 i

N

(2)

where the summation is over times ti ranging from 0 to tf - t, tf is the final time that the x and y molecular positions are monitored (i.e., 48 ns) and N is the number of terms included in the summation. Determining Dxy from the slope of the MSD over the first 16 ns (one-third of the trajectory) ensured statistical convergence. Error bars for ensemble averaged Dxy (i.e., averaged over all molecules) were determined from the bootstrap method.31,32 3. Results and Discussion The melting characteristics of TIP4P ice Ih have been studied by several groups.8,9,23,33 As mentioned in the Introduction, MD simulations by Karim and Haymet suggest a surface melting temperature less than 270 K, while a stable surface structure was observed at 240 K for their simulation time of 160 ps.8 In agreement with these studies, Kroes showed that, within simulation times of ≈ 100-300 ps, there was no calculable self-diffusion below 210 K, indicating that the surface is solidlike over these times.9 Reorientational correlation times are, however, closer to that of liquid than bulk ice. Nada and Furukuwa observed crystallization of TIP4P(0001) liquid water in contact with ice within 1 ns at 230 K.23 This suggests a melting point larger than 230 K. On the other hand, recent free energy studies of Vlot et al. yield a melting point of 214 ( 6 K.33 Although the calculated melting points depend on the theoretical method and/or analysis technique employed, all studies suggest a temperature larger than 210 K. The present investigations also show that the ice does not melt (within 48 ns) at the temperatures investigated. As is discussed below, the ice structure does not show a continuous change over time that would be associated with melting. These simulations are, thus, of a solid-phase ice slab. However, the simulation (TIP4P) temperatures are probably not equal to actual (experimental) temperatures. MD results obtained at 180 K, for example, which is ≈ 35-70 K below the TIP4P melting point, may not be directly comparable to experimental results obtained at this temperature, but rather to results obtained at higher temperatures. Comparison between simulation and experiment is also complicated by other factors, e.g., the different time scales of experimental and simulation studies. Comparison between the present simulation results and experimental data is, thus, qualitative. 3.1. Out-of-Plane Dynamics. The probabilities of finding molecules at various depths below the ice surface are shown in Figure 1. The histograms are similar to those obtained by Kroes for shorter trajectories.9 A distinct bilayer structure exists at all temperatures, but there is a merging of adjacent bilayers with increasing temperature. The third bilayer does not merge with the fourth bilayer (not shown) at any of the temperatures, and at 180 K the second and third bilayers do not merge.

1592 J. Phys. Chem. B, Vol. 104, No. 7, 2000

Figure 1. Histograms of the water z (out-of-plane) coordinates for ice surfaces at 180, 190, 200, and 210 K. Only the top three bilayers are shown.

Figure 2. Out-of-plane motion (panel a) and in-plane x coordinate motion (panel b) of a surface water molecule for a surface temperature of 180 K. The dot-dashed line in panel a demarks the boundary between adjacent bilayers.

The out-of-plane motion of two surface water molecules, one at 180 K and the other at 210 K, are shown in Figures 2a and 3a, respectively. The regions of the two uppermost bilayers, defined from Figure 1, are also shown. Similar dynamics are observed for 190 and 200 K. The molecules move rapidly between adjacent bilayers, and the frequency of interbilayer crossing increases with temperature. This is consistent with the merging of the upper two bilayers seen in Figure 1. The molecule at 180 K is trapped in the second bilayer between ≈ 8 and 48 ns. Figures 2a and 3a show that molecules “hop” from one z position (or bilayer) to another, remain at the new z position for a relatively long time, and then hop to a new position. This indicates that out-of-plane translation occurs by repeated barrier crossings. As will be discussed below, this is probably due to repeated breaking and formation of hydrogen bonds (H-bonds). 3.2. In-Plane Dynamics. Under equilibrium thermal conditions the in-plane structure of the top bilayer is not the hexagonal crystal structure, but consists primarily of heptagons, hexagons, and pentagons. This surface reconstruction, which was reported by Kroes,9 occurs within 0.5 ns of the initial surface thermalization for all temperatures considered here. The surface did not return to the crystal structure at any stage during the trajectories, and when the energy of typical reconstructed

Bolton and Pettersson

Figure 3. Same as for Figure 2 but for a surface temperature of 210 K.

Figure 4. In-plane diffusion coefficients, Dxy, for all molecules in the top three bilayers at 210 K. Dxy is plotted as a function of the timeaveraged z position for each molecule.

surfaces was minimized they maintained their reconstructed geometries. The minimum energies of the reconstructed surfaces straddled the crystal structure minimum energy and were within 0.15 kJ mol-1 molecule-1 of the crystal structure energy. The TIP4P ice surface thus has multiple local energy minima. By restricting the motion of the surface molecules it was observed that the rapid reconstruction occurs only when a molecule in the first bilayer moves above the surface. This creates a deficit of molecules in the original crystal bilayer which subsequently reconstructs. Figures 2b and 3b show the time dependence of the x coordinates for the molecules plotted in Figures 2a and 3a, respectively. Translation is faster at the higher temperature and occurs via repeated hopping from one x position to another. The x positions of the molecules remain constant for a relatively long time between successive hops. Similar dynamics is seen in the y direction and at 190 and 200 K. This indicates that, similar to out-of-plane translation, in-plane motion occurs by repeated barrier crossings, probably due to successive breaking and formation of H-bonds. The in-plane diffusion coefficients, Dxy, for all molecules in the top three bilayers at 210 K are shown in Figure 4. The plot shows the dependence of Dxy on the time averaged z position of each molecule, i.e., the vertical position in the ice surface averaged over the 48 ns trajectory. Similar positive correlations between Dxy and z are obtained for the lower temperatures, but

Long-Time Ice Ih Surface Dynamics

Figure 5. Ln(Dxy) versus the inverse surface temperature.

Dxy decreases with temperature. As suggested by George and Livingston,13 molecules that are, on average, in the upper bilayers diffuse faster than those in the lower bilayers. An average in-plane diffusion coefficient was obtained for molecules in the first bilayer (zj g -2.4 Å in Figure 4) for each surface temperature. They are 1.7 × 10-7, 3.8 × 10-7, 7.9 × 10-7, and 1.6 × 10-6 cm2 s-1 for 180, 190, 200, and 210 K, respectively. Figure 5 is a plot of ln(Dxy) against the inverse temperature. The excellent fit of the simulated data to a straight line indicates that the same mechanism determines the diffusion rate at all temperatures in this range. An Arrhenius analysis yields an in-plane diffusion barrier E0 ) 23.2 ( 2.9 kJ mol-1 and a preexponential D0 ) 0.77 cm2 s-1 (ln(D0) ) -0.26 ( 1.99). Errors are for a 95% confidence interval. The barrier height is slightly less than the energy required to break a H-bond (26.1 kJ mol-1 for the TIP4P model),21 which supports the proposal that in-plane translation occurs via repeated breaking and formation of H-bonds. Comparison with PreVious Studies. The simulated in-plane diffusion barrier height is within the range determined by Batista and Jo´nsson (15.4-27.0 kJ mol-1)11 and is larger than the value determined from an effective central force model (10.6 kJ mol-1).10 It differs from the value of 17.4 kJ mol-1 obtained from the microcanonical simulations by Batista and Jo´nsson, but their value is for the motion of water molecules on a surface where the top bilayer molecules are in their crystal lattice sites, and our value is for molecules that are part of a reconstructed top bilayer. Their barrier may have more relevance at lower temperatures and where there is no reconstruction of the top surface bilayer. The present simulations yield the same barrier height as that determined by NMR studies for ice temperatures between 253 and 273 K.17 This indicates that the diffusion mechanism proposed in this work is valid at these higher ice temperatures. The simulated diffusion rates are consistent with Kroes’s upper bound of Dxy e 10-6 cm2 s-1 for temperatures below 210 K.9 To compare with experiment, the simulated E0 and D0 were used to calculate Dxy ) 1.9 × 10-9 cm2 s-1 at 140 K and Dxy > 10-4 cm2 s-1 for temperatures above 250 K. Although the value at 140 K is consistent with the upper bound of 10-9 cm2 s-1 measured by George and co-workers,15,16 it must be remembered that their measurement includes water molecules that have penetrated into the surface, whereas the simulated rate is for molecules that have remained, on average, in the first bilayer. The Dxy obtained from the NMR study, which are between 1.35 × 10-9 and 3.08 × 10-9 cm2 s-1 for temperatures

J. Phys. Chem. B, Vol. 104, No. 7, 2000 1593 between 253 and 273 K, are also for molecules that penetrated the surface. The dependence of Dxy on surface penetration, and on the difference in simulation (TIP4P) and actual temperature, prevents direct comparison between the simulated and experimental results, and the simulated diffusion rate is expected to be larger than the experimental value. The similarity between the experimental and simulated inplane diffusion barrier, which is temperature independent, indicates that repeated breaking and formation of H-bonds in the upper bilayer is the mechanism for in-plane diffusion. Thus molecules that are temporarily in the lower bilayers do not rapidly diffuse in the xy plane, but move to the top bilayer before in-plane translation occurs. This behavior is seen when comparing the data in panels a and b of Figures 2 and 3, respectively. At 180 K the molecule is in the second bilayer between 2-3 and 8-48 ns, and at 210 K the molecule repeatedly enters the second bilayer, e.g., between 19-25, 32-35 and 45-48 ns. During these times there is no overall translation in the x direction. The same is true for the y direction and for 190 and 200 K. Substantial in-plane surface diffusion only occurs when the molecules are in the top bilayer. The Arrhenius preexponential, D0 ) 0.77 cm2 s-1, reflects the frequency of possible diffusion events and is influenced by the total time that molecules are in the top bilayer. As shown in Figure 4, molecules that are, on average, in the third bilayer (zj e -6.4 Å) do not move to the top bilayer within 48 ps and have thus not diffused over long distances. Experiments are performed over far longer times and molecules may exchange between many bilayers. The total time that each molecule spends in the top bilayer depends on the frequency that it returns to this bilayer (which depends on the ice surface thickness and the out-of-plane diffusion rate) and the time of the experimental measurement. To illustrate the dependence of the in-plane diffusion preexponential factor on the ice thickness, we consider an ice surface consisting of Nb bilayers, where there is sufficient time for all surface molecules to repeatedly return to the top bilayer so that in-plane diffusion can occur. Assuming that molecules spend, on average, the same time in each bilayer (so that there is no accumulation of molecules in any bilayer), the probability that a molecule is in the top bilayer is

P)

1 Nb

(3)

and, based on the simulation result, the preexponential factor for in-plane diffusion is

D0 )

1 × 0.77 cm2 s-1 Nb

(4)

Based on Brown and George’s15 out-of-plane diffusion barrier height (69.5 kJ mol-1) and preexponential (9.7 × 107 cm2 s-1), the out-of-plane diffusion rate at 140 K is 8.3 × 10-19 cm2 s-1. Hence, molecules diffuse just 7.7 Å in their experimental time of 1 h. This indicates that their measured in-plane diffusion rate should be for molecules that are, on average, in the uppermost surface bilayers. This would yield an experimental diffusion rate (determined to be less than 10-9 cm2) that is only slightly lower than that extrapolated from the simulated rates (1.9 × 10-9 cm2), given that the effect of the different simulation (TIP4P) and experimental temperatures is small. 3.3. Model for Out-of-Plane Dynamics. The dynamics illustrated in Figures 1, 2a, and 3a suggests a simple model for out-of-plane motion. (It should be remembered that motion

1594 J. Phys. Chem. B, Vol. 104, No. 7, 2000

Bolton and Pettersson

TABLE 1: Probabilities, pij, for Molecules Moving from Bilayer i to j, and the Average Number of Molecules, Ni, in Bilayer ia p11) 0.95 p12) 0.08 a

p21) 0.06 p22 ) 0.91 p23 ) 0.03

p32 ) 0.03 p33 ) 0.97 p34 ) 0.00

N1 ) 92.0 N2 ) 100.4 N3 ) 95.6

The results are obtained from the 48 ns simulation at 210 K.

occurring on time scales longer than 48 nsse.g., bulk diffusions is not observed in the simulations and will not be included in the model.) Molecules are rapidly exchanged between the upper two to three bilayers, and there is no exchange with the lower (bulk) bilayers. If pij is the probability that a molecule in bilayer i, i ) 1,2,3, at time t moves to bilayer j, j ) 1,2,3, at a later time t + ∆t, then the number of molecules in bilayer i at t + ∆t is

Ni(t + ∆t) ) Ni(t) [1 ) Ni(t)pii +

pij] + ∑ Nj(t)pji ∑ i*j j*i

(5)

Nj(t)pji ∑ j*i

(6)

Equation 6 is obtained from eq 5 from the normalization criterion 3

pij ) 1, ∑ j)1

i ) 1,2,3

(7)

This model can readily be extended to many bilayers, or when bilayer and interbilayer regions can be distinguished. The probabilities, pij, were obtained from the simulations by determining the relative frequency that a molecule hops from bilayer i to j, i.e.,

pij )

no. of hops i f j no. of hops i f k, k ) 1,2,3

(8)

where the number of hops was summed over the 48 ns trajectory. Interbilayer hopping was monitored by obtaining the timeaveraged z position of each molecule over 60 fs windows, and determining if this average position changed between successive intervals (i.e., ∆t ) 60 fs in eq 5). Effects of short-lived hopss that would be associated with molecular vibrations rather than diffusionsare thus eliminated. The results presented below are not sensitive to moderate changes in ∆t. An ensemble-averaged pij, which was used in the model, was obtained by averaging over all molecules. The results for 210 K are shown in Table 1. The probabilities for remaining in a bilayer, pii, i ) 1,2,3, are larger than those for hopping between bilayers, and there is a larger probability for exchange between bilayers one and two than between two and three. There is no direct exchange between bilayers one and three, and the probability of hopping to the fourth bilayer is zero. Similar results were obtained for the lower surface temperatures although the probability of exchange between bilayers decreased with temperature. At 180 K there is no exchange between the upper two bilayers and bilayer three. The number of molecules in each bilayer, Ni, shown in the table is the time-averaged value obtained from the trajectory at 210 K. Similar numbers are obtained at lower temperatures. The number of molecules in the third bilayer is very similar to that for the crystal lattice, but four molecules have, on average, moved from bilayer one to two. Although this exchange may be energetically favorable, the tabulated numbers depend on the criterion used to distinguish molecules in bilayers one and two. Based on the histograms in Figure 1, a depth of -2.4 Å

Figure 6. Probabilities, pi, of finding a molecule, initially in the top bilayer, in bilayer i as a function of time. The results are for 210 K.

was used to delineate bilayers one and two. However, it is clear that there is substantial merging of these bilayers at all temperatures. The model can be used with the data in Table 1 to follow, for example, the out-of-plane motion of one molecule initially in the top bilayer at 210 K. The results, which show the ensemble averaged motion, are plotted in Figure 6. Due to the rapid exchange of molecules between bilayers at 210 K, a molecule initially in the top bilayer has a similar probability of being in any of the three bilayers within 10 ns. At long times a steady state is reached, as seen by the constant probabilities, pi, i ) 1,2,3. The number of molecules in each bilayer, Ni ) piNtot, where Ntot ) 288 is the total number of water molecules in the three bilayers, is thus constant. The fraction of molecules in each bilayer is similar to that determined directly from the MD simulations and listed in Table 1. 3.4. Relevance to Stratospheric Chemistry. The uncertainty in the difference between actual (e.g., experimental) and simulation (TIP4P) temperatures prevents extrapolation of the MD results to stratospheric temperatures. Nonetheless, it is expected that the rate of in- and out-of-plane diffusion in the outer surface bilayers will be larger than for the bulk (with both surface and bulk rates decreasing with temperature). Surface reconstruction, enhanced rotational motion of the surface molecules7 and the enhanced dynamics between the upper surface bilayers may affect surface-catalyzed heterogeneous chemistry occurring on PSCs. Previous studies showed that large surface heptagons and octagons enhance the rate of surface penetration of Ar atoms,34 and that penetration can occur after the Ar thermalizes on the surface. Recent simulations also showed that it is possible for molecular HCl to directly penetrate the top bilayer,35 which may enhance the efficiency of HCl adsorption and ionization. 4. Conclusions Molecular dynamics simulations of the TIP4P ice Ih surface between 180 and 210 K showed that the top surface bilayer reconstructs to a stable noncrystalline geometry, that in-plane diffusion of molecules that are, on average, in this bilayer is faster than that of molecules in lower bilayers, and that there is rapid exchange of molecules between the upper two to three bilayers. An Arrhenius analysis of the in-plane diffusion rate for the molecules that are, on average, in the top bilayer yielded a

Long-Time Ice Ih Surface Dynamics barrier height of 23.2 ( 2.9 kJ mol-1 and a preexponential factor of 0.77 cm2 s-1. The barrier height is consistent with experimental data at higher temperatures. Since the simulated and experimental barrier heights are the same, it is proposed that molecules must move to the upper surface layers before rapid in-plane diffusion can occur. This results in an inverse dependence of the in-plane diffusion rate on the ice thickness and/or experimental measurement time, and suggests a means of validating the proposed mechanism in future experiments. Acknowledgment. K.B. and J.P. are grateful for financial support from the Adlerbertska Forskningsfonden and the Swedish Natural Science Research Council, respectively. References and Notes (1) Robinson, G. W.; Zhu, S.-B.; Singh, S.; Evans, M. W. Water in Biology, Chemistry and Physics; World Scientific: River Edge, NJ, New Jersey, 1996. (2) Tabazadeh, A.; Turco, R. P. J. Geophys. Res. 1993, 98, 727. (3) Banham, S. F.; Horn, A. B.; Koch, T. G.; Sodeau, J. R. Faraday Discuss. 1995, 100, 321. (4) Haynes, D. R.; Tro, N. J.; George, S. M. J. Phys. Chem. 1992, 96, 8502. (5) Gertner, B. J.; Hynes, J. T. Science 1996, 271, 1563. (6) Gertner, B. J.; Hynes, J. T. Faraday Discuss. 1998, 110, 301. (7) Kroes, G.-J. Comments At. Mol. Phys. 1999, 34, 259. (8) Karim, O. A.; Haymet, A. D. J. J. Chem. Phys. 1988, 89, 6889. (9) Kroes, G.-J. Surf. Sci. 1992, 275, 365. (10) Hale, B. N.; Kiefer, J.; Ward, C. A. J. Chem. Phys. 1981, 75, 1991. (11) Batista, E. R.; Jo´nsson, H. in Proceedings of the IIX Workshop on Computational Materials Science; Ruggerone, P., Fiorentini, V., Eds.; to be published. (12) Chaix, L.; van den Bergh, H.; Rossi, M. J. J. Phys. Chem. A 1998, 102, 10300.

J. Phys. Chem. B, Vol. 104, No. 7, 2000 1595 (13) George, S. M.; Livingston, F. E. Surf. ReV. 1997, 4, 771. (14) Livingston, F. E.; Whipple, G. C.; George, S. M. J. Phys. Chem. B 1997, 101, 6127. (15) Brown, D. E.; George, S. M. J. Phys. Chem. 1996, 100, 15460. (16) Livingston, F. E.; Whipple, G. C.; George, S. M. J. Chem. Phys. 1998, 108, 2197. (17) Mizuna, Y.; Hanafusa, N. J. Phys. (suppl 3) 1987, 48, C1-511. (18) Eisenberg, D.; Kauzmann, W. The Structure and Properties of Water; Oxford University Press: New York, 1969. (19) Delzeit, L.; Devlin, M. S.; Rowland, B.; Devlin, J. P.; Buch, V. J. Phys. Chem. 1996, 100, 10076. (20) Bolton, K.; Svanberg, M.; Pettersson, J. B. C. J. Chem. Phys. 1999, 110, 5380. (21) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 29, 926. (22) Materer, N.; Starke, U.; Barbieri, A.; Van Hove, M. A.; Somorjai, G. A.; Kroes, G.-J.; Minot, C. J. Phys. Chem. 1995, 99, 6267. Materer, N.; Starke, U.; Barbieri, A.; Van Hove, M. A.; Somorjai, G. A.; Kroes, G.-J.; Minot, C. Surf. Sci. 1997, 381, 190. (23) Nada, H.; Furukawa, Y. J. Phys. Chem. B 1997, 101, 6163. (24) Poole, P. H.; Sciortino, F.; Essmann, U.; Stanley, H. E. Phys. ReV. E 1993, 48, 3799. (25) Tanaka, H. J. Chem. Phys. 1998, 108, 4887. (26) Wilson, M. A.; Pohorille, A.; Jenniskens, P.; Blake, D. F. Orig. Life EVol. Biosp. 1995, 25, 3. (27) Hayward J. A.; Reimers, J. R. J. Chem. Phys. 1997, 106, 1518. (28) Svanberg, M. Mol. Phys. 1997, 92, 1085. (29) DePristo, A. E.; Metiu, H. J. Chem. Phys. 1989, 90, 1229. (30) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1990. (31) Efron, B. SIAM ReV. 1979, 21, 460. (32) Diaconis, P.; Efron, B. Sci. Am. 1993, 248, 96. (33) Vlot, M. J.; Huinink, J.; van der Eerden, J. P. J. Chem. Phys. 1999, 110, 55. (34) Bolton, K.; Pettersson, J. B. C. Chem. Phys. Lett., accepted. (35) Al-Halabi, A.; Kleyn, A. W.; Kroes, G. J. Chem. Phys. Lett. 1999, 307, 505.