Molecular Perspective on Water Vapor Accommodation into Ice and Its Dependence on Temperature

Accommodation of vapor-phase water molecules into ice crystal surfaces is a fundamental process controlling atmospheric ice crystal growth. Experimental studies investigating the accommodation process with various techniques report widely spread values of the water accommodation coefficient on ice, αice, and the results on its potential temperature dependence are inconclusive. We run molecular dynamics simulations of molecules condensing onto the basal plane of ice Ih using the TIP4P/Ice empirical force field and characterize the accommodated state from this molecular perspective, utilizing the interaction energy, the tetrahedrality order parameter, and the distance below the instantaneous interface as criteria. Changes of the order parameter turn out to be a suitable measure to distinguish between the surface and bulk states of a molecule condensing onto the disordered interface. In light of the findings from the molecular dynamics, we discuss and re-analyze a recent experimental data set on αice obtained with an environmental molecular beam (EMB) setup [KongX.; J. Phys. Chem. A2014, 118 ( (22), ), 3973−397924814567] using kinetic molecular flux modeling, aiming at a more comprehensive picture of the accommodation process from a molecular perspective. These results indicate that the experimental observations indeed cannot be explained by evaporation alone. At the same time, our results raise the issue of rapidly growing relaxation times upon decreasing temperature, challenging future experimental efforts to cover relevant time scales. Finally, we discuss the relevance of the water accommodation coefficient on ice in the context of atmospheric cloud particle growth processes.


Evaluating flux models with molecular beam data
The two flux models discussed in the main text -one accounting solely for evaporative losses from the surface layer (1P), and a second also accounting for bulk accommodation surface layer desorption (2P) -were evaluated against the molecular beam data collected by Kong et al. (2014). Both 1P and 2P contain a fitted pre-factor A in the interval [0, 1000], while 2P includes also a fitted bulk accommodation rate " in the interval [0, 500] molecules/s. Examples of the detected intensity of evaporated D2O as modelled by these models for temperatures T = 170 K and T = 200 K are shown in Fig. S1 alongside their respective measurement time series. Figure S1: Evaporation only 1P (cyan) and evaporation and bulk accommodation 2P models (red) flux models for detected D2O intensity due to evaporation from the ice slab. Raw detected data is shown by the grey line and its averaged time series into 48 bins, to which the models are fit, are indicated by black dots. Panel a) is for temperature T = 170 K and panel b) for T = 200 K.  Figure S1 for all temperatures and including repeat fits for 1P and 2P flux models to: raw data (6000 timestamps) dark blue and orange, data binned into 48 bins yellow and purple, and 100 bins green and light blue, respectively. Figure S2 shows the modelled D2O detection intensity for all temperatures found by fitting the 1P and 2P frameworks to the various time-averaged binnings given in the legend. The associated fit statistics and fitted parameters are given in Fig. S3. Differences between 1P and 2P schemes are more pronounced at lower temperatures T = 170, 180, 182 and 184 K, and S4 negligible for higher temperatures. Upon visual inspection, the resolution of the temporal averaging has no impact on the modelled intensity or fit accommodation rates (Fig. S3a).
Though, as expected, this is apparent in the R 2 ( Fig. S3c) which increases with decreased bin number due to noise in the measurement data. This modelling was performed using a saturation vapour pressure of D2O over ice 1 as it returned the best R 2 statistics for the 48-bin procedure (see Fig. S4 for comparison with other saturation vapour pressures).   The 10 th , 25 th , 50 th , 75 th , and 90 th percentiles of the probability density of order parameter and short-range interaction energy values for the first layer (-3 Å < d <= 0 Å) and the third layer (-9 Å < d <= -6 Å) were calculated for all temperatures. Tables S1 and S2 contain the 10 th , 25 th , and median values for the first and third layer, respectively. Table S1: Definition of different "surface-accommodated" (AC) sub-spaces in E-q-d space. The values are presented for the top-most layer (-3 Å < d <= 0 Å). The basis for the sub-space definitions is a subdivision into 3 Å layers (about one molecular diameter) parallel to the instantaneous interface, and the numbers in the table are the 10-percentile, the 25-percentile, and the median of energy and order parameter distributions, respectively.    Table S2 for the definition of Bulk AC sub-spaces). For a molecule to be considered inside a sub-space, the criteria for both tetrahedrality parameter and interaction energy must be fulfilled.   Table S3 shows the characteristic time scales for the decay of the time correlation functions.

Time correlation functions of condensing molecules
We note that the short-time decay exhibits a non-monotonous trend with temperature while the long-time decay takes increasingly longer time with decreasing temperature.
Characteristic long-time decay scales of E and q are found to agree well.

Frequency analysis of exchange processes between accommodated and nonaccommodated states
Using the criteria for accommodation defined in the paper, we find that even bulk ice molecules fluctuate between accommodated and non-accommodated states. The analysis of the exchange frequency shows that bulk molecules exhibit a higher frequency than molecules in the surface layer, see Tab. S4 and Fig. S8. We interpret this finding as a result of stronger bonds at crystal sites than at the surface, in analogy to the frequency of a harmonic oscillator depending on the force constant. This picture also explains the trend with temperature since indeed the mean square displacement of molecules about crystal sites due to thermal motion is directly proportional to temperature, assuming a harmonic potential. Due to time binning during analysis, the maximum frequency detectable here is fex,max = 100 ns -1 .  The approximation of harmonic potentials becomes worse close to the surface and at higher temperatures when anharmonic effects have to be expected. Comparing rough estimates of frequencies from the first passage time analysis with the frequency of exchange processes we conclude that they are of the same order of magnitude with the general trend of increasing frequencies with increasing temperature. At the highest temperature of T = 250 K, exchange frequencies of bulk molecules decrease compare to T = 230 K. This is speculated to be due to surface effects gradually reaching down to what in this analysis was defined as "bulk" with increasing temperature. S11 Figure S8: Exchange frequency between AC and NAC states depending on the distance to the instantaneous interface.