Glycine in aerosol water droplets : a critical assessment of Köhler theory by predicting surface tension from molecular dynamics simulations

Introduction Conclusions References


Introduction
Clouds, particularly marine stratocumulus in mid to high latitudes, have been shown to exert a large net negative cooling effect on the Earth (Ramanathan et al., 1989). The radiative or reflective (albedo) properties of clouds are strongly dependent on the number concentration of airborne watersoluble particles, known as cloud condensation nuclei (CCN) (Twomey, 1974). However, this so-called indirect effect Correspondence to: H.Ågren (agren@theochem.kth.se) of aerosol particles on clouds still constitutes a large uncertainty in the models (Solomon et al., 2007) projecting future climate.
Given that the number of CCN over remote marine areas generally shows remarkable stability, the existence of a source that continually restores the CCN population is implied. Blanchard (1971) and Blanchard and Syzdek (1988) have long advocated that a significant proportion of the remote oceanic aerosol is derived from bubble bursting. The bubbles result from entrainment of air induced by wind stress at the air-water interface, which produces primary aerosol particles in CCN sizes. In this process, bubbles scavenge sea-salt, debris and high molecular weight soluble organic surface-active compounds as they rise through the water prior to their injection into the atmosphere.
The discovery of proteinaceous material in Antarctic cloud water sample (Saxena et al., 1983) led us to the assumption that proteins and bacterial enzymes are present in the bursting bubbles. Support for this assumption comes from the presence of free and combined amino acids in marine aerosol and precipitation, which are greatly enriched relative to sea water concentrations because, in the atmosphere, proteins are subject to degradation to amino acids by bacterial enzymes. Typical levels of dissolved amino acids found in marine rain (Mopper and Zika, 1987) showed enrichment by two orders of magnitude over typical seawater values. The amino acids were found predominantly in their L-optical isomers and were therefore probably of biological origin.
The uncertainty in explaining the maintenance of the number of CCN over remote marine areas apparently is a critical element in improving our ability to assess the potential role Published by Copernicus Publications on behalf of the European Geosciences Union. of marine biogenic CCN in climate forcing. Relevant studies include chemical composition, concentration, shape and size of the cloud activating nuclei but also particular properties like surface tension, density, and hydrophilicity of the nuclei.
A well-established model that guides the so-called critical supersaturation, under which the droplet is activated, is given by the Köhler theory (Köhler, 1936). This theory accounts for surface tension and droplet curvature and also the change in water activity and vapor pressure in the presence of solute substances.
The Köhler theory has been widely used to describe the nucleation and the equilibrium growth of cloud droplets since it was first introduced in 1936. Two balancing terms, the Kelvin term and the Raoult term, compose the theory and determine the growth behavior of the cloud droplets. The surface tension of the droplet appears in the Kelvin term and has an important impact on the critical supersaturation. It has been suggested that a decrease in cloud droplet size, which may increase the droplet population and cloud albedo, can result from a decrease of the surface tension of the droplets (Facchini, 1999). The size distribution of the cloud droplets has significant effect on the cloud albedo or reflectivity, and consequently on the absorption of terrestrial infrared radiation and reflection of solar irradiation (Sun and Ariya, 2006). Therefore, the surface tension of aerosol droplets is of great importance and the knowledge of the surface tension could help in understanding climatic issues.
Theoretical prediction of liquid surface tension has been common due to the fast development of the molecular dynamics (MD) and the Monte Carlo (MC) simulation techniques. However, most studies focus on the surface tension of planar interfaces, where the surface tension could be obtained from the diagonal components of the pressure tensor and the length of the simulation box (e.g., Alejandre et al., 1995;Chen and Smith, 2007). Only a few MD or MC simulation studies on liquid drops have been reported (Thompson et al., 1984;Zakharov et al., 1997;van Giessen and Blokhuis, 2009;Sampayo et al., 2010). Thompson et al. studied simple Lennard-Jones liquid drops by using the thermodynamic route and the mechanical route (Thompson et al., 1984); the thermodynamic route requires the knowledge of the surface tension of planar surface and assumes the first-order truncation of the Tolman equation (Tolman, 1949), while the mechanical route is based on the Irving-Kirkwood pressure tensor (Irving and Kirkwood, 1950). Thompson's results suggest a positive Tolman length, which means that the surface tension increases monotonically as the droplet radius increases. Differently, a recent MD study by van Giessen and Blokhuis suggests a negative Tolman length (van Giessen and Blokhuis, 2009). Sampayo et al. also claimed that the second-order contribution in liquid droplets cannot be neglected and the surface tension-radius curve has a weak maximum which indicates a negative Tolman length (Sampayo et al., 2010). Very recently Julin et al. (2010) and Block et al. (2010) studied the curvature-dependent surface tension of simple Lennard-Jones liquid and obtained negative Tolman length by extrapolating to planar limit.
In our previous paper (Li et al., 2010), we reported a molecular dynamics study of the surface-active cis-pinonic acid in small water clusters under atmospheric conditions, and showed that the amount of decrease in surface tension of a droplet is not only related to the concentration of the surfactant, but also dependent on the size of the droplet. In particular, we found that cis-pinonic acid molecules are remarkably concentrated on the droplet surface and, with the increase of the size of the droplet, a lower concentration of cis-pinonic acid is needed to reduce the surface tension by a certain amount. In this paper, we discuss the curvature dependence of the surface tension of pure water droplets and water droplets containing glycine -an amino acid that appears in aerosol particles and behaves as an effective CCN material (Kristensson et al., 2010). Glycine in water is selected as a test system, motivated by that amino acid-water clusters are produced in marine environments by sea spray of phytoplankton and bacteria, and commonly assumed to form aerosols. The choice is also motivated by that activation properties of amino acid-water clusters have been analyzed in laboratory tests (Kristensson et al., 2010). In this work, the surface tension is studied over a range of cluster sizes, from sub-nanometers to the limiting case of a planar interface. By carefully analyzing the dependence of the surface tension on the parameterization of molecular dynamics we find it possible to interpolate the size dependence of the surface tension over the so-called Aitken mode particles, critical for water droplet formation. In our previous study (Li et al., 2010) the surface tension of water droplets has been calculated by using the rigid SPC/E (extended simple point-charge) water model (Berendsen et al., 1987), and the computed values are between 50 and 60 mJ m −2 , lower than the experimental value of 72 mJ m −2 . A more flexible water model which provides more degrees of freedom is expected to give larger surface tension and therefore better agreement with experimental measurements (Alejandre et al., 1995;Zakharov et al., 1997). In this work, we employed the so-called semi-flexible SPC/E water model with constrained OH bond lengths and flexible HOH angle and followed the procedure described in our previous paper (Li et al., 2010) to study the surface tension for water droplets. The radial distributions of molecules in the droplets are analyzed and the curvature dependence of the surface tension was discussed. The surface tension change induced by the presence of glycine is also discussed in the framework of Köhler theory.

