A new modeling tool for the diffusion of gases in ice or amorphous binary mixture in the polar stratosphere and the upper troposphere

To elaborate stratospheric ozone depletion processes, measurements of diffusion coefficients of selected gas phase molecules (i.e. HCl, CH 3OH, HCOOH and CH3COOH; Katsambas et al., 1997; Kondratyev and Varotsos, 1996; Varotsos et al., 1994, 1995) in ice in the temperature range 170–195 K have been analyzed with respect to the mechanisms and rates of diffusion. It is argued that the diffusion in ice of these compounds is governed by a vacancy – mediated mechanism, i.e. H 2O vacancies are required to diffuse to lattice sites adjacent to these compounds prior to the diffusion of the corresponding molecule into the vacancy sites. In addition, we show that the diffusion coefficients of these compounds exhibit a specific interconnection, i.e. a linear relationship holds between the logarithm of the pre-exponential factor, Do, and the activation energy E. The physical meaning of this interconnection is discussed.


Introduction
It is generally accepted that H 2 O ice particles play a prominent role in the heterogeneous chemistry in the polar stratosphere that leads to ozone depletion (Solomon, 1988;Tabazadeh and Turco, 1993;Kondratyev et al., 1994;Kondratyev andVarotsos, 1995, 2001a, b;Reid et al., 1998;Schulz et al., 2000Schulz et al., , 2001;;Varotsos, 2002).Moreover, such particles, apart from their involvement to various mechanisms in the lower troposphere (Asimakopoulos et al., 1992;Cartalis and Varotsos, 1994;Ferm et al., 2005Ferm et al., , 2006) also contribute to radiative forcing as well as to chemical conversions in the upper troposphere where they exist in form of cirrus clouds and contrails.Since the rates of heterogeneous chemical conversions depend on both surface and bulk Correspondence to: C. A. Varotsos (covar@phys.uoa.gr)concentrations, diffusion measurements in or on ice have attracted major interest.
When in a solid a single diffusion mechanism is operative, the diffusion coefficient, D, is often found to obey an Arrhenius type behaviour, i.e.
where the activation energy E and the pre-exponential factor D o are essentially temperature independent and k B is the Boltzmann constant.In this case, the slope of the linear lnD vs. 1/T plot leads to the value of E which -as shown by a careful thermodynamical study (Varotsos andAlexopoulos, 1979, 1986) -coincides with the enthalpy of activation H act of the diffusion process.Whilst such findings apply to the rate of diffusion they do not provide any information on the diffusion mechanism.
Early self-diffusion studies in ice have been performed using microtome sectioning and scintillation or mass spectrometric probes (H 18  2 O, D 2 O, T 2 O) near the ice melting point (see Petrenko and Whitworth, 1999 and references therein).These studies led to a diversity of views concerning the mechanism of H 2 O isotopic diffusion in ice (Petrenko and Whitworth, 1999).Opinions have been expressed, for instance, that for T >223 K self-diffusion occurs by an interstitial-mediated mechanism, while for T <223 K a vacancy-mediated process is operative (Goto et al., 1986;Hondoh et al., 1989).
During the last decade, new infrared-laser resonant desorption (LRD) techniques have been employed to derive diffusion coefficients by depth profiling of ice films (Livingston andGeorge, 1999, 2001;Livingston et al., 1997Livingston et al., , 1998Livingston et al., , 2000Livingston et al., , 2002;;Krasnopoler and George, 1998).These measurements revealed two different rate categories for bulk diffusion in ice.Whilst compounds such as HCl, CH 3 OH, HCOOH and CH 3 COOH display diffusion kinetics similar to the kinetics of H 2 O self-diffusion, species such as Na and NH 3 exhibit diffusion kinetics appreciably slower than for H 2 O self  Fig. 1.Arrhenius plots for temperature dependent diffusion coefficients for various species in ice as obtained by the LRD depthprofiling technique (Livingston et al., 2001(Livingston et al., , 2002) ) and chemical titration (Nehme, 2006).The figures on the individual lines are activation energies (in kcal/mol) as derived by the authors.
diffusion.Livingston et al. (2002) argued that the diffusion process of the compounds of the first category occurs via a vacancy-mediated mechanism in which H 2 O vacancies may be required to diffuse to lattice sites adjacent to the impurity compounds prior to the diffusion of this compound into the vacancy sites.The present paper focuses on the heterodiffusion data of compounds of the first category.We will draw special attention to the fact that the diffusion data of these species are interconnected in a specific manner.