Surface tension for planar interface
The surface tension of planar liquid-gas interface is usually divided into two contributions, the original part σ o and the dispersion correction σ d . The original part can easily be calculated from the diagonal components of the pressure tensor as described in the previous paper by Alejandre et al., 1995.
where P xx , P yy and P zz denote the diagonal elements of the pressure tensor and L z is the box length along z-axis.
Due to the use of the truncated Lennard-Jones potential, the dispersion correction to the surface tension has to be evaluated. To do this, the number density of the liquid along the z-axis is fitted to the hyperbolic tangent function where ρ α and ρ β are densities for the liquid phase and the gas phase, respectively. Then the dispersion correction is calculated as (Blokhuis et al., 1995) where we have used the substitution r' = r −1 to perform the double integration on finite intervals. For liquids consisting of multiple components, the average value of εσ 6 is used

Surface tension evaluation of the spherical interface
In our previous paper (Li et al., 2010) we have demonstrated that the surface tension of an isolated spherical water cluster could be evaluated from the work of formation and the radius of the equimolar dividing surface, where the work of formation can be obtained from the Irving-Kirkwood pressure tensor, according to McGraw and Laaksonen (1996) and Zakharov et al. (1997). In this paper we implement the same procedure to calculate the surface tension for spherical interfaces. First, the radial number density is fitted to the hyperbolic tangent function where ρ α and ρ β are densities for the liquid phase and the gas phase, respectively. Then the radius of the equimolar dividing surface is calculated from the density profile The normal component of the Irving-Kirkwood pressure tensor is calculated according to the formula given by Thompson et al. (1984) where k B is the Boltzmann constant, T is temperature, S is the area of the spherical surface of radius r, and f k is the normal component of the force between a pair of molecules acting across the surface S. Afterward the work of formation is obtained by performing integration over the P N (r) function as described by McGraw and Laaksonen (1996) and Zakharov et al. (1997) where R β is the radius of a sphere in the vapor region and P β = P N (R β ) is the vapor pressure. Finally, the approximate surface tension can be obtained from the work of formation and radius of equimolar surface The above equation uses the so-called effective surface tension to approximate the real surface tension, since it is difficult to determine the difference in the pressure of the two phases ( p) or the radius of the surface of tension (R s ) from molecular dynamics simulations. As described in the papers by Brodskaya et al. (1996) and Zakharov et al. (1997), the effective surface tension is equal to (γ e − R e dγ e /dR e ), where the dividing surface is chosen to be the equimolar surface, the radius of which can be readily determined from molecular dynamics simulations. The effective surface tension serves as an excellent approximation to γ e and has the same asymptotic behavior as γ e at large R e , especially at the sizes of the cloud condensation nuclei (CCN), which are of interest in atmospheric research. Our idea is to use the effective surface tension to describe the real surface tension semiquantitatively, and our previous paper (Li et al., 2010) has demonstrated that the effective surface tension can properly reproduce the decreased surface tension in droplets containing cis-pinonic acid.

Computational details
The GROMACS program package (version 4.0.7) (Hess et al., 2008;van der Spoel et al., 2005) was used to carry out the molecular dynamics simulations. For water molecules the semi-flexible SPC/E (extended simple point-charge) water model (Berendsen et al., 1987) was employed, with the two OH bond lengths constrained at 1.0Å by the LINCS (linear constraint solver) algorithm (Hess et al., 1997;Hess, 2008) and the HOH angle kept flexible. The HOH angle was modeled by a harmonic potential with an equilibrium angle of 109.47 • and a force constant of 383.0 kJ mol −1 rad −2 . For zwitterionic glycine molecules the OPLS (optimized potentials for liquid simulations) all-atom force field (Jorgensen et al., 1996) was used, with the bonds containing H atoms constrained by the LINCS algorithm (Hess et al., 1997;Hess, 2008). All simulations were performed in canonical (constant-NVT) ensembles, with temperature maintained at 298 K by the Nosé-Hoover thermostat (Nosé, 1984;Hoover, 1985). Periodic boundary conditions were applied during the simulations. The van der Waals interaction was modeled by the Lennard-Jones potential truncated at the cutoff radius of 10Å. The Coulomb interaction was evaluated by the Ewald summation, where the real part of the potential is truncated at 10Å and the reciprocal interaction is calculated by the SPME (smoothed particle mesh Ewald) (Darden et al., 1993;Essmann et al., 1995) method with a mesh of reciprocal vectors of 1.2Å in every direction and a spline of order 4. The time step of the simulations was set to 2 fs.
To model the surface tension of spherical and planar interfaces, twelve systems consisting of either pure water or glycine/water mixtures were constructed, as listed in Table 1. Each system was initially set up from pre-equilibrated cubic simulation boxes, by using the "genbox" utility of the GROMACS program package (van der Spoel et al., 2005). In systems 1-10, the pre-equilibrated simulation boxes were enlarged in three dimensions so that the lengths of the boxes were 6.0 nm longer than the diameters of the spherical clusters, in order to minimize the interaction between the cluster and its periodic images. Here the diameter of the spherical cluster was estimated from the volume of the initial cubic box. In systems 11-12, the boxes were elongated three times along the z-axis to create the planar liquid-gas interfaces. Each system was subject to energy minimization and then 2 ns simulation to reach equilibrium, and another 2 ns simulation was carried out and the trajectories are saved every 0.2 ps for further analysis and calculations.
After the simulations were finished and the trajectories recorded, surface tension calculations were carried out as described in Sect. 2. The surface tension values for planar interfaces were evaluated from the diagonal components of the pressure tensor, while for spherical interfaces the surface tension data were obtained from the normal component of the Irving-Kirkwood pressure tensor. The curvature dependent surface tension was modeled by a polynomial function Table 1. Number of molecules and length of the box in the twelve systems. In systems 1-10 the interfaces are spherical and L z = L x = L y , while in systems 11-12 the interfaces are planar and L z = 3L x = 3L y . of the inverse of the radius of the equimolar dividing surface, and the glycine-induced surface tension change was extracted from the calculated data. The densities and radial distributions were also analyzed. To get insight into the effect of glycine on the growth of aerosol droplets, the Köhler theory was used to model the activation behavior of the droplets. The Köhler theory (Köhler, 1936) describes the nucleation and the equilibrium growth of cloud droplets. Two terms comprise the theory; the Kelvin term relates to the curvature of the droplet, which decreases the ability of droplet nucleation, a necessary step prior to further growth of the droplets, while the Raoult term increases the ability for droplet growth because of the solution of the aerosol constituents inside the droplet. Here the Kelvin term is related to the surface tension of the droplet, and an increase in the surface tension would result in an increase of the critical supersaturation.