Theoretical considerations of hetero-diffusion in ice
In Fig. 1 we present the temperature dependence of kinetic data for bulk diffusion in ice for selected gas phase compounds as measured by infrared LRD depth profiling (Livingston et al., 2002) as well as by chemical titration (for CH 3 COOH only; Nehme, 2006)   E have been reported.The corresponding Arrhenius parameters are summarized in Table 1.All these parameters-except those for CH 3 COOH in the higher temperature range-are also summarized in Table 1 of Livingston et al. (2002), where the values for CH 3 OH are reported for the first time, while those for HCl come from Livingston et al. (2001).As for the corresponding values for HCOOH and (the lower temperature data of) CH 3 COOH Livingston et al. (2002) state that they come from earlier measurements by Livingston and co-workers.Table 1 also  According to the theory of diffusion in solids, the temperature dependence of the diffusion coefficient D for a single diffusion mechanism in a given crystal is given by (Varotsos and Alexopoulos, 1986): where ( G act ) is the activation Gibbs energy for the diffusion process, (f ) a numerical constant which depends on the mechanism, (α) the lattice constant and (ν) the attempt frequency.To a first approximation the attempt frequency depends on the mass of the diffusing species according to (Varotsos and Alexopoulos, 1986): where m m , m i denote the mass of the matrix (m) and the diffusing species (j ), respectively, and ν m , ν j for the corresponding attempt frequencies.
The Gibbs energy G act can be decomposed into an entropy S act = −d G act /dT | p ) and an enthalpy H act = G act +T d G act /dT | p ) term whereupon Eq. ( 2) turns into: which can be alternatively written as: with As a consequence the temperature dependence of the diffusion coefficient follows a simple exponential behaviour, for which the pre-exponential factor D o is associated with a dynamic (frequency) term as well as a thermodynamic (entropy) term.For the same matrix therefore, low values of D o hence reflect low attempt frequencies or low entropies of activation.It should be noted that S act is always positive because the diffusing species is more loosely bound or less ordered at the transition state than at the individual lattice sites.
where (B) is the isothermal compression modulus, ( ) the mean volume per atom and (c act,j ) a constant independent of temperature and pressure.The additional subscript j clarifies here that the value of this constant may be different for various diffusing species j .
The validity of this model has been tested for various defect processes and different categories of solids and the results were compiled by Varotsos and Alexopoulos (1986).
Moreover, it was later found (Varotsos et al., 1999) that this model also holds even for more complex defect processes that e.g. are responsible for the emission of electric pulses observed upon pressurizing a solid, which may explain the electric signals detected before major earthquakes (Varotsos and Alexopoulos, 1984b;Varotsos et al., 1996Varotsos et al., , 2003Varotsos et al., , 2005aVarotsos et al., , b, 2006a, b), b).On the basis of Eq. ( 7), the corresponding entropy ( S act,j ) and enthalpy ( H act,j ) terms are found to be given by S act,j = -c act,j (βB where (β) is the thermal volume expansion coefficient.Thus, by inserting Eq. ( 8) into Eq.( 6), the pre-exponential (D o ) value for the diffusion of species j (labeled hereafter D oj ) is given by (when Eq. ( 3) is also considered): Together with Eq. ( 9), this relation reveals that for the diffusion of various species j into the same matrix material (m), when plotting ln{D oj /[f α 2 ν m (m m /m j ) 1/2 ]} versus H act,j a straight line should result, the slope of which is solely governed by the elasticity and expansivity properties of the host material.Since in the present case we only consider water ice as a matrix and data for other hosts are not easily available, the variation of diffusion data with these properties has yet not been tested.
In the present case the D oj values vary by several orders of the magnitude, i.e., the D o value for HCOOH exceeds that of HCl by as much as seven orders of magnitude (see Table 1; the error for D, e.g., for the case of HCl, is around 30% as given in Table 2 of Huthwelker et al., 2006).In view of this large variability, the mass effect associated with the frequency term ν m (m m /m j ) 1/2 which emerged from Eq. ( 3) can, to a first approximation, be discarded.This explains why the plot of the lnD oj versus H act,j or E is almost a straight line.
The apparent success of this analysis is, except for other potential sources, mainly due to the invariance of the diffusion mechanism.In all cases considered here this mechanism is identical, i.e. diffusion is governed by a vacancy controlled mechanism.This in turn leads to the conclusion that diffusion data that do not fall onto the line delineated here correspond to different diffusion mechanisms.This, for instance, applies to NH 3 and Na.
The following comments, however, are worthwhile to be added as far as the diffusion data analyzed here.It has been argued (Domine and Xueref, 2001) that the diffusion coefficients reported by Livingston et al. (2002) were too large to be interpreted as molecular diffusion in crystalline ice.This is so because in most solids, diffusion coefficients are of the order of 10 −10 to 10 −12 cm 2 /s near the melting point while Livingston et al. (2002) obtained values in that range 100 K below the melting point.Domine and Xueref ( 2001) have shown convincingly using IR spectroscopy that the ice samples used by Livingston et al. (2002) were not real ice, but were in fact amorphous binary mixture with a degree of disorder such that diffusion could be expected to be fast.These arguments by Domine and Xueref (2001) seem to be well supported by the following fact: Livingston et al. (2002) obtained diffusion coefficient values near 180K higher than those reported earlier by other authors in the range 240-260 K (cf.since diffusion is thermally activated, the diffusion coefficients should decrease upon decreasing the temperature).It seems therefore reasonable to accept that the ice samples used by Livingston et al. (2002) were in fact amorphous binary mixture rather than real ice.The interesting fact that even in this case the present analysis seems to be well obeyed (since the plot lnD o versus H act , in fig 2 is almost a straight line) stems from the following: Eq. ( 11) is based on the cB model, i.e., Eq. ( 7), which is of thermodynamic origin and has been shown (Varotsos and Alexopoulos, 1986) to be valid for amorphous binary mixtures as well (cf.the B value of such a mixture can be estimated in terms of B−values of the end constituents, see Varotsos, 1980).

Discussion
We stress that measurements of diffusion in ice are nowadays difficult and many measurements have been debated in the past.The debate holds even when restricting ourselves to the diffusion data of a given gas molecule and tackles several issues like, polycrystallinity, amorphous ice, hydrates etc.As an example, let us focus hereafter on the data for the HCl diffusion.The data available before 2006 have been reviewed by Huthwelker et al. (2006) and have shown a large scatter (as illustrated in their Fig.20).The data vary by more than a factor of 10 9 .After a careful inspection and study of this scatter, Huthwelker et al. (2006) concluded that when investigating the diffusion of HCl in ice, the thermodynamics of the system must be considered.In particular, the following cases should be distinguished: First, HCl partial pressures high enough to melt the ice will yield extremely high diffusion constants.Second, hydrate formation can strongly affect the transport of HCl through ice and leads to intermediate values for the diffusion constant.Third, measurements in the ice stability region yield the lowest values for the diffusion coefficient.At 200 K, the available data for these conditions refer only to the near-surface region, i.e., the upper microm-eter range of the ice surface.Interestingly, it was guessed (Huthwelker et al., 2006), that these diffusion constants may represent the diffusivity into small micrometer-sized atmospheric ice particles.These values, however, may not represent the bulk diffusivity of HCl in ice (because hypothesized surface melting effects may affect the transport speed of HCl in the region close to the ice surface).This cannot be confirmed, because no direct measurements of the HCl diffusion in the bulk of single crystals are available for temperatures lower than 243 K.
We now proceed to data published after 2006.Using fast thermal desorption spectroscopy, a novel technique recently developed, Lu et al. (2007) investigated the kinetics of H/D isotopic exchange in 3 µm thick polycrystalline H 2 O ice films containing D 2 O layers at a temperature 271±1.5 K.These data, in combination with a theoretical analysis of the diffusion in polycrystalline medium, gave an estimate of the range of possible values for water self-diffusion coefficients and the grain boundary widths characteristic of the ice samples used.Their analysis showed that the diffusivity of D 2 O along the grain boundaries must be at least two orders of magnitude lower than that in bulk water.Based on these results, Lu et al. (2007) argued that, in the limit of low concentration of impurities, polycrystalline ice does not undergo grain boundary premelting at temperatures up to 271 K.More recently, Lu et al. (2009) extended such fast thermal desorption spectroscopy studies, in the temperature range 255 to 272 K, to a few micrometer thick, pure polycrystalline ice film as well as to ice films doped with HCl.For low impurity concentration, the effective water self-diffusion coefficients were found to exhibit an Arrhenius dependence with activation energy almost equal to that published by Livingston et al. (2002).This activation energy for water self-diffusion along grain boundaries inferred from H/D exchange data indicates that the diffusion mechanism at these interfaces is similar to that in single crystal ice.Thus, according to Lu et al. (2009), grain boundaries in a pure polycrystalline ice sample can be viewed simply as bulk ice characterized by a high density of interstitial and vacancy defects.The water selfdiffusion coefficients in single crystal ice are approximately two orders of magnitude lower than those in their pure polycrystalline ice samples.The potential effects of polycrystallinity, including the case that the character of ice samples (single-vs.poly-crystals, degree of crystallinity and orientational homogeneity) may play a role in species mobility, has been earlier discussed by Rempel and Wettlaufer (2003).Upon adding 0.04% (by mass) of HCl to polycrystalline ice samples, the effective water self-diffusion coefficient measured by Lu et al. (2009), increases by one to two orders of magnitude.For temperatures lower than 265 K, the activation energy in the Arrhenius plot is almost equal to that found by Livingston et al. (2002), but it nearly doubles above 265 K.In other words, above 265 K an upward curvature of the Arrhenius plot is found in the case of HCl doped ice, which is likely to be caused by onset of grain boundaries premelting.
In general, an upward curvature of the Arrhenius plot is usually attributed to the coexistence of at least two operating diffusion mechanisms.In such a case, however, we do not expect to find a behavior similar to that observed in Fig. 1.This is so, because Eq. ( 11), which predicts a linear relationship between the logarithm of D oj and the activation enthalpy H act,j , has been derived under the condition of a single operating mechanism (cf.For the case of two mechanisms, for example, Eq. ( 2) does not hold since it must be replaced by another one having a sum of two exponential terms).This condition seems to be fulfilled for the diffusion process of HCl, CH 3 OH, HCOOH and CH 3 COOH depicted in Fig. 1, because as mentioned in Sect. 1, Livingston et al. (2002) argued that this process occurs solely via a vacancy-mediated mechanism.A second prerequisite for the linearity when plotting the logarithm of D oj vs. H act,j is, as mentioned, that the diffusion of various species i should take place into the same matrix material.Only in this case, the elasticity and expansivity properties (i.e, the quantities β, B and dB dT ) in the right side of Eq. ( 11) remain constant since they refer to the same host material.As a result, since single crystals of ice have elasticity and expansivity properties different from those of polycrystalline ice, we expect that diffusion data from these two type of crystals (such as those compiled in Table 2 and Fig. 20 of Huthwelker et al, 2006) should not lead to a linearity when plotting the logarithm of D oj vs.
H act,j In other words, in order to observe such a linearity attention should be paid to select carefully diffusion data that refer to identically prepared laboratory ices.This is important in view of the following (see Marchand et al., 2006, and references there in): The acute and complex dependence of the phase, morphology and microstructure of laboratory ice, on preparation methods (vapour condensation or crystallization from the melt) and growth conditions (flux/pressure, temperature, nature of the heterogeneous substrate etc.) and the great difficulty to experimentally quantify defect densities (vacancies, interstitials, dislocations etc.) are all factors that severely limit meaningful comparisons between the results from these different studies.In combination with our limited ability to accurately characterize the morphology and microstructure of this delicate material, these considerations contribute to the large discrepancies between the results obtained from laboratory ice samples prepared by different methods (recall that the aforementioned data for HCl diffusion as compiled in Fig. 20 of Huthwelker et al. (2006), exhibit wide scatter extending over several orders of magnitude).These are the reasons, why the linearity of the logarithm of D oj vs. H act,j plot as predicted by Eq. ( 11), has been checked here by relying on the data by Livingston and coworkers that have been obtained from identically prepared laboratory samples (irrespective if they were amorphous binary mixture and not real ice, as mentioned in Sect.2).

Conclusions
It is well-recognised that H 2 O ice particles play a key role in the ozone depletion in the polar stratosphere and contribute to radiative forcing (e.g.Reid et al., 1994;Gernandt et al., 1995;Efstathiou et al., 2003;Cracknell and Varotsos, 1994, 1995, 2007;Tzanis and Varotsos, 2008;Alexandris et al., 1999).Temperature dependent measurements of diffusion coefficients for HCl, CH 3 OH, HCOOH and CH 3 COOH in ice, as studied in recent measurements (Livingston et al., 2002;Nehme, 2006) show that their D o -values differ by several orders of magnitude.By contrast, the absolute values of the diffusion coefficients do not differ significantly (see Fig. 1), indicating that the transport of material occurs with more or less comparable rates for the same mechanism.Despite this variation, we show that the (natural) logarithm of D o varies almost linearly with the activation energy E for diffusion which for linear ln D vs. 1/T plots coincides with the activation enthalpy, H act .If values for the pre-exponential factors D o could be estimated (i.e. by considering the entropy change during the diffusion process) then the correlation derived here could be of substantial predictive power for diffusion coefficients of other species in ice.The present analysis holds even if we consider convincing arguments clarifying that the samples used by Livingston et al. (2002) were in fact amorphous binary mixtures and not real ice.
Edited by: W. E. Asher

Fig. 2 .
Fig. 2. Relationship between the pre-exponential factor D o of the diffusion coefficient and the activation energy E.
includes the temperature range over which these data have been obtained as well as the absolute values of D for a temperature of 180 K.As can be seen from this compilation the pre-exponential factors D o are vastly different; small values of D o correspond to smaller activation energies E and vice versa.On the other hand, the differences in the absolute values of D are much less pronounced.For further analysis of these data we plot in Fig.2the (natural) logarithm of D o versus the activation energy E.An inspection of this figure reveals that, for the diffusion of the species selected here, a near linear relation holds.A least squares fit to a straight line leads to a slope 2.502 and an intercept −20.365 (when D o is measured in cm 2 /s and E in kcal/mol) with a correlation coefficient 0.98.The maximum uncertainty in this slope is at the most around 10%, if one considers the experimental errors in the values of E and D o extracted from the diffusion data.

Table 1 .
Summary of rate parameters for the diffusion coefficients of selected compounds in ice.