Spherical interface simulations
The calculated results for systems 1-6 are listed in Table 2, including the liquid density, the radius for equimolar dividing surface, the work of formation, and the surface tension. A decrease in density is found as the size of water droplet grows, and all the densities are larger than the standard density of the SPC/E water model (0.998 g cm −3 (Berendsen et al., 1987), equivalent to 33.361 nm −3 ). This indicates a compression in water droplet induced by the surface tension, although such effect is very small. As the droplet grows larger, the effect of surface tension on the density of bulk phase becomes smaller. The radial number densities of water molecules are shown in Fig. 1, which can readily be fitted to the hyperbolic tangent Atmos. Chem. Phys., 11, 519-527, 2011 www.atmos-chem-phys.net/11/519/2011/ function to derive the radius for the equimolar dividing surface as described in Sect. 2. It is worth noting that the density profile can also be fitted to the error function, and the results are shown in parentheses in Tables 2, 3 and 4. It is found that the use of the error function leads to negligible changes in the fitted densities as well as in R e and σ . Therefore we rely on the data obtained by fitting to the hyperbolic tangent function.
The calculated normal components of the Irving-Kirkwood pressure tensors in systems 1-6 are shown in Fig. 2. In previous simulations for pure Lennard-Jones droplets (Thompson et al., 1984), it has been hypothesized that the P N (r) profile could also be fitted to a hyperbolic tangent function like density profile, however, this is not the case for water droplets we studied here. As can be seen from Fig. 2, each P N (r) curve has a peak near the surface, which could be tentatively attributed to the presence of the strong surface tension. In the bulk phase of each water droplet, the P N curve is flatter, and the P N value in bulk phase decreases as the size of the droplet increases. This again shows that the effect of surface tension on the bulk phase becomes smaller as the droplet grows larger.
Listed in Table 3 are the computational results for systems 7-10. Clearly the densities in these systems are smaller than those in systems 1-6, due to the presence of the larger glycine molecules. The density decreases from system 7 to system 9, while in system 10 the density increases. A plausible interpretation to this phenomenon is that the glycine molecules are more uniformly distributed in larger droplets and that the decrease in the density would not be as pronounced as in Table 4. Calculated liquid density, fitting parameter ξ , original part and dispersion correction of surface tension for planar interfaces in systems 11-12. Shown in parentheses are the densities obtained by fitting to the error function. smaller droplets. Figure 3 shows the number densities for glycine molecules and water molecules in the droplets. The peaks in the 0 < r <0.5 region do not mean much due to the lack of statistics. It can be seen that in all droplets the density of glycine vanishes at the surface, indicating that the zwitterionic glycine molecules are very hydrophilic and prefer the bulk phase of the droplet. The P N (r) profiles for systems 7-10 are also calculated, as shown in Fig. 4. Comparing to Fig. 2, the P N (r) curves for systems 7-10 also have similar peaks near the surface. This is in accordance with the fact that the glycine molecules stay inside the droplet and have minor effect on the surface of the droplet. Consequently, the P N (r) curve in bulk phase of the droplet is significantly affected by the presence of glycine molecules. The glycine-induced P N (r) is also shown in Fig. 4. Here P N (r) is the difference between the P N (r) values of the glycine-containing droplet and the pure water droplet. It can be seen that the presence of glycine molecules has an undulate effect on the P N (r) curve. Although glycine molecules mainly reside in the bulk of the droplet, a small influence can be found in the surface region of the P N (r) curve. Inside the droplet where glycine molecules assemble, the value of P N (r) goes negative and therefore tends to reduce the surface tension. As the droplet size increases, the P N (r) curve becomes closer to zero inside the droplet and  preserves the positive peak at the surface. This explains the trend of the change of the glycine-induced surface tension in droplets with different sizes.

Planar interface simulations
The calculated densities and surface tension values for planar interfaces are listed in Table 4. The density in system 11 is smaller than those in systems 1-6, and the density in system 12 is smaller than those in systems 7-10, proving the compressing effect of the surface tension on spherical interfaces. The calculated surface tension for pure water (system 11) is obtained as σ o + σ d = 65.63 mJ m −2 . This value is smaller than that reported by Alejandre et al. (1995) since they used limited number of lattice vectors in the reciprocal sum when evaluating the long-range electrostatic interaction, which is expected to overestimate the surface tension (Chen and Smith, 2007). On the other hand, our calculated surface tension for the planar water-air interface is larger than that reported by Chen and Smith (2007) in that we used the semi-flexible SPC/E model which could give larger surface tension and better agreement with experimental value than the rigid model. In the presence of glycine molecules (system 12) the surface tension becomes 66.17 mJ m −2 . Since the concentration of 10 glycine molecules in 1000 water molecules is 0.56 mol L −1 , the slope of surface tension with respect to concentration can be estimated as dσ /dC=0.96 mJ m −2 L mol −1 , in good agreement with experimental measurements: dσ /dC=0.92 (Pappenheimer et al., 1936) or 1.12 (Bull and Breese, 1974).

Curvature dependence of the surface tension
Regarding the planar interface as spherical interface with infinite radius, the relationship between the surface tension and the inverse of droplet radius could be drawn by fitting to a polynomial function. Tolman suggested a linear relationship between σ and 1/R by neglecting higher-order terms (Tolman, 1949), however, our results suggest a quadratic relationship, with non-negligible contribution from the 1/R 2 term. The two fitted curves in Fig. 5 are σ W = 65.45 + 50.66/R e -100.00/ R 2 e and σ GLY = 66.27 + 56.92/R e -132.71/R 2 e , respectively, where R e is the radius of the droplet (in nm) and σ is the surface tension (in mJ m −2 ). Note that the second equation is for a glycine concentration of 0.56 mol L −1 . As shown in Fig. 5, the two curves intersect at R e = 3.57 nm, meaning that glycine will suppress the surface tension for droplets with smaller radius while increase the surface tension for droplets with larger radius. Assuming the deviation of the surface tension is linearly dependent on the concentration of glycine, as observed in experimental measurements (Pappenheimer et al., 1936;Bull et al., 1974) where C is the concentration of glycine (in mol L −1 ). These formulas are suitable for calculating the surface tension of aerosol droplets containing low-concentration of glycine.

Surface tension correction to the Köhler theory
The Köhler theory is formulated in terms of the Raoult effect and the Kelvin effect where p is the vapor pressure of water over a droplet, p 0 is the water vapor pressure over a flat surface of pure water, a w is the water activity of the solution, M w is the molar mass of water, σ is the surface tension, R is the gas constant, T is the temperature, ρ is the density of the solution, and D p is the droplet diameter. By applying the surface tension formula to the Köhler equation, the Köhler curve for a glycine dry particle with diameter of 50 nm could be drawn, as shown in Fig. 6. All parameters in the Köhler equation were taken from the paper by Kristensson et al. (2010). The black line denotes the Köhler curve obtained by using the surface tension of pure water, that is, σ = σ W = 72.4. The red dashed line denotes the Köhler curve obtained by using the planar surface tension correction: σ = σ W + (66.27-65.45)· C/0.56, σ W = 72.4. The blue line denotes the Köhler curve obtained by using the curvature-dependent surface tension correction, σ = σ W + (0.82 + 6.26/R e -32.71/R 2 e )· n √ C/0.56, σ W = 72.4 + 50.66/R e -100.00/R 2 e . Here the value of σ W has been shifted upward so that σ W is equal to the experimental value of 72.4 when R e is infinity. All the surface tension values are in mJ m −2 . From Fig. 6 it can be seen that for the three cases the critical supersaturation (the supersaturation where the Köhler curve is at its maximum) are 0.546%, 0.548%, and 0.553%, respectively, while the experimental value reported by Kristensson et al. (2010) is between 0.6% and 0.7%. Clearly the surface tension correction gives rise to a better agreement with experimental observation. Nevertheless, such correction is relatively small, and other factors affecting the critical supersaturation cannot be ruled out. As suggested by Kristensson et al. (2010), the discrepancy between the experimental results and the Köhler model arises from a combination of factors including surface tension, non-ideality, and uncertainty in the density of glycine. Further study is necessary to resolve the underestimation of the critical supersaturation by the Köhler model.

Conclusions
In summary, atmospheric aerosol droplets, representative for remote marine areas, consisting of either pure water or glycine-water mixtures, were studied by molecular dynamics simulations. The surface tension for each droplet was estimated from the normal component of the Irving-Kirkwood pressure tensor. Different from the simple Lennard-Jones liquid droplet, the P N (r) profile for water cluster possesses a peak at the surface, indicating that the surface molecules are under strong tension. This is in accordance with the fact that water has a large surface tension and a low vapour pressure under ambient conditions. The radial distributions of glycine molecules in water droplets were also studied, and it was found that the radial densities vanish at the surface. In other words, the zwitterionic glycine molecules are very hydrophilic and tend to stay in the bulk phase of the droplet. The P N (r) profile of the droplet was perturbed by the presence of glycine molecules, resulting in glycineinduced change of surface tension. To obtain the curvature dependence of the surface tension, the planar interface was also investigated as the limiting case. The surface tension for the planar water-air interface was calculated to be 65.63 mJ m −2 , close to the experimental value of 72 mJ m −2 . The surface tension change induced by glycine is computed as dσ /dC = 0.96 mJ m −2 L mol −1 , in agreement with experimental measurements.
By interpolating the surface tension over the nano-sized droplets and the limiting case of the planar interface, a quadratic relationship is found between the surface tension and the inverse of the radius of the droplet, indicating a strong curvature dependence of the droplet surface tension. In this case the Aitken mode particles with larger sizes could be covered and the Köhler equation could be improved by incorporating surface tension corrections. By applying our results on the activation process of glycine dry particle with diameter of 50 nm, the critical supersaturation is improved from 0.546% to 0.553%. The correction is not so pronounced, since glycine is not a strong surfactant as, for instance, cispinonic acid and its effect on the surface tension is expected to be small. Nevertheless, our approach provides an efficient way of evaluating the surface tension for aerosol droplets, which could be of use in improving the accuracy of the Köhler equation in predicting the critical supersaturation of the activation process of cloud droplets in general and with specific interest to the remote marine aerosol derived from biology at the ocean surface.