Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Simulations of water adsorption and structure on kaolinite surfaces Croteau, Timothé 2010

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2010_spring_croteau_timothe.pdf [ 3.81MB ]
Metadata
JSON: 24-1.0060989.json
JSON-LD: 24-1.0060989-ld.json
RDF/XML (Pretty): 24-1.0060989-rdf.xml
RDF/JSON: 24-1.0060989-rdf.json
Turtle: 24-1.0060989-turtle.txt
N-Triples: 24-1.0060989-rdf-ntriples.txt
Original Record: 24-1.0060989-source.json
Full Text
24-1.0060989-fulltext.txt
Citation
24-1.0060989.ris

Full Text

Simulations of Water Adsorption and Structure on Kaolinite Surfaces by Timoth´e Croteau B.Sc., McGill University, 2001 M.Sc., The University of British Columbia, 2004 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Chemistry)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) February 2010 c Timoth´e Croteau 2010  Abstract Grand canonical Monte Carlo calculations are used to determine water adsorption and structure on kaolinite surfaces, with and without the presence of trench-like defects, as a function of relative humidity (RH), at 235 K and 298 K. Both basal planes (the Al- and Si-surfaces), as well as two edge-like, defect free surfaces are considered. The trenches simulated are rectangular in geometry, and have a fixed depth and varying width. The general force field CLAYFF is used together with the SPC/E and TIP5P-E models for water. At both 235 K and 298 K, the edges, Al-surface, and trenches adsorb water at sub-saturation, in the atmospherically relevant pressure range. The Si-surface remains dry up to saturation. Both edges and the Al-surface adsorb water up to monolayer coverage. Adsorption on the Al-surface exhibits properties of a first-order process with evidence of collective behavior, whereas adsorption on the edges is essentially continuous and appears dominated by strong water lattice interactions. Only next to the Al-surface, were hexagonal rings observed in the water layer. However, they did not match hexagonal ice Ih. The results obtained using trenches show that the granularity of the surfaces can play a major role in the adsorption of multiple layers of water over a large range of RH. Our calculations suggest that water adsorption in trenches, and possibly in other similar defects, can offer an explanation of the large water coverages reported experimentally. Related to ice, the very dense, proton ordered, ferroelectric structures found in the trenches at 235 K do not correspond to any recognizable form of bulk ice. We speculate how these structures might aid ice nucleation and growth, and suggest how this possibility could be further explored with simulations and experiments. ii  Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iii  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  vi  List of Tables  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  x  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xi  Dedication  Statement of Co-Authorship . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Mineral Dust . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Importance of Dust . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Role of Mineral Dust in Atmospheric Chemistry . . . . . . 1.2.2 Role of Mineral Dust in Climate: Direct and Indirect Effect 1.3 Kaolinite Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Thesis Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Review of Simulation . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . .  . . . . . . . . .  . . . . . . . . .  . . . . . . . . .  . . . . . . . . .  1 1 2 2 3 4 5 8 9  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2 Adsorption and Structure of Water on Kaolinite 2.1 Introduction . . . . . . . . . . . . . . . . . . . . 2.2 Model and Method . . . . . . . . . . . . . . . . . 2.3 Results and Discussion . . . . . . . . . . . . . . . 2.3.1 Adsorption Isotherms . . . . . . . . . . . 2.3.2 Water Structure and Lattice Match . . . 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . .  Surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  13 13 15 20 20 23 26  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3 Simulation of Water Adsorption on Kaolinite under Atmospheric Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 iii  Table of Contents 3.2 3.3  3.4  Model and Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Water Adsorption at 298 and 235 K . . . . . . . . . . . . . . 3.3.2 Orientations and Binding Energies of Single Water Molecules 3.3.3 Water Structure at Monolayer and Sub-Monolayer Coverage . 3.3.4 Water Density and Hydrogen Bonding at Monolayer Coverage Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  33 39 39 42 48 50 52  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4 Water Adsorption on Kaolinite Surfaces Containing Trenches . . . . 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Model and Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Water Adsorption and Structure at 298 K . . . . . . . . . . . . . 4.3.2 Binding Energies . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.3 Comparison of Calculated Adsorption Isotherms with Experiment 4.4 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . .  . . . . . . . .  58 58 60 66 66 72 76 78  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5 Structural Analysis of Water Adsorbed on Kaolinite Particles Containing Trenches at 235 K . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.2 Model and Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.3.1 Water Adsorption at 235 K . . . . . . . . . . . . . . . . . . . . . . 87 5.3.2 Water Density Profiles . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.3.3 Structural Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.3.4 Orientational Ordering . . . . . . . . . . . . . . . . . . . . . . . . 96 5.3.5 Hydrogen Bonding . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.4 Relationship between Observations and Ice Nucleation in the Atmosphere 105 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6 Conclusions and Perspective . . . . . . . . . . . . . . . . . . . . . . . . . . 110 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113  I  Appendices  114  A Simulation Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 A.1 Molecular Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 A.1.1 Long-Range Forces . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 iv  Table of Contents B Simulation Methodology . . . . . . . . B.1 Monte Carlo Simulation Method . . B.1.1 Statistical Mechanics . . . . B.1.2 Grand Canonical Monte Carlo  . . . . . . . . . . . . . . . . . . . . . Simulation  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  122 122 122 125  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129  v  List of Tables 1.1  Mineral composition of North African dust collected at Sal Island (in the Cape Verde Islands), Barbados, and Miami. . . . . . . . . . . . . . . . . .  2  3.1  Binding energies and hydrogen-bond numbers for single water molecules obtained using NV T Monte Carlo simulations at 10 K. The numbers in brackets represented one standard deviation. . . . . . . . . . . . . . . . . . 44  4.1  Dimensions of the simulation cells in terms of the lattice parameters, a = 5.1535 ˚ A, b = 8.9419 ˚ A, and c = 7.3906 ˚ A. The y dimension and trench ˚ ˚ depth are 4b (35.77 A) and 4.5a (23.19 A), respectively, for all systems. . . 66 Average configurational energies at 298 K for trench 2. The numbers in brackets are estimated standard deviations. . . . . . . . . . . . . . . . . . . 76  4.2 5.1  Calculated dipole order parameters (obtained for the water molecules included either inside or outside the trench) and hydrogen bond numbers (with an angle condition of either 35o or 40o ). The associated standard deviations are included in brackets. . . . . . . . . . . . . . . . . . . . . . . 101  A.1 Nonbonded parameters for kaolinite (Al2 Si2 O5 (OH)4 ) using the CLAYFF Force Field and water models. The symbols and subscripts wo, wh, and lp represent water oxygen, water hydrogen, and lone pair, respectively. . . . . 116  vi  List of Figures 1.1 1.2  Structure of kaolinite (Al2 Si2 O5 (OH)4 ). . . . . . . . . . . . . . . . . . . . . SEM image of kaolinite particles covering a large quartz grain. . . . . . . .  2.1 2.2  A sketch of the simulation cell. . . . . . . . . . . . . . . . . . . . . . . . Blown up views of the (a) Al- and (b) Si- surfaces, (c) (100) unprotonated and (d) protonated edges. . . . . . . . . . . . . . . . . . . . . . . . . . . Adsorption isotherms obtained at 235 K using the SPC/E and TIP5P-E water models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Top view of a monolayer on the Al-surface obtained at 235 K with the TIP5P-E model and a 2-D projection in the x-y plane of the hexagonal arrangement of the water molecules and the underlying hydroxyl groups forming the Al-surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . .  2.3 2.4  3.1  3.2  3.3  3.4  3.5  3.6  4.1 4.2  The surface-water interaction energy of a single, randomly oriented, SPC/E water molecule, as a function of the perpendicular distance from the center of the highest hydroxyl hydrogen atom on the Al-surface. . . . . . . . . . Adsorption isotherms at 298 K for the four surfaces considered, obtained with a slab separation of 30 ˚ A. Results for the SPC/E (left) and TIP5P-E (right) models are shown. . . . . . . . . . . . . . . . . . . . . . . . . . . Adsorption isotherms obtained for the unprotonated edge, the protonated edge, and the Al-surface using the SPC/E and TIP5P-E models, at 298 K and 235 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Snap shots of a water molecule at 10 K for the Si-surface, Al-surface, unprotonated edge, and protonated edge using the SPC/E and TIP5P-E water models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Snap shots obtained at 235 K using the TIP5P-E model for low water coverages on the unprotonated edge and the protonated edge, and for monolayer coverages on the unprotonated edge and the Al-surface. . . . . . . . . . . The water density (g/cm3 ) profile and the average number of hydrogen bonds per water molecule at 235 K, obtained for monolayer coverage using the SPC/E and TIP5P-E models. . . . . . . . . . . . . . . . . . . . . . .  6 7  . 17 . 18 . 21  . 25  . 37  . 40  . 43  . 46  . 49  . 51  The simulation cell. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 The trenches considered at 298 K. Note that the systems are replicated along the x axis to match the dimension of the largest case. From top to bottom, the trenches are 1, 2, 3 and 4. Simulations were performed for trench 1 without (1A) and with (1B) the inter-trench spacing filled (clay atoms inside the black rectangle). The black arrows indicate the trenches. . 63 vii  List of Figures 4.3  4.4  4.5  4.6  4.7  4.8  5.1  5.2 5.3  5.4  5.5  The interaction energy averaged over the orientations of a single water molecule as a function of distance from the bottom of trench 1A (left) and 1B (right). The red and blue curves correspond to a water molecule moving up along both trench walls, while the green curve corresponds to a water molecule moving up the middle of the trench. . . . . . . . . . . . . . . . Adsorption isotherms obtained at 298 K for trenches 1A (red triangles up), 1B (orange stars), 2 (black circles), 3 (green squares), and 4 (maroon plus signs) with a slit separation of 50 ˚ A, as functions of RH. Results for an atomistically smooth unprotonated edge surface [36] (blue triangles down) are included for comparison. The results plotted as purple stars were obtained for trench 1B at a slit separation of 100 ˚ A. . . . . . . . . . Snapshots of trench 1A (left) and 1B (right) as functions of RH at 298 K. The corresponding chemical potentials are µ = −90 kJ/mol (top), −70 kJ/mol (middle), and −59 kJ/mol (bottom). The water oxygen atoms are blue, lattice oxygen atoms are red, hydrogen atoms are white, aluminum atoms are grey, and silicon atoms are brown. . . . . . . . . . . . . . . . . Snapshots of trench 2 as functions of RH at 298 K. The corresponding chemical potentials are µ = −90 kJ/mol (top), −70 kJ/mol (middle), and −59 kJ/mol (bottom). The atoms are denoted as in Fig. 4.5. . . . . . . . Snapshots of trench 3 as functions of RH at 298 K. The corresponding chemical potentials are µ = −90 kJ/mol (top), −64 kJ/mol (middle), and −59 kJ/mol (bottom). The atoms are denoted as in Fig. 4.5. . . . . . . Snapshots of trench 4 as functions of RH at 298 K. The corresponding chemical potentials are µ = −64 kJ/mol (top), −62 kJ/mol (middle), and −59 kJ/mol (bottom). The atoms are denoted as in Fig. 4.5. . . . . . . . Snapshots of trenches 1A (a) and 1B (b) using SPC/E at 235 K and ∼ 60 % RH. The water oxygen atoms are blue, lattice oxygen atoms are red, hydrogen atoms are white, ”lone pairs” are orange, aluminum atoms are grey, and silicon atoms are brown. . . . . . . . . . . . . . . . . . . . . . . Snapshots of trench 1B, 3 and 4 using SPC/E (∼ 60 % RH) and TIP5P-E (∼ 77 % RH) at 235 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . The water density (g/cm3 ) profile obtained at 235 K along the x axis (trench width) for trenches 1B (top panel), 3 (middle panel) and 4 (bottom panel) using the SPC/E (solid black curves) and TIP5P-E (dashed red curves) water models at ∼ 60 % and ∼ 77 % RH, respectively. Note the different scales on the x axis for all three trenches. . . . . . . . . . . . . . The water density (g/cm3 ) profile obtained at 235 K along the y axis for trenches 1B, 3, and 4 (from top to bottom) using the SPC/E (solid black curves) and TIP5P-E (dashed red curves) water models at ∼ 60 % and ∼ 77 % RH, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . The water density (g/cm3 ) profile obtained at 235 K along the z axis (trench depth) for trenches 1B, 3, and 4 (from top to bottom) using the SPC/E (solid black curves) and TIP5P-E (dashed red curves) water models, at ∼ 60 % and ∼ 77 % RH, respectively. . . . . . . . . . . . . . . . . . . . .  . 67  . 69  . 71  . 73  . 74  . 75  . 89 . 90  . 92  . 93  . 94 viii  List of Figures 5.6  Images of the individual layers along the x axis for trench 1B using SPC/E at 235 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Images of the individual layers along the z and y axis for trench 1B using SPC/E at 235 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Top and bottom views of the first layer above the top unprotonated edge at 235 K using SPC/E for trench 1B. . . . . . . . . . . . . . . . . . . . . 5.9 Profiles of the average dipole order parameter vector components obtained along the z axis starting from the bottom of trench 1B at 235 K using SPC/E at ∼ 60 % RH (black curves) and TIP5P-E at ∼ 77 % RH (red curves). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Profiles of the average magnitude of the dipole order parameter obtained along the z axis starting from the bottom of the trench at 235 K for trenches 1B (red curve), 3 (black curve), and 4 (blue curve) using SPC/E (bottom panel) and TIP5P-E (top panel) at ∼ 60 % and ∼ 77 % RH, respectively. 5.11 Profiles of the average dipole order parameter vector components obtained along the z axis for 1 water molecule near the Si-surface (bottom panel), middle of the trench (middle panel) and near the Al-surface (top panel) for trench 1B at 235 K using the SPC/E water model. The black, red and blue curves represent the x, y and z components of the order parameter, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  . 97 . 98 . 99  . 102  . 103  . 104  A.1 Schematics of the three-site SPC/E and five-site TIP5P-E water models. . 117  ix  Acknowledgements When I started my PhD I had a good idea of the task at hand, what would be the different stages, what would be the easy and the hard parts of the process. But I wasn’t really prepared for the real difficulties, the loneliness, the lack of motivation. Many thanks to my supervisors, Professors Gren Patey and Allan Bertram who guided and supported me all the way through my program here at UBC. Without them the completion of this thesis would not have been possible. At home, I was fortunate enough to have the most wonderful person to live with. My girlfriend, Virginie, put so much effort in supporting me that she would deserve to have her name on the front page of this thesis. Without her love and constant attention I would have been long gone. I also want to thank the members of my lab for their help. I want to thank my relatives and friends who encouraged me all along these past three years. Finally, thanks to the Westgrid facility and its managers who made the large number of computations possible and a priority.  x  Dedication To Virginie  xi  Statement of Co-Authorship • Chapters. 2, 3, 4, and 5 are coauthored with Prof. A. K. Bertram and Prof. G. N. Patey. The design of the research programs were developed by T. Croteau and Professors Bertram and Patey. The research and results are due to T. Croteau with guidance and suggestions from Professors Bertram and Patey. The final manuscripts were written by T. Croteau with revisions and additions by Professors Bertram and Patey.  xii  Chapter 1 Introduction 1.1  Mineral Dust  Mineral dust particles are abundant in the atmosphere, with estimated total global dust emissions to the atmosphere ranging from 800-2000 Tg/year (megatonne/year) [1–6]. Arid and semi-arid regions located around the globe are major contributors to the global dust emissions. In fact, most sources are found in the Northern hemisphere along the west coast of North Africa, the Middle East, Central and South Asia, and China [7], forming a region known as the “dust belt”. Outside of this region large-scale dust activity is considered remarkably small. Large dust storms occurring over the desert and arid regions of the “dust belt” are responsible for bringing and transporting mineral dust into the atmosphere. Monitoring of dust events has shown that these particles can be transported over large distances with atmospheric lifetimes on the order of days for small particles 0.1 − 10µm in diameter [1, 8]. Table 1.1 lists the mineral composition found in aerosolized dust coming from North Africa and collected in Sal Island in the Cape Verde Islands, Barbados and Miami [9]. Mica and quartz are usually major contributors while kaolinite, calcite, chlorite and other clay mineral are found in a lesser concentration. This thesis focuses on kaolinite particles which make up between 5 − 10 % of the aerosolize mass of mineral dust samples [9].  1  1.2. Importance of Dust  Table 1.1: Mineral composition of North African dust collected at Sal Island (in the Cape Verde Islands), Barbados, and Miami [9]. Mineral Sal Island Barbados Miami Mica 53.8% 64.3% 62.1% Quartz 19.6% 13.8% 14.2% Kaolinite 6.6% 8.3% 7.1% Calcite 8.2% 3.9% 6.9% Plagioclase 5.4% 4.1% 4.5% Chlorite 4.3% 4.1% 4.2% Montmorillonite ≤ 5% ≤ 5% ≤ 5%  1.2 1.2.1  Importance of Dust Role of Mineral Dust in Atmospheric Chemistry  Mineral dust particles can play an important role in the Earth’s climate and atmospheric chemistry. It is believed that the presence of mineral dust particles in the troposphere may provide a surface for heterogeneous reactions to occur with trace atmospheric gases. This can affect the concentration of trace atmospheric species such as N2 O5 , nitric acid and ozone. Removal of gas phase atmospheric species such as nitrogen oxides can have important implications as they are involved in the complex system of reactions leading to the formation of ozone [10]. Processes occurring on surfaces could possibly explain some of the reaction mechanisms that cannot yet be answered simply based on a gas-phase reaction point of view. In addition, it has been shown that the interaction of mineral dust with trace atmospheric species is influenced by the relative humidity. For example, an increase in uptake and rate of heterogeneous hydrolysis of N2 O5 on quartz [11], an increase in uptake of nitric acid on Na-Montmorillonite [12] and dust samples [10, 13, 14], and a decrease in uptake and dissociation of ozone on alumina and hematite surfaces [11] was reported as the relative humidity was increased. These results show that it is necessary to understand the interaction of water with mineral surfaces if we want to understand these reactions under atmospheric conditions.  2  1.2. Importance of Dust  1.2.2  Role of Mineral Dust in Climate: Direct and Indirect Effect  Mineral particles can affect climate on a local and global scale by adsorbing and scattering radiation. This can occur either directly by the particles themselves or indirectly through the modification of the reflective properties of clouds. The effect that these processes have on climate is highly uncertain. The net radiative forcing induced by mineral aerosol was estimated to range from -0.60 to +0.40 Wm−2 (watts per square metre) by the Intergovernmental Panel on Climate Change (IPCC) [15]. Thus, the impact can either be a cooling or a warming effect. The role played by mineral particles on climate is highly uncertain due to the lack of understanding of the above chemical processes and the physical properties of the particles [8]. Two important properties of mineral particles related to the indirect effect that need to be better understood are their cloud condensation and ice nucleation abilities. For example, in the atmosphere mineral dust particles can take up water and act as cloud condensation nuclei [16–22]. In the troposphere, they can also act as ice nuclei to form ice clouds. In both cases, either more precipitation (rain, snow or hail) or more clouds can be produced depending on atmospheric conditions and the type and abundance of aerosol particles present, which serve as cloud condensation nuclei (CCN) or ice forming nuclei (IN) [16, 23, 24]. This could result in an increase (more precipitation) or decrease (more clouds) of the amount of solar radiation reaching the surface of the Earth thus indirectly affecting climate. Understanding the interaction of H2 O with mineral dust surfaces is important for understanding the role of mineral dust in atmospheric chemistry and climate. Note that the role played in ice cloud formation by the atmospherically abundant mineral dust particles is thought to be important [25]. As mentioned above, we will focus on kaolinite particles. These particles are known to be especially effective ice nuclei based on laboratory studies [26]. We expect that the information determined for kaolinite will also help understand 3  1.3. Kaolinite Structure other clay minerals.  1.3  Kaolinite Structure  Kaolinite [Al2 Si2 O5 (OH)4 ] is a clay mineral that has a 1:1 layered structure. Each layer is composed of a single silica tetrahedral sheet [Si2 O5 ] bonded to a single alumina octahedral sheet [Al2 (OH)4 ] called silicate and gibbsite sheets, respectively. In the remainder of this thesis we will refer to these two sheets in a different way. Based on their intrinsic nature and in agreement with an establish convention, the silicate and gibbsite sheets are named the Si- and Al-surfaces, respectively. The stacking of layers in the kaolinite mineral leads to a triclinic geometry with a, b, and c parameters of 5.1535 ˚ A, 8.9419 ˚ A, and 7.3906 ˚ A, respectively; α = 91.926◦; β = 105.046◦; γ = 89.797◦ ; and the space group is C1 [27]. The layers (basal planes) are attached together by hydrogen bonds between the hydroxyl hydrogens extending from the Al-surface and bridging oxygens on the Si-surface. The plane between the layers is called a cleavage plane. The large number of hydrogen bonds connecting the layers does not allow water to penetrate the cleavage plane, and kaolinite is described as a non-expanding or non-swelling clay for that reason. Fig. 1.1 shows the atomic structure of kaolinite and the hydrogen bonding between the layers. In addition to both basal (001) planes (Al- and Si-surfaces), kaolinite particles also contain several possible “edges” (other planes) which are believed to make a significant contribution to the surface area of kaolinite samples [28–30]. The flaky structure of a macroscopic sample is shown in Fig. 1.2, where the flat surfaces correspond to the basal surfaces. The edges are not always perpendicular to the basal surface. The edges of the flakes are somewhat irregular. Kaolinite particles are composed of several layers stacked together which leads to thick edges. Well crystallized kaolinite particles have a well defined hexagonal geometry with elongation in one direction, whereas poorly crystallized samples have less distinct six-sided flakes [31]. The Si-surface is known to be more hydrophobic than the Al-surface which is termi4  1.4. Thesis Objectives nated by hydroxyl groups and attracts water more readily [32]. Similarly, the edges of kaolinite are terminated by oxygen atoms or hydroxyl groups, which can easily be protonated or deprotonated depending on pH [30]. Two extreme edge cases are considered here, one which has “bare” oxygens exposed (unprotonated edge) and another where all oxygens are protonated (protonated edge). Both edge constructions are very hydrophilic.  1.4  Thesis Objectives  The overall goal of this thesis is to use computer simulations to better understand the interaction of water with mineral surfaces, and obtain new insights regarding their ability as cloud condensation and ice nuclei, as well as their role in heterogeneous chemistry. Using the grand canonical Monte Carlo method we look at the adsorption and structure of water on kaolinite surfaces at 235 K, a temperature relevant for ice nucleation in the atmosphere. We carry out these calculations for four different surfaces: the Al-surface (gibbsite), Si-surface (silicate) and two edges. This information is then used to gain insight into ice nucleation on kaolinite. We want to know if water adsorbs on defect-free kaolinite surfaces at 235 K and under typical atmospheric relative humidities. The effect of surface structure is also investigated, as it is believed that adsorbed water molecules will favorably adopt a structure similar to hexagonal ice when the surface structure itself matches that of ice Ih. We also determine water adsorption isotherms at 235 and 298 K on kaolinite surfaces, and the results are compared with recent laboratory studies. The water affinity of the different kaolinite surfaces is investigated by looking at the orientation and binding energies of single adsorbed water molecules at 10 K. This includes determination of the binding energies and geometric arrangements of single water molecules on the surface, as well as water-surface hydrogen bond analyses. Water structure in the monolayers is analyzed using density and hydrogen-bond number profiles. In addition to studying atomistically smooth kaolinite surfaces, we performed water 5  1.4. Thesis Objectives  Figure 1.1: Structure of kaolinite [Al2 Si2 O5 (OH)4 ]. The dotted lines show the hydrogen bonds formed between the gibbsite and silicate layers. The aluminum, silicon, oxygen, and hydrogen atoms are grey, brown, red, and white, respectively.  6  1.4. Thesis Objectives  Figure 1.2: SEM image of kaolinite particles covering a large quartz grain [33]. adsorption studies on kaolinite trenches, which are represented as infinitely long rectangular cavities in one dimension and of fixed width and depth. The same method as for flat surfaces is used, namely, the grand canonical Monte Carlo simulation method both at 235 and 298 K. The goal here is to find evidence for the presence of large quantities of water on clay minerals and make comparisons with experiment. More specifically, we want to know if the water coverages in the case of trenches are larger than the coverages observed for atomistically flat surfaces. We also want to determine what kind of effects does the size of the trenches have on adsorption, and if our results can help explain the high coverages observed in recent laboratory studies [34]. Analysis of the structures adopted by the water molecules and the binding energies inside the trenches is also of interest as it may well be the key to understanding ice nucleation on clay minerals.  7  1.5. Review of Simulation  1.5  Review of Simulation  Molecular simulations are a well suited technique for studying the interactions of water with mineral surfaces. It allows the acquisition of useful thermodynamical and statistical data unavailable from experiments. For example, the hydrophobic and hydrophilic nature of the Si- and Al-surfaces of kaolinite was confirmed in several ab initio studies [35–39]. Note that the term hydrophobic here refers to the rather weak interaction of a water molecule with the Si-surface compared to other surfaces such as the Al-surface. Previous theoretical studies of ice nucleation on surfaces have shown the possible formation of hexagonal rings but rarely the creation of an ordered multilayer ice crystal. For example, Shevkunov [40, 41] studied the structure of water on the surface of a β-AgI crystal as well as in microcracks using a Monte Carlo approach. In his simulations, the water molecules organized themselves into monolayer islands on a single surface or into multilayered films in a microcrack with clear features of hexagonal ice. Similar patterns were obtained by Taylor and Hale [42] also on a AgI substrate, although, they found that the ice-like character did not extend beyond the first water layer. Using molecular dynamics to study different modes of water aggregation on the surface of CaF2 and BaF2 , Wassermann et al. [43] found that at high coverage a double layer ice-Ih structure was formed on the BaF2 surface. Recently, Liang Hu and Michaelides [39] reported that water monomers bind strongly onto the Al-surface and form hexagonal rings. However, what they found does not agree with the idea that the good ice nucleation ability of kaolinite is related to the good crystallographic match between kaolinite surfaces and hexagonal ice surfaces [44]. In fact, they note a mismatch between the ice Ih lattice and the 2D water structure imposed by the substrate. They also found that multilayer ice growth is not favored, and describe the monolayer covered kaolinite surface as hydrophobic. However, as these results were not obtained in an open system it is not certain if these structures are stable under normal atmospheric conditions. A relevant exception is the work of Delville [45], who focused on the Si-surface. Other simulation studies of water on various related 8  1.6. Thesis Overview surfaces include metals [46, 47], and model hexagonal surfaces [48]. Although the edges contribute significantly to the total surface area of clay minerals and other dust particles very few simulation studies actually investigated the effect of surface irregularities. One example is the study of Cailliez et al who used grand canonical Monte Carlo simulations (GCMC) to investigate the effect of surface defects on the hydration thermodynamic behavior of silicalite-1 zeolite [49]. They found that as you increase the number of surface defects the hydrophilic character of the surface also increases. Similarly, Shevkunov [50] found that crystal defects with size 15-20˚ A stimulated more bulk condensation than small ones when looking at the adsorption of water on an irregular Ag-I crystal using a Monte Carlo method.  1.6  Thesis Overview  This thesis is presented in six Chapters and one Appendix. The current chapter dealt with background information relevant to mineral dust particles, their atmospheric implications, and the need to understand the interaction of water with kaolinite and other surfaces. Chapters 2 and 3 focus on water structure and water adsorption on atomistically flat kaolinite surfaces, respectively. Chapter 4 describes water adsorption at 298 K for the case of kaolinite trenches. Chapter 5 considers possible ice nucleation at 235 K in kaolinite trenches. Finally, a summary of the results and possible future directions are discussed in Chapter 6. Appendix A describes in detail the simulation method and models used.  9  Bibliography [1] Bauer, S. E.; Balkanski, Y.; Schulz, M.; Hauglustaine, D. A.; Dentener, F. J. Geophys. Res.-Atm. 2004, 109 (D02304). [2] Lunt, D. J.; Valdes, P. J. J. Geophys. Res.-Atm. 2002, 107 (D23). [3] Tegen, I.; Harrison, S. P.; Kohfeld, K. E.; Prentice, I. C.; Coe, M.; Heimann, M. J. Geophys. Res.-Atm. 2002, 107 (D21), 4576. [4] Tegen, I.; Miller, R. J. Geophys. Res.-Atm. 1998, 103 (D20), 25975. [5] Woodward, S. J. Geophys. Res.-Atm. 2001, 106 (D16), 18155. [6] Zender, C. S.; Bian; H.; Newman, D. J. Geophys. Res. 2003, 108 (D14), 4416. [7] Prospero, J. M.; Ginoux, P.; Torres, O. ; Nicholson, S. E.; Gill, T. E. Rev. Geophys. 2002, 40 (1), 1002. [8] Husar, R. B.; et al. J. Geophys. Res. 2001, 106, 18317. [9] Glaccum, R. A.; Prospero, J. M. Marine Geology 1980, 37 (3-4), 295. [10] Underwood, G. M.; Li, P.; Al-Abadleh, H.; Grassian, V. H. J. Phys. Chem. A 2001, 105, 6609. [11] Mogili, P. K.; Kleiber, P. D.; Young, M. A.; Grassian, V. H. J. Phys. Chem. A 2006, 110, 13799. [12] Mashburn, C. D.; Frinak, E. K.; Tolbert, M. A. J. Geophys. Res.-Atm. 2006, 111 (D15). [13] Hanisch, F.; Crowley, J. N. Phys. Chem. Chem. Phys. 2001, 3 (12), 2474. [14] Laskin, A.; Wietsma, T. W.; Krueger, B. J.; Grassian, V. H. J. Geophys. Res. 2005, 110, D10208. [15] Forster, P.; Ramaswamy, V.; Artaxo, P.; Berntsen, T.; Betts, R.; Fahey, D. W.; Haywood, J.; Lean, J.; Lowe, D. C.; Myhre, G.; Nganga, J.; Prinn, R.; Raga, G.; Schulz, M.; Dorland, R. V. Changes in Atmospheric Constituents and in Radiative Forcing., in Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change 2007, edited by S. Solomon, D. Qin, M. Manning, Z. Chen, M. Marquis, K.B. Averyt, M. Tignor, and H.L. Miller, pp. 129-234, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA. 10  Chapter 1. Bibliography [16] Cziczo, D. J.; Murphy, D. M.; Hudson, P. K.; Thomson, D. S. J. Geophys. Res.-Atm. 2004, 109(D4210), 4517. [17] DeMott, P. J.; Cziczo, D. J.; Prenni, A. J.; Murphy, D. M.; Kreidenweis, S. M.; Thomson, D. S.; Borys, R.; Rogers, D. C. Proc. Nat. Acad. Sci. 2003, 100 (25), 14655. [18] Eastwood, M. L.; Cremel, S.; Gehrke, C.; Girard, E.; Bertram, A. K. J. Geophys. Res.-Atm. 2008, 113. [19] Dymarska, M.; Murray, B. J.; Sun, L.; Eastwood, M. L.; Knopf, D. A.; Bertram, A. K. J. Geophys. Res.-Atm. 2006, 111 D4. [20] Zuberi, B.; Bertram, A. K.; Cassa, C. A.; Molina, L. T.; Molina, M. J. Geophys. Res. Let. 2002, 29 (10) [21] Twohy, C. H.; Kreidenweis, S. M.; Eidhammer, T.; Browell, E. V.; Heymsfield, A. J.; Bansemer, A. R.; Anderson, B. E.; Chen, G.; Ismail, S.; DeMott, P. J.; Van den Heever, S. C. Geophys. Res. Let. 2009, 36. [22] Twohy, C. H.; Poellot, M. R. Atm. Chem. Phys. 2005, 5, 2289. [23] Facchini, M.C.; Mircea, M.; Fuzzi, S.; Charlson, R. J. Nature 1999, 401(6750), 257. [24] Levin, Z.; Ganor, E.; Gladstein, V. J. Appl. Meteorol. 1996, 35(9), 1511. [25] Cantrell, W.; Heymsfield A. J. Bull. Amer. Meteor. Soc. 2005, 86(6), 795. DeMott, P. J.; Cziczo, D. J.; Prenni, A. J.; Murphy, D. M.; Kreidenweis, S. M.; Thomson, D. S.; Borys, R.; Rogers, D. C. Proc. Nat. Acad. Sci. 2003b, 100, 25, 14655. [26] Zimmerman, F.; Ebert, M.; Worringen, A.; Schtz, L.; Weinbruch, S. Atmos. Environ. 2007, 41, 8219. [27] Bish, D. L. Clays Clay Miner. 1993, 41, 783. [28] Brady, P. V.; Cygan, R. T.; Nagy, K. L. J. Colloid and Interface Science 1996, 183, 356. [29] Zbik, M.; Smart, R. S. C. Clays and Clay Minerals 1998, 46(2), 153. [30] Wang, Y.-S.; Siu, W.-K. Can. Geotech. J. 2006, 43, 587. [31] Grim, R. E. Clay Mineralogy 1953, edited by McGraw-Hill Book Company, New York, USA. [32] Tunega, D.; Gerzabek, M. H.; Lischka, H. J. Phys. Chem. B 2004, 108, 5930. [33] Barthelmy, D. Kaolinite image (http://webmineral.com/specimens/picshow.php?id=1283) 2009, copyright Weatherford Labs, http://www.weatherfordlabs.com. [34] Schuttlefield, J. D.; Cox, D.; Grassian, V. H. J. Geophys. Res. 2007, 112 (D21303). 11  Chapter 1. Bibliography [35] Tunega, D.; Gerzabek, M. H.; Lischka, H. ıJ. Phys. Chem. B 2004, 108, 5930. [36] Tunega, D.; Benco, L.; Haberhauer, G.; Gerzabek, M. H.; Lischka, H. J. Phys. Chem. B 2002, 106, 11515. [37] Tunega, D.; Haberhauer, G.; Gerzabek, M. H.; Lischka, H. Langmuir 2002, 18, 139. [38] Hu, X. L.; Michaelides, A. Surf. Sci. 2007, 601, 5378. [39] Hu, X. L.; Michaelides, A. Surf. Sci. 2008, 602, 960. [40] Shevkunov, S. V. Colloidal Journal 2005, 67, 497. [41] Shevkunov, S. V. Colloidal Journal 2007, 69, 360. [42] Taylor, J. H.; Hale, B. N. Phys. Rev. B 1993, 47, 9732. [43] Wassermann, B.; Reif, J.; Matthias, E. Phys. Rev. B 1994, 50, 2593. [44] Pruppacher, H. R.; Klett, J. D. Microphysics of Clouds and Precipitation, edited by Kluwer Academic, Dordrecht, 2nd rev. ed., 1997, p. 329-331. [45] Delville, A. J. Phys. Chem. 1995, 99, 2033. [46] Cerd´a, J.; Michaelides, A.; Bocquet, M.-L.; Feibleman, P. J.; Mitsui, T.; Rose, M.; Fomin, E.; Salmeron, M. Phys. Rev. Lett., 2004, 93, 116101-1. [47] Michaelides, A.; Morgenstern, K. Nature Mater, 2007, 6, 597. [48] Nutt, D. R.; Stone, A. J. Langmuir 2004, 20, 8720. [49] Cailliez, F; Stirnemann, G.; Boutin, A.; Demachy, I.; Fuchs, A. H. J. Phys. Chem. C 2008, 112, 10435. [50] Shevkunov, S. V. Colloidal Journal 2006, 68, 370.  12  Chapter 2 Adsorption and Structure of Water on Kaolinite Surfaces: Possible Insight into Ice Nucleation from Grand Canonical Monte Carlo Calculations∗ 2.1  Introduction  Ice clouds play an important role in the Earth’s radiation budget by reflecting and absorbing solar and terrestrial radiation [1]. Ice particles, and hence ice clouds, form when ice nucleates on or in aerosol particles present in the atmosphere. This can occur either homogeneously (in a liquid aerosol particle) or heterogeneously (on the surface of a solid aerosol particle such as mineral dust or soot). Heterogeneous ice nucleation has been studied for many years, yet it remains poorly understood, especially on a molecular level [2]. Mineral dust particles, which are abundant in the atmosphere, are thought to play an important role in ice cloud formation in the atmosphere [3]. Common minerals found in aerosolized dust include quartz, illite, muscovite, chlorite, kaolinite, and calcite [4]. In the following we will focus on kaolinite particles which make up between 5 − 10 % of the aerosolized mass of mineral dust [4]. In addition, kaolinite is known to be an especially effective ice nucleus based on laboratory studies [5–7]. Kaolinite is a type of clay particle and information determined on kaolinite may also lead to a better understanding of other A version of this chapter has been published. T. Croteau, A. K. Bertram, G. N. Patey (2008), Adsorption and Structure of Water on Kaolinite Surfaces: Possible Insight into Ice Nucleation from Grand Canonical Monte Carlo Calculations, J. Phys. Chem. A, 112:10708-10712. ∗  13  2.1. Introduction clays. Kaolinite particles are composed of different surface types, two basal (001) planes and several possible “edges” (other planes). The basal (001) planes in kaolinite are referred to as the Al-surface and Si-surface. The Si-surface is known to be hydrophobic, while the Al-surface, which is terminated by hydroxyl groups, attracts water [8]. Similarly, the edges of kaolinite are terminated by oxygen atoms or hydroxyl groups which can easily be protonated or deprotonated depending on pH [9]. The edges used here are termed protonated and unprotonated edges and are surfaces parallel to the (100) plane. These are hydrophilic surfaces that may play a role in ice nucleation or at least in water adsorption. The reason kaolinite is a good ice nucleus remains unclear. In has been speculated that the good ice nucleation ability of kaolinite is related to the good crystallographic match between kaolinite surfaces and hexagonal ice surfaces. This implies that kaolinite surfaces have sites that arrange the water molecules in a structure that has similar lattice constants to hexagonal ice. The absorbed water structures may then be good templates for ice embryo formation, since the resulting embryos may have minimal strain due to the good match between the adsorbed water at the kaolinite surface and hexagonal ice. In order to quantify the crystallographic match between kaolinite and hexagonal ice, Pruppacher and Klett [2] calculated an apparent crystallographic misfit or disregistry using δ≡  nao,N − mao,i , mao,i  where ao,N is the crystallographic lattice parameter of a particular face of the nucleus, ao,i is the corresponding constant in the ice lattice and n, m are integers chosen such that δ is minimal. Pruppacher and Klett found a good match (δ = 1.1%) between the basal face of hexagonal ice and the Al-surface of kaolinite if n and m are 3 and 2, respectively. These calculations employed a hypothetical “pseudohexagonal” surface with an ao lattice parameter of 2.98 to represent the Al-surface, and an ao lattice parameter of 4.52 for 14  2.2. Model and Method hexagonal ice. In the present study we investigate the adsorption and structure of water on kaolinite surfaces at 235 K, a temperature relevant for ice nucleation in the atmosphere. Calculations are carried out for four different surfaces, the Al-surface, the Si-surface and protonated and unprotonated edges. The information obtained is used to gain insight into ice nucleation on kaolinite. Key questions we address are the following: 1) Does water adsorb on defect-free kaolinite surfaces at 235 K and under typical atmospheric relative humidities? 2) Do the kaolinite surfaces promote adsorbed water molecules to adopt a structure similar to hexagonal ice? If so, these water structures could serve as good building blocks for ice embryo formation since they would impose minimal strain. Recently, Hu and Michaelides [10, 11] used DFT calculations to examine the structure of water on the Al-surface of kaolinite, and the possible implications for ice nucleation. We compare with these previous calculations. Several other groups have also looked at water adsorption on different surfaces in the context of ice nucleation. These calculations are briefly discussed as well.  2.2  Model and Method  In the present calculations, we consider two different water models, the extended single point charge SPC/E [12] and TIP5P-E [13]. Kaolinite is described with the Clay Force Field (CLAYFF) recently developed by Cygan et al. [14]. SPC/E is a rigid water model with point charges located at the three atomic nuclei, and all are embedded in a LennardJones (LJ) sphere. This model has been widely used in the past and is known to give accurate results in the liquid state. However, it was found [15] that the SPC/E model does not accurately predict the melting point of ice. In comparison, the TIP5P model of Mahoney and Jorgensen [16] gives a value of 274 K for the melting point. Therefore, the TIP5P-E [13] model, which is a recent reparametrization of the original TIP5P [16] model to account properly for long-range interactions, is also used. In this model, equal 15  2.2. Model and Method positive and negative partial charges are placed at the center of both hydrogen nuclei and at oxygen “lone pairs” locations, respectively, and are embedded in a LJ sphere centered at the oxygen nucleus. The use of two different water models is useful as it gives an indication of the sensitivity of our results system to the water model employed. CLAYFF is an atomistic model again constructed of LJ spheres and point charges. In its general form the CLAYFF model can be flexible, but for practical purposes we employ a rigid lattice set in the stable kaolinite structure as obtained by Bish [17]. The total interaction energy for our model is given by the sum of Coulombic and van der Waals terms 1 U= 4πǫo  i=j  qi qj + rij  Do,ij i=j  Ro,ij rij  12  Ro,ij −2 rij  6  ,  (2.1)  where qi is the charge on site i, ǫo is the permittivity of free space, and we have used the notation of Ref. [14]. The summations are over all interaction sites in the system including water molecules and the kaolinite lattice. All short-range interactions in the van der Waals term are represented by the LJ potential. The parameters involved in the kaolinite-water interactions are given in Ref. [14]. The simulation cell employed (Fig. 2.1) is similar to that used in earlier work involving two slabs and long-range electrostatic interactions [18, 19]. The basic cell contains one of the four surfaces displayed in Fig. 2.2, and is repeated periodically in the x and y dimensions. The (x, y) cell dimensions are (30.921 ˚ A, 35.7676 ˚ A) for the Al- and Sisurfaces, and (29.524 ˚ A, 35.7676 ˚ A) for the edges. These “single-layer” slabs consist of 816 atoms (sites) which corresponds to the translation of 48 kaolinite unit cells obtained from the ideal formula Al2 Si2 O5 (OH)4 . The short-range interactions were calculated employing a spherical cutoff of half the box length, and the total Coulombic energy was obtained using the usual threeA was found to be dimensional Ewald method [18–20]. An empty space (Fig. 2.1) of 107 ˚ sufficient to avoid unphysical interactions with images in the z direction. Most simulations were performed with the two slabs separated by 30 ˚ A. Some were 16  2.2. Model and Method  1 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 1111111111111111111 0000000000000000000  EMPTY SPACE  WATER  Slab  z (001) (010) y (100) x  1111111111111111 0000000000000000 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 EDGE 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111  1111111 0000000 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000000000000000 1111111111111111111 0 1 0000000 1111111 0000000000000000000 1111111111111111111 0 1 0000000 1111111 0000000000000000000 1111111111111111111 0 1 0000000 1111111 0000000000000000000 1111111111111111111 0 1 0000000000000000000 1111111111111111111 0 1 0000000000000000000 1111111111111111111 0 1 0000000000000000000 1111111111111111111 0 1 000000 111111 0000000000000000000 1111111111111111111 0 1 000000 111111 0000000000000000000 1111111111111111111 0 1 000000 111111 0000000000000000000 1111111111111111111 0 1 000000 111111 0000000000000000000 1111111111111111111 0 1 000000 111111 0000000000000000000 1111111111111111111 0 1 000000 111111 0000000000000000000 1111111111111111111 0 1 000000 111111 0000000000000000000 1111111111111111111 0 1 000000 111111 0000000000000000000 1111111111111111111 0 1 000000 111111 0000000000000000000 1111111111111111111 0 1 0 1 000000 111111 0000000000000000000 1111111111111111111 0 1 0 1 000000 111111 0 1 0 1 000000 111111 0 1 0 1 000000 111111 0 1 0 1 000000 111111 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 1111111111111111111 0000000000000000000  EMPTY SPACE  Figure 2.1: A sketch of the simulation cell.  17  2.2. Model and Method  (a)  (b)  (c)  (d)  Figure 2.2: Blown up views of the (a) Al- and (b) Si- surfaces, (c) (100) unprotonated and (d) protonated edges. The oxygen, hydrogen, silicon and aluminum atoms are red, white, brown and grey, respectively.  18  2.2. Model and Method performed at a separation of 60 ˚ A to verify that water adsorption at one surface is not significantly influenced by the presence of the other. Thus, before the “filling transition” at saturation (corresponding to bulk condensation of the vapor) the system behaves as two surfaces adsorbing independently, and results for both surfaces can be averaged to improve the statistics. Calculations were also performed with thicker kaolinite slabs, and adding an extra layer was shown to have no significant effect on water adsorption. Edge configurations (Fig. 2.2c and 2d) were constructed by cleaving the Al-surface (Fig. 2.2a) parallel to the (100) plane, and replicating in the xy plane. The configuration thus obtained represents one limiting case with pH above the isoelectric point (IEP) such that the surface is completely deprotonated [9]. At the other extreme, well below the IEP, all of the surface oxygens would be protonated such that the surface is positively charged. The protonation state of the edges under atmospheric conditions is not known, and is likely variable, therefore, we examine both extremes. We find in fact that the water adsorption ability of an edge is not particularly sensitive to the protonation state and is similar for both limiting cases considered (see Fig. 2.3). The protonated edge is modelled by placing protons (with partial charges appropriate to hydroxyl groups [14]) directly above every “bare” oxygen atom with an internuclear distance of 1 ˚ A. To recover global neutrality in our sample, the extra positive charges are countered by distributing an equal negative charge over all of the atoms below the top Al-Si layer. This amounts to adding a small charge (−0.05 e, e is an elementary charge) to each of the 544 atoms below the surface. Again, changing the slab thickness (i.e., reducing the negative charge added per subsurface atom) had no significant effect on the results. All results reported were obtained employing grand canonical Monte Carlo calculations [20]. In this method the thermodynamic state is fixed by specifying the volume, V , the temperature, T , and the water chemical potential µ. The number of water molecules in the system, NW , is allowed to fluctuate. A Monte Carlo step consisted of a particle insertion, deletion, translation or rotation. In all simulations, the system was equilibrated  19  2.3. Results and Discussion for at least 5 × 107 MC steps. Following equilibration, averages were collected for an additional 2 × 108 or more MC steps.  2.3 2.3.1  Results and Discussion Adsorption Isotherms  Adsorption isotherms in the form of the average number of water molecules, < NW >, versus the chemical potential, µ, were obtained for all four surfaces. The 235 K adsorption isotherms cast in a more familiar form are plotted in Fig. 2.3. In this plot, the “water coverage” is < NW > AW /A, where AW is the surface area occupied by a water molecule (taken to be 9 ˚ A2 ) and A is the total surface area of the particular face considered (A = Lx Ly , where Lx and Ly are the x and y dimensions of the simulation cell). The pressure, P , of the water vapor at equilibrium corresponding to particular chemical potentials was estimated using the ideal gas relationship, P = eβµ /βΛ3 , where Λ = (βh2 /2πm)1/2 is the thermal de Broglie wavelength, β = 1/kB T , kB is the Boltzmann constant, and m is the mass of a water molecule. P0 is the saturated vapor pressure chosen such as to set P/P0 = 1 at the observed filling transition. It is important to remark that the < NW > versus µ isotherms directly obtained in our simulations are accurate and are in no way influenced by the pressure estimates used in the plot. Accurate vapor pressures could in principle be obtained from simulation, but the ideal gas equation is sufficiently accurate to give a qualitative indication of the pressures involved. For the Si-surface, we found no adsorption before saturation (condensation), and hence no results for this surface are included in Fig. 2.3. This is not surprising since the Sisurface is known to be hydrophobic, and our result is consistent with earlier work [21]. An interesting observation is that both the protonated and unprotonated edges take up water at lower pressures than does the Al-surface, with the protonated edge having the greatest affinity for water. Nevertheless, the relative pressure range, P/P0 , at which adsorption 20  2.3. Results and Discussion  water coverage  2.5  2  1.5  1  0.5  0  0  0.2  0.4  0.6  0.8  1  0.2  0.4  0.6  0.8  1  P/Po Figure 2.3: Adsorption isotherms obtained at 235 K using the SPC/E (left) and TIP5PE (right) water models. Results are shown for the unprotonated edge (black circles), protonated edge (red squares), and Al-surface (green triangles).  21  2.3. Results and Discussion occurs on the Al-surface (Fig. 2.3) is significant for both models, TIP5P-E predicting the largest range. These findings are important because they strongly suggest that edges can contribute significantly to kaolinite’s ability to take up water, and perhaps to the heterogeneous chemistry where water is needed for hydrolysis [22, 23]. The results shown in Fig. 2.3 were obtained by increasing the chemical potential from a low value (empty simulation cell) to filling. All surfaces were checked for hysteresis. This was done by starting simulations from a configuration obtained at the last point before filling and decreasing the chemical potential until all water had evaporated. For both edges the curves obtained (not shown) are relatively smooth with little hysteresis. This indicates that the strong surface-water interactions are largely responsible for the adsorption, and there isn’t much suggestion of collective behavior. The situation is different for the Al-surface. Here one notes that the surface coverage (see Fig. 2.3) “jumps” quite quickly from essentially nothing to a monolayer for both SPC/E and TIP5P-E. For the SPC/E at 235 K, this is followed by a second layer and some additional thickening before filling. This behavior was verified by performing simulations with the surface separation increased to 60 ˚ A (not shown). Additionally, for the Al-surface significant hysteresis (not shown) is observed. These observations indicate that collective behavior amongst the water molecules is important for adsorption to occur on the Al-surface. In other words, water condensation onto the Al-surface exhibits features characteristic of first-order processes, whereas condensation onto the edges appears essentially continuous. A detailed analysis of the water-surface interactions for different surfaces and configurations will be given in a forthcoming article. Here we simply report that the average adsorption energies (in kJ/mol) at 150 K for a single water are (-15.7, -15.6), (-43.1, -41.4), (-69.4, -69.3), (-89.7, -72.0) for the Si-surface, the Al-surface, the unprotonated edge and the protonated edge, respectively, where the first number of each pair is for the SPC/E model and the second number is the TIP5P-E result. We note that the SPC/E and TIP5P-E values are in reasonably good agreement except for the protonated edge  22  2.3. Results and Discussion where the model discrepancy is significant. Also, it is obvious that the single molecule binding energies are in accord with the adsorption isotherms discussed above.  2.3.2  Water Structure and Lattice Match  Fig. 2.4a shows a typical monolayer obtained on the Al-surface at 235 K using TIP5P-E. The water molecules are quite uniformly spread, and one can discern loose “networks” with the water molecules clearly interacting with surface atoms and with each other. It is observed that only for the Al-surface can strong correlations amongst the water molecules produce some hexagonal patterns which are easily discernable in Fig. 2.4a. The same patterns are found using SPC/E, but they are more evident for TIP5P-E. Well defined six-membered ring patterns are not observed on the edges. Instead, on the unprotonated edge, the water molecules primarily interact with the surface in very specific sites along the interlayer junctions (center of Fig. 2.2c) where two kaolinite sheets come into contact and form multiple hydrogen bonds between the basal Al- and Si-surfaces. On the protonated edge, the water molecules prefer to bind first along the hydroxyl hydrogen lines created by protonation (slightly off center in Fig. 2.2d). As more water is adsorbed a network is created across the surface for both edges. The hexagonal structure we obtain for kaolinite (Fig. 2.4b) does not match ice Ih. The ao lattice parameter that would be applicable to the hexagonal structures adsorbed on kaolinite is roughly 5.15, significantly different than 4.52, the ao lattice parameter for hexagonal ice. To estimate the effect of this difference between the ao lattice parameters we use the following equation to calculate the strain on the ice embryo due to mismatch [2] ′  ao,i − ao,i ǫ≡ , ao,i ′  where ao,i and ao,i are the lattice parameters of the ice in the strain-free and strained ′  conditions, respectively. Assuming ao,i = 4.52 and ao,i = 5.1535, we obtain a strain of 14.0%. Based on Turnbull and Vonnegut [24], a strain of more than 5% would lead to 23  2.3. Results and Discussion an additional depression in the nucleation temperature of at least 40o C, so it is hard to imagine that the hexagonal rings observed in our simulations will be beneficial for ice nucleation due to the possible strain imposed on the ice embryo. The crystallographic misfit calculations discussed by Pruppacher and Klett [2] (see Introduction) lead to very different conclusions than our GCMC calculations. Their calculations suggest that the actual mismatch between kaolinite and hexagonal ice is only 1.1%, which should result in a very small difference between the ao lattice parameters of the adsorbed water and hexagonal ice. However, these calculations assume that the water molecules will arrange themselves on the Al-surface with a 3:2 ratio of the ao lattice parameters. However, the hexagonal structure we obtain for kaolinite (Fig. 2.4b) does not match ice Ih, but rather fits the roughly hexagonal arrangement of the surface hydroxyl groups on the Al-surface of kaolinite in a 1:1 ratio. Recently, Hu and Michaelides [10, 11] used DFT calculations to study water binding to the Al-surface of kaolinite. They reported that water formed 2-D ice-like layers with a stability matching that of ice Ih, but they noted that there was a lattice mismatch between the substrate and ice. They also noted that multilayer ice growth is not favored, being considerably unstable compared to bulk ice. Our calculations which are applicable to 235 K and relative humidities covering the range important in the atmosphere, are consistent with the picture that emerges from the DFT calculations (presumably 0 K results). Other previous theoretical studies of ice nucleation on surfaces have shown the possible formation of hexagonal rings. For example, Shevkunov [25, 26] studied the structure of water on the surface of a β-AgI crystal as well as in microcracks using a Monte Carlo approach. In his simulations, the water molecules organized themselves into monolayer islands on a single surface or into multilayered films in a microcrack with clear features of hexagonal ice. Similar patterns were obtained by Taylor and Hale [27] also on a AgI substrate, although, they found that the ice-like character did not extend beyond the first water layer. Using molecular dynamics to study different modes of water aggregation on  24  2.3. Results and Discussion  o  ao~ 5.15 A  _  = OH = water oxygen  (a)  (b)  Figure 2.4: (a) Top view of a monolayer (112 water molecules) on the Al-surface obtained at 235 K with the TIP5P-E model. (b) 2-D projection in the x-y plane of the yellow hexagonal pattern overlayed in (a) and the underlying hydroxyl groups forming the Alsurface. The ao lattice parameter of the rings matches that of the lattice (5.1535) but not ice Ih (4.52). The water oxygens and hydrogens are blue and white, respectively. The lattice oxygen, hydrogen, silicon and aluminum atoms are red, white, brown and grey, respectively.  25  2.4. Conclusion the surface of CaF2 and BaF2 , Wassermann et al. [28] found that at high coverage an ice Ih layer only formed on the BaF2 substrate and not on the CaF2 substrate. Hexagonal patterns have also been observed experimentally, and studied with DFT calculations for water on metal surfaces [29, 30]. Clearly, water structure on a surface is very sensitive to the underlying substrate, which is consistent with the results presented here.  2.4  Conclusion  Earlier attempts to understand water adsorption and ice nucleation on kaolinite have focussed almost exclusively on the Al- and Si-surfaces. Our calculations demonstrate that edges can also make an important contribution to water uptake, and might well play an important role in heterogeneous surface chemistry that involves water. As expected, the Si-surface is hydrophobic and exhibits no tendency to adsorb water at pressures below saturation. On the Al-surface the adsorption displays first-order characteristics, with evidence of collective behavior provided by the hysteresis observed for adsorption-desorption simulations. Adsorption on the edges is much more continuous, and appears to be dominated by the very strong water- surface interactions. Despite some differences in adsorption both SPC/E and TIP5P-E give similar monolayer structures. Detailed examination of the structure of the water film on the Al-surface reveals some tendency to form hexagonal patterns, but the observed ao lattice parameter is significantly different from that of ice Ih. In short we did not observe water structures that closely match the structure of hexagonal ice on any of the surfaces investigated. Hence, our calculations do not confirm previous speculations that kaolinite is a good ice nucleus in part because the crystallographic properties of kaolinite promotes water structures on the surface that closely resemble hexagonal ice. If these structures do form, they are extremely rare and are not the thermodynamically preferred state at the surface. Given this, one may suspect that “active sites” such as cracks, pores, steps, etc., are important for ice nucleation on kaolinite. We 26  2.4. Conclusion are currently investigating this possibility (see Chapters 4 and 5).  27  Bibliography [1] Forster, P.; Ramaswamy, V.; Artaxo, P.; Berntsen, T.; Betts,R.; Fahey, D. W.; Haywood, J.; Lean, J.; Lowe, D. C.; Myhre, G.; Nganga, J.; Prinn, R.; Raga, G.; Schulz, M.; Dorland, R. V. Changes in Atmospheric Constituents and in Radiative Forcing., in Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change 2007, edited by S. Solomon, D. Qin, M. Manning, Z. Chen, M. Marquis, K.B. Averyt, M. Tignor, and H.L. Miller, pp. 129-234, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA. [2] Pruppacher, H. R.; Klett, J. D. Microphysics of Clouds and Precipitation, edited by Kluwer Academic, Dordrecht, 2nd rev. ed., 1997, p. 329-331. [3] Cantrell, W.; Heymsfield A. J. Bull. Amer. Meteor. Soc. 2005, 86(6), 795. DeMott, P. J.; Cziczo, D. J.; Prenni, A. J.; Murphy, D. M.; Kreidenweis, S. M.; Thomson, D. S.; Borys, R.; Rogers, D. C. Proc. Nat. Acad. Sci. 2003b, 100, 25, 14655. [4] Glaccum, R. A.; Prospero, J. M. Mar. Geol. 1980, 37(3-4), 295. [5] Zimmerman, F.; Ebert, M.; Worringen, A.; Schulz, L.; Weinbruch, S. Atmos. Environ. 2007, 41, 8219. [6] Zuberi, B.; Bertram, A.K.; Cassa, C. A.; Molina, L. T.; Molina, M. J. Geophys. Res. Lett. 2002, 29, 142. [7] Dymarska, M; Murray, B. J.; Sun, L; Eastwood, M.; Knopf, D. A.; Bertram, A. K. J. Geophys. Res. 2006, 111, D04204, doi:10.1029/2005JD006627. [8] Tunega, D.; Gerzabek, M. H.; Lischka, H. J. Phys. Chem. B 2004, 108, 5930. [9] Wang, Y.-H.; Siu, W.-K. Can. Geotech. J. 2006, 43, 587. [10] Hu, X. L; Michaelides, A. Surf. Sci. 2007, 601, 5378. [11] Hu, X. L; Michaelides, A. Surf. Sci. 2008, 602, 960. [12] Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. [13] Lisal, M.; Nezbeda I.; Smith, W. R. J. Phys. Chem. B 2004, 108, 7412. [14] Cygan, R. T.; Liang, J.-J.; Kalinichev, A. G. J. Phys. Chem. B 2004, 108, 1255. [15] Vega, C.; Sanz, E.; Abascal, J. L. F. J. Chem. Phys. 2005, 122, 114507. 28  Chapter 2. Bibliography [16] Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910. [17] Bish, D. L. Clays Clay Miner. 1993, 41, 783. [18] Shelley, J. C.; Patey, G. N. Mol. Phys. 1996, 88, 385. [19] Yeh, I.-C.; Berkowitz, M. L. J. Chem. Phys. 1999, 111, 3155. [20] Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids (Clarendon, Oxford, 1989). [21] Delville, A. J. Phys. Chem. 1995, 99, 2033. [22] Riemer, N.; Vogel, H.; Vogel, B.; Schell, B.; Ackermann, I.; Kessler, C.; Hass, H. J. Geophys. Res. 2003, 108(D4), 4144. [23] Dentener, F. J.; Crutzen, P. J. J. Geophys. Res. 1993, 98, 7149. [24] Turnbull, D.; Vonnegut, B. Industr. Eng. Chem. 1952, 44, 1292. [25] Shevkunov, S. V. Colloidal Journal 2005, 67, 497. [26] Shevkunov, S. V. Colloidal Journal 2007, 69, 360. [27] Taylor, J. H.; Hale, B. N. Phys. Rev. B 1993, 47, 9732. [28] Wassermann, B.; Reif, J.; Matthias, E. Phys. Rev. B 1994, 50, 2593. [29] Cerd´a, J.; Michaelides, A.; Bocquet, M.-L.; Feibleman, P. J.; Mitsui, T.; Rose, M.; Fomin, E.; Salmeron, M. Phys. Rev. Lett., 2004, 93, 116101. [30] Michaelides, A.; Morgenstern, K. Nature Mater., 2007, 6, 597.  29  Chapter 3 Simulation of Water Adsorption on Kaolinite under Atmospheric Conditions∗ 3.1  Introduction  Mineral dust particles are abundant in the atmosphere, with estimated total global emissions of 800-1500 Tg/year [1–5]. The components of aerosolized mineral dust found in the atmosphere include illite, kaolinite, montmorillonite, quartz and calcite. Kaolinite represents a significant component, comprising approximately 5-10% of aerosolized mineral dust [6]. In the atmosphere mineral dust particles can take up water and act as cloud condensation nuclei or ice nuclei [7–13]. Thus, mineral dust particles can influence climate by changing the frequency and properties of clouds [8, 10, 14, 15]. To completely understand these processes one must first understand the interaction of water with mineral dust particles. Mineral dust particles can also provide a surface on which heterogeneous reactions can occur in the atmosphere. These reactions could be a sink of gas-phase atmospheric species, and influence the chemical composition of the mineral dust particles [16]. Examples of reactions that may occur on mineral dust particles are the hydrolysis of N2 O5 , and the uptake of HNO3 . Recent studies [17–23] have shown that these reactions can be strongly influenced by the presence of water on mineral surfaces. Hence, in order to understand these reactions under atmospheric conditions, an understanding of the interaction of water with mineral surfaces is also necessary. A version of this chapter has been published. T. Croteau, A. K. Bertram, G. N. Patey (2009), Simulation of Water Adsorption on Kaolinite under Atmospheric Conditions, J. Phys. Chem. A, 113:7826-7833. ∗  30  3.1. Introduction One method of studying the interaction of water with mineral surfaces is to use computational methods. The interaction of water with kaolinite has been investigated in several ab initio studies, which have confirmed the hydrophobic and hydrophilic nature of the Siand Al-surfaces, respectively [24–28]. The hydrophobicity of the Si-surface refers to its rather weak interaction with a water molecule, that, in some instances, can lead to water clustering instead of layering. The recent computational studies of Hu and Michaelides [27, 28] are of particular relevance for ice nucleation studies. These authors have examined water on the Al-surface of kaolinite employing density functional theory (DFT). They found that a two-dimensional (2D) water layer with a stability comparable to that of ice Ih could be formed on the Alsurface, but they note a mismatch between the ice Ih lattice and the 2D water structure imposed by the substrate. Moreover, they note that multilayers of water on the Al-surface are less stable than bulk ice. More recently, we used grand canonical Monte Carlo (GCMC) calculations to study water adsorption and structure on kaolinite surfaces as a function of relative humidity (RH) at 235 K, with a focus on ice nucleation [29]. We found monolayer formation but the lattice parameters do not match that of ice Ih. Indeed, the strain produced on an ice embryo by the lattice mismatch strongly suggests that an atomically smooth, defect-free Al-surface would not be a good ice nucleus, contrary to some previous speculation [30]. There have been a number of simulation studies of water on various related surfaces, including the Al- and Si-surfaces of kaolinite [31–33], CaF2 and BaF2 [34], β − AgI [35–37], metals [38, 39], and model hexagonal surfaces [40]. However, we note that the vast majority of this work has employed closed systems, and does not address adsorption questions. A relevant exception is the work of Delville [33], who focused on the Si-surface. The present paper describes a significant extension of our initial studies of water adsorption on kaolinite at 235 K, which focused on ice nucleation [29]. Here we investigate in detail water adsorption on kaolinite at both 298 and 235 K, again employing the GCMC  31  3.1. Introduction method, which is well suited to address adsorption questions. The calculations at 298 K allow us to compare qualitatively with recent adsorption experiments [41], and those at 235 K are motivated by ice nucleation experiments carried out at that temperature [9, 10]. Four kaolinite surfaces are considered. These include the hydrophilic Al-surface, the hydrophobic Si-surface, and two so-called edge configurations (defined below). Some of the topics addressed in the present paper that were not included in our preliminary report are as follows. We examine the modeling/simulation techniques employed, and identify potential problems that can arise. For example, we note that arbitrary truncation of a kaolinite slab (a lattice of point charges) can lead to strong unphysical fields at long distances from the surface of interest. These fields can seriously distort the results obtained, but we show that such effects can be essentially eliminated by carrying out simulations with two slabs arranged such that the long-range, truncation-dependent fields cancel exactly. We also examine the convergence of the Ewald sums necessary to take account of the long-range Coulombic interactions. We determine water adsorption isotherms at 298 and 235 K for all four surfaces, and the results are compared with recent laboratory studies. We study the orientation and binding energies of single adsorbed water molecules at 10 K to understand the differences in water uptake by the four kaolinite surfaces. This includes determination of the binding energies and geometric arrangements of single water molecules on the surface, as well as water-surface hydrogen bond analyses. Water structure in the monolayers is analyzed using density and hydrogen-bond number profiles. The remainder of this paper is organized in three parts. The models and simulation issues are discussed in Section II, the results are presented and discussed in Section III, and our conclusions are summarized in Section IV.  32  3.2. Model and Method  3.2  Model and Method  The model used for the clay lattice is CLAYFF, a general force field developed by Cygan et al. [42]. In this model the short-range interactions of all atoms are represented by Lennard-Jones (LJ) potentials. For the lattice oxygens, the parameters assigned are identical to those of the single point charge (SPC) model for water [43]. Partial charges for all lattice atoms were derived by periodic density functional theory calculations based on the electronic structures of simple oxides and hydroxides (e.g., quartz (α-SiO2 ), corundum (α-Al2 O3 ), boehmite (γ-AlO(OH)), and others) [42]. The structure of kaolinite obtained by Bish [44] was used to improve the optimization of the interaction parameters. The versatility of this model, mainly due to the use of a nonbonded description for the interatomic interactions of the lattice, gives a general force field suitable for application to a wide variety of different clays. To keep the computational requirements practical, we use a rigid lattice corresponding to the most stable kaolinite configuration [44]. We employ the extended single point charge SPC/E [45] and TIP5P-E [46] models to represent the water-water and water-lattice interactions. The SPC/E is a widely used three-site model with the Coulombic interactions described by three point charges (one for each atomic nucleus) embedded in a LJ sphere centered at the oxygen. TIP5P-E is a recent reparameterization of the original TIP5P model [47] to better treat the long-range interactions when using Ewald sums. It is a five-site model with positive charges located at the hydrogen nuclei and negative charges at “lone pair” positions, all embedded in a LJ sphere centered at the oxygen nucleus. The water models are rigid with oxygenhydrogen bond lengths of 1 ˚ A and 0.9572 ˚ A, and bond angles of 109.47o and 104.52o , for the SPC/E and TIP5P-E models, respectively. For TIP5P-E, the oxygen “lone pairs” are located at 0.7 ˚ A, with a “bond angle” of 109.47o . Both models give a reasonably good description of water structure and properties [45, 48] under ambient conditions. We note that the TIP5P-E model gives a more accurate prediction of the melting point of ice, 274 K compared with 215 K for the SPC/E model [49]. However, this difference does not 33  3.2. Model and Method appear to have a significant influence on water adsorption on kaolinite. Two different water models are used to ensure as far as possible that model-dependent effects are not significantly influencing our observations and conclusions. It should be noted that both water models employed are rigid point-charge models that do not allow for water dissociation. While water dissociation might play a significant role in surface chemistry, we would not expect it to have a large influence on the adsorption isotherms. The configurational energy for our system, U, is given by the sum of Coulombic and LJ interactions e2 U= 4πǫo  i=j  qi qj + rij  Do,ij i=j  Ro,ij rij  12  Ro,ij −2 rij  6  ,  (3.1)  where the sums are over all lattice and water interaction sites, and rij is the distance between sites i and j. The first term represents the Coulombic interactions, where qi is the point charge on site i, e is an elementary charge, and ǫo is the dielectric permittivity of free space. The second term represents the LJ contribution and the energy and distance parameters for the different site-site interactions, Do,ij and Ro,ij , respectively, are given in Ref. [42]. For the cross interactions, the Lorentz-Berthelot rules, Do,12 =  Do,11 Do,22  and Ro,12 = (Ro,11 + Ro,22 )/2, are used. Kaolinite [Al2 Si2 O5 (OH)4 ] is a clay mineral that has a layered structure. These layers consists of octahedral aluminum (Al-surface) and tetrahedral silicon (Si-surface). Kaolinite consists of many of these layers attached together by hydrogen bonds between the hydroxyl groups extending from the Al-surface and bridging oxygens on the Si-surface. The large number of hydrogen bonds connecting the layers makes them almost inseparable. Thus kaolinite is described as a non-expanding or non-swelling clay. The threedimensional simulation cell illustrated in Fig. 2.1 is constructed such that the Al- and Si-surfaces are parallel to the (001) plane. The rectangular lattices periodically replicated in the xy plane are composed of 48 kaolinite unit cells for a total of 816 atoms per surface. The x and y dimensions of the cell are 30.921 ˚ A and 35.7676 ˚ A, respectively. 34  3.2. Model and Method In addition to the Al- and Si-surfaces, we consider two additional surfaces that commonly occur in kaolinite samples [50–52] and are sometimes called “edges”, a term that we will adopt here. Two extreme edge cases are considered, one which has “bare” oxygens exposed (unprotonated edge) and another where all oxygens are protonated (protonated edge). The unprotonated edge is constructed by cleaving the Al-surface along the (100) plane such as to have oxygen atoms exposed (see Fig. 2.1). The surface is then rotated 90o counterclockwise to bring the exposed oxygens into the (001) plane, and expanded to have x and y dimensions of 29.524 ˚ A and 35.76764 ˚ A, respectively, and contain 816 atoms. A protonated edge is constructed by placing hydrogen atoms (a total of 64), each with charge 0.4250e, 1 ˚ A above every bare oxygen atom. To counterbalance this extra charge and maintain a neutral surface, an opposing negative charge was equally distributed on all sites below the first layer of Al and Si atoms. Specifically, a negative charge of 0.05e is placed on all 544 sites below the surface. This procedure is obviously somewhat arbitrary, but changing the surface thickness (thus reducing the charge per subsurface site) did not have any significant effect on the results obtained. We remark that many “edge” surfaces could be constructed by cleaving kaolinite along different planes. The particular edges considered here have been discussed by previous authors [50–52] and are believed to make a significant contribution to the surface area of kaolinite samples. Care must be taken in simulating slabs that are periodically infinite in two dimensions (x and y), but arbitrarily truncated in the third (z). Although our samples are overall electrically neutral, arbitrary truncation of the discrete point charge lattice can induce local charge imbalance at the resulting surface, and significant net polarization. This can create unphysical fields that depend strongly on the slab thickness, and can act at large distances from the surface. Such effects would not occur for an infinitely thick sample because charges in neighboring layers preserve local charge neutrality. As an illustration of the problem, the interaction energy of a water molecule (arbitrary orientation) with a single surface is shown in Fig. 3.1 (top). We note that the interaction energy does  35  3.2. Model and Method not decay quickly and is about −3 kJ/mol at a distance of 25 ˚ A from the surface. This unphysical interaction strongly depends on the slab thickness, is different for different surfaces, and can significantly influence the results obtained. For example, for the Aland Si- surfaces one can observe small isolated water “towers” on the surface, such as those previously reported by Delville [33] for the Si-surface. We believe that such effects are artefacts of the long-range, unphysical field resulting from truncation of the charge lattice. One way to solve this problem is to add an opposing surface having the same size and layer structure (see Fig. 2.1), such that the unphysical fields produced by truncation simply cancel. We note that this technique was previously employed by Warne and Cosgrove [31] for similar reasons. Results for a single water molecule in a two-slab arrangement are also shown in Fig. 3.1, and we see that the interaction now approaches zero at about 10 ˚ A from the surface. The strong dependence on slab thickness (not shown) found for a single slab vanishes when two slabs are used. Furthermore, unusual structures such as isolated water towers without any obvious physical explanation no longer occur. The configurational energies were obtained using a method previously applied to slab geometries with long-range electrostatic interactions [53, 54]. We wish to calculate the long-range interactions in a slab of finite thickness. One way to do this is to periodically repeat the cell shown in Fig. 2.1 in three dimensions and apply the usual 3D Ewald method [1], but leaving enough empty space to prevent the undesired images in the z direction from having a significant physical influence [53, 54]. For our systems an empty gap of 107 ˚ A was found to be sufficient for this purpose, increasing the gap further had no noticeable effect. Comparing adsorption results, it was judged sufficient to employ 6858 wave vectors in the Ewald sums, as no significant change was observed with larger ˚−1 and a reciprocal numbers. The other Ewald parameters employed were α = 0.175 A space cut-off of 2.0 ˚ A−1 . The LJ and real space interactions were cut off at half the length of the x side of the simulation cell.  36  3.2. Model and Method  -2  Energy (kJ/mol)  -2.5 -3  -3.5 -4 0  -0.5 -1 -1.5 -2  0  5  10  15  20  25  z (Å) Figure 3.1: The surface-water interaction energy of a single, randomly oriented, SPC/E water molecule, as a function of the perpendicular distance from the center of the highest hydroxyl hydrogen atom on the Al-surface. The top and bottom panels show the results obtained for simulation cells with one and two slabs, respectively. The blue (circles), red (squares), and green (triangles) lines correspond to 6858, 15624, and 226980 k-space wave vectors, respectively.  37  3.2. Model and Method The simulations were carried out employing the GCMC method [1]. This method allows the equilibrium properties of an open system to be obtained at fixed chemical potential, µ, volume and temperature. In the present calculations the number of water molecules, NW , fluctuates. A Monte Carlo step is defined as an attempted insertion, deletion, or displacement (rotation or translation) of a water molecule, and all three possibilities were attempted with equal probability. The acceptance rate for insertion (deletion) moves was ∼ 0.01%, which is comparable to that found in GCMC simulations of bulk water [56]; the maximum size of the displacement moves was adjusted to give an acceptance rate of ∼ 50%. To ensure convergence, the systems were equilibrated for at least 5 × 107 Monte Carlo (MC) steps. Following equilibration, averages were collected for 2×108 or more MC steps and include, energy, water content, density profile, and hydrogenbond number. This procedure was repeated for a wide range of chemical potentials and at different temperatures to produce adsorption isotherms and related properties for all four surfaces. For the two-slab geometry discussed above (Fig. 2.1), the number of layers included in each kaolinite slab did not significantly influence the results, and a single layer (constructed as described above) was sufficient for this investigation. However, since we are interested in modeling adsorption on single surfaces rather than in slits of finite width, care must be taken to place the slabs sufficiently far apart such that the water adsorption on each slab is not influenced by the presence of the other. Most simulations were carried out with the slabs separated by 30 ˚ A, but some test runs were performed with a separation of 60 ˚ A. These confirmed that, at least before the filling transition at saturation (bulk condensation of the vapor), the surfaces adsorbed independently. The results reported were obtained by averaging values for both surfaces, which improves the statistics. Some simulations were also performed varying the the x and y cell dimensions by factors of 0.5 and 1.5, and again these tests revealed no significant system-size dependence.  38  3.3. Results and Discussion  3.3 3.3.1  Results and Discussion Water Adsorption at 298 and 235 K  Adsorption isotherms for all four surfaces obtained at 298 K are displayed in Fig. 3.2 as a function of chemical potential. The water coverage plotted on the vertical axis is the average number of water molecules adsorbed on a single surface. For atmospheric comparisons it is often useful to convert the chemical potential scale into pressure/saturated vapor pressure (P/Po). Since the saturated vapor pressure of water at 298 K and 235 K is low (< 0.1 atm), the vapor pressures can be estimated using the ideal gas relationship, P = eβµ /βΛ3 , where Λ is the thermal de Broglie wavelength. The saturated vapor pressure Po for our models can be found by setting P/Po = 1 at the filling transitions observed in the simulations. Some P/Po estimates are given below. In contrast to the Al-surface and the edges, for both water models considered, the Si-surface remains completely dry up to the filling transition. For this reason no further adsorption results are presented for the Si-surface. Qualitatively similar adsorption behavior is observed for both water models. Both edges adsorb readily even at low chemical potentials (vapor pressures), the protonated edge having the greatest affinity for water. Adsorption begins at very low chemical potentials for these surfaces, i.e., below µ ≈ −75.0 kJ/mol (P/Po ≈ 9.6 × 10−4 ) and −90.0 kJ/mol (P/Po ≈ 2.3 × 10−6) for the unprotonated and protonated edges, respectively. As the chemical potential is increased, the water uptake increases slowly, partially covering the surface at well defined binding sites (see discussion below) until a monolayer is obtained at µ ≈ −63.0 kJ/mol (P/Po ≈ 0.12)(SPC/E). The film then thickens and, finally, a transition to a completely filled state occurs. The strong water affinity of the edges demonstrates that they can make an important contribution to water uptake, and could well play an important role in heterogeneous surface chemistry that involves water. In contrast, adsorption on the Al-surface occurs over a much narrower range of chem-  39  3.3. Results and Discussion  600  500  <Nw>  400  300  200  100  0 -96 -90 -84 -78 -72 -66 -60 -54 -90 -84 -78 -72 -66 -60 -54 -48  µ (kJ/mol)  Figure 3.2: Adsorption isotherms at 298 K for the four surfaces considered, obtained with a slab separation of 30 ˚ A. Results for the SPC/E (left) and TIP5P-E (right) models are shown. The black (circles), red (squares), green (diamonds), and blue (triangles) curves correspond to the unprotonated edge, protonated edge, Al-surface, and Si-surface, respectively.  40  3.3. Results and Discussion ical potentials (which is slightly larger for TIP5P-E) suggesting a weaker water affinity on that surface. On the chemical potential scale shown in Fig. 3.2, it is only just prior to filling, at µ ≈ −59.2 kJ/mol (SPC/E) and µ ≈ −56.5 kJ/mol (TIP5P-E), that monolayer coverage is achieved. Note, that we take monolayer coverage to be a density of 1 water molecule per 9 ˚ A2 of surface area (assuming a smooth surface), or about 100 molecules on a single surface in our simulation cell. In terms of the relative vapor pressure, P/Po , the adsorption range on the Al-surface is of atmospheric relevance. The P/Po range over which adsorption occurs is approximately 0.55 − 1 at 298 K, and 0.2 − 1 at 235 K, which covers the range of interest for atmospheric science [14]. For all surfaces considered, adsorption isotherms at 235 K exhibit trends similar to those described above for 298 K. A comparison is made in Fig. 3.3, where results obtained at 298 and 235 K are plotted for both edges and the Al-surface. For all three surfaces, the overall shape of the adsorption isotherms is similar at both temperatures. As noted above, the Al-surface takes up water just before filling occurs, and both edges adsorb water readily at low chemical potentials. The thickness achieved by the water films before filling does show some variation with temperature, particularly for the SPC/E model. At 235 K, the edges acquire slightly thicker films before the filling transition occurs, at µ ≈ −54.0 kJ/mol for SPC/E and µ ≈ −49.0 kJ/mol for TIP5P-E. On the Al-surface, there is a significant difference between the SPC/E and TIP5P-E models. The TIP5P-E model shows no increase in film thickness before filling when the temperature is reduced from 298 to 235 K. However, the SPC/E model exhibits an increase in thickness, forming essentially a bilayer just prior to filling. To check that this increased thickness is not associated with the two-slab geometry, simulations were performed with the Al-surfaces separated by 60 ˚ A. The results obtained are also shown in Fig. 3.3 (bottom panel) and verify that the additional thickening observed at 235 K is not related to interactions with the second surface. Some associated desorption isotherms for the SPC/E model are also plotted in Fig.  41  3.3. Results and Discussion 3.3. These were obtained by carrying out simulations at decreasing chemical potentials, starting with an initial configuration obtained at a chemical potential that is just below the filling transition on the adsorption curve. Simulations were carried out at decreasing chemical potentials until all water molecules evaporate. We noted above that adsorption on the edges is relatively smooth before the filling transition. This, together with the fact that the desorption curves reveal relatively little hysteresis, suggests that water adsorption on the edges is essentially a continuous process, dominated by water-surface interactions. In contrast, water adsorption on the Al-surface is not smooth, and “jumps” rather sharply from essentially nothing to monolayer coverage. Additionally, the desorption curves show a large hysteresis at both temperatures. These observations indicate that water adsorption on the Al-surface has strong “first-order” characteristics, with important collective behavior amongst the water molecules. Recently Schuttlefeld et al. [41], investigated water uptake on kaolinite minerals using ATR-FTIR spectroscopy coupled with quartz crystal microbalance measurements at 296 K. The data shows that kaolinite particles take up water continuously from 0 to 90% RH (P/Po × 100). At an RH of 50% they observed approximately 6-25 monolayers of water uptake and at 80% RH they observed approximately 10-40 monolayers of water uptake, depending on the source of the kaolinite. This is in contrast to our simulations at 298 K, where we do not observe more than approximately monolayer coverage below 100% RH for any of the surfaces considered. The simulations suggest that the water uptake observed in the laboratory experiments cannot be explained by water uptake on defect-free surfaces.  3.3.2  Orientations and Binding Energies of Singly Adsorbed Water Molecules  The differences in water uptake by the four kaolinite surfaces can be largely understood by considering the low-temperature, water-surface binding energies. Results obtained in NV T Monte Carlo simulations at 10 K involving only a single water molecule are 42  3.3. Results and Discussion  250 200 150 100 50  <Nw>  200 150 100 50  300 200 100 0  -62  -60  -58  -56  -54  -58  µ (kJ/mol)  -56  -54  -52  -50  -48  Figure 3.3: Adsorption isotherms obtained for the unprotonated edge (top panel), the protonated edge (middle panel), and the Al-surface (bottom panel) using the SPC/E (left) and TIP5P-E (right) models, at 298 K (black curves with filled circles) and 235 K (red curves with filled squares). The green curve with empty squares was obtained at 235 K with a slab separation of 60 ˚ A. All other curves were obtained with a slab separation ˚ of 30 A. The blue and yellow dashed curves are desorption isotherms.  43  3.3. Results and Discussion Table 3.1: Binding energies and hydrogen-bond numbers for single water molecules obtained using NV T Monte Carlo simulations at 10 K. The numbers in brackets represented one standard deviation. Surface Type Water-Surface < Uwater-surface > Hydrogen Bonds (kJ/mol) Water model SPC/E TIP5P-E SPC/E TIP5P-E Si-surface 1 0 -21.64(0.03) -21.03(0.04) Al-surface 3 1 -46.36(0.04) -45.64(0.02) unprotonated edge 1 1 -73.45(0.02) -71.42(0.03) protonated edge 2 2 -94.14(0.02) -78.51(0.04) summarized in Table 3.1. The standard deviations included in Table 3.1 were obtained by dividing Monte Carlo runs of 106 steps into ten equal blocks, and assuming that the block averages are independent estimates of the water-surface interaction energy. Other error estimates using widely spaced single configuration energies gave similar values. We note that the average water-surface energies obtained for the SPC/E and TIP5P-E models are in good agreement except for the protonated edge, where the SPC/E interaction is more attractive by ∼ 15.6 kJ/mol. For both models the single water molecule binding energies correlate well with the water adsorption curves shown in Figs. 3.2 and 3.3. The weakest interaction is obtained for the Si-surface which does not adsorb water before the chemical potential reaches bulk saturation, and the strongest interactions are for the edges which adsorb at very low chemical potentials (vapor pressures). In this case, the very strong water-surface interactions can support very low coverages both at 298 and 235 K. For the Al-surface, which adsorbs near bulk saturation but in the atmospherically relevant range, the interaction energy lies between the two extremes. For the Al-surface, the direct water-surface interaction is not strong enough for significant sub-monolayer coverage, and collective water-water interactions are necessary to stabilize the monolayer. This explains the relative sharpness of the transitions on this surface (Figs. 3.2 and 3.3), and the hysteresis noted above. It is instructive to consider the surface-single-water-molecule interactions in more detail. The geometric arrangements obtained at 10 K are shown in Fig. 3.4, indicating  44  3.3. Results and Discussion clearly the clay atoms that form hydrogen bonds with the water molecule. We define a hydrogen bond using geometric criteria similar to those applied in earlier work [57, 58]. Specifically, a water-water hydrogen bond is said to exist if and only if the water-oxygenwater-oxygen distance is less than 3.41 ˚ A, the water-hydrogen-water-oxygen distance is less than 2.38 ˚ A, and the O-O-H angle is less than 35o . This definition can be generalized to water-surface hydrogen bonds in an obvious manner. The lengths and angle chosen are clearly arbitrary to some extent, and using different parameters might in some cases influence the number of hydrogen bonds counted. In the present case, the number of water-surface hydrogen bonds is somewhat sensitive to the angle used, but 35o appears to be a reasonable choice yielding hydrogen-bond numbers consistent with those previously reported for the Al-surface (see below). We first consider the Si-surface which has the weakest interaction with water. We note that the energies obtained for both water models are very similar despite the differing orientations [Figs. 3.4(a) and3.4(b)]. For the SPC/E model [Fig.3.4(a)] there is a strong hydrogen bond between the water hydrogen and the closest bridging oxygen of the surface. The water hydrogens also have strong attractive interactions with other neighboring bridging oxygens. However, the structure shown in Fig. 3.4(a) leads to highly repulsive interactions between the water hydrogens and silicon atoms of the surface, leading to a net attraction that is relatively weak. For the TIP5P-E model the most attractive and repulsive interactions are modified, and come from the lone-pair charges interacting with surface silicon and oxygen atoms, respectively. This gives the structure shown in Fig. 3.4(b), where both water hydrogens lie flat in the xy plane, and there are no hydrogen bonds with the surface. As noted above, somewhat surprisingly given the structural differences, both water models give essentially the same binding energy. A water molecule binds much more strongly to the Al-surface, and the energies we obtain (Table 3.1) compare well with previous DFT calculations [28]. Note that when the zero point energy (not present in our classical model) is removed, the DFT binding energy  45  3.3. Results and Discussion  (a)  (b)  (c)  (d)  (e)  (f )  (g)  (h)  Figure 3.4: Snap shots of a water molecule at 10 K for the Si-surface (a,b), Al-surface (c,d), unprotonated edge (e,f), and protonated edge (g,h) using the SPC/E (a,c,e,g) and TIP5P-E (b,d,f,h) water models. The water oxygens are blue, lattice oxygens are red, hydrogens are white, aluminums are grey, silicons are brown, lone-pairs are orange, and the surface atoms involved in hydrogen bonding with the water molecule are indicated in green.  46  3.3. Results and Discussion is ∼ −43.4 kJ/mol, which is close to our values (−46.4 kJ/mol for SPC/E and −45.6 kJ/mol for TIP5P-E). The SPC/E model promotes a very attractive interaction between a downward pointing water hydrogen and a hydroxyl oxygen forming a donor hydrogen bond [Fig. 3.4(c)]. Two additional acceptor hydrogen bonds are formed between the water oxygen and two surface hydroxyl hydrogens giving a total of three. The number of hydrogen bonds and the structure shown in Fig. 3.4(c) is in agreement with the DFT result reported by Hu and Michaelides [28]. The TIP5P-E model gives a somewhat different structure with the strongest attractive interactions again coming through the lone-pair charges. One lone pair forms an acceptor hydrogen bond with a hydroxyl hydrogen, and interacts favorably with neighboring underlying aluminum atoms. This gives a structure with the water hydrogens lying flat in the xy plane as shown in Fig. 3.4(d). We note that, despite the orientational differences, both models give similar net binding energies. Also, the monolayer structures obtained at higher temperatures, which are clearly influenced by water-water as well as water-surface interactions, are very similar [see Fig. 3.5(d)]. Both low-temperature structures [Figs. 3.4(c) and 3.4(d)] actually occur in the monolayers obtained at 298 and 235 K, for both water models. However, we note that it is only for the TIP5P-E model that the water-surface interaction is attractive for both configurations, likely explaining the slightly larger adsorption range observed for TIP5P-E. For both edges, the low energy water structures occur near the junction between the Si- and Al-surfaces [Figs. 3.4(e) and 3.4(f)]. For the unprotonated edge, a very favorable interaction is formed between the water oxygen (SPC/E) or lone-pairs (TIP5P-E) and a hydroxyl hydrogen from the Al-surface, thus forming an acceptor hydrogen bond with both models. Note that, although there are a large number of bare oxygen atoms exposed on the unprotonated edge, no donor hydrogen bonds are observed. For the protonated edge, the most attractive interactions are formed between the water hydrogens and hydroxyl and bridging oxygens, leading to two donor hydrogen bonds for both models [Figs. 3.4(g) and 3.4(h)]. For both types of edge, the delocalization of the  47  3.3. Results and Discussion partial negative charge from the center of the water oxygen nucleus (SPC/E) to the lonepairs locations (TIP5P-E) has minimal impact on the structures obtained. However, it does give a significant change in the binding energy on the protonated edge, increasing from −94.1 kJ/mol (SPC/E) to −78.5 kJ/mol (TIP5P-E), whereas the energies obtained for the unprotonated edge are very similar. For the protonated edge the difference is explained by the fact that for the TIP5P-E model the water-surface interaction is more sensitive to the presence of neighboring, non-hydrogen-bonding surface atoms.  3.3.3  Water Structure at Monolayer and Sub-Monolayer Coverage  Snapshots obtained for sub-monolayer coverage on both edges, and for monolayer coverage on the unprotonated edge and the Al-surface are shown in Fig. 3.5. These results were obtained for the TIP5P-E model at 235 K. Images obtained for the SPC/E model at 235 K, and for both models at 298 K, are very similar to those given in Fig. 3.5, and hence are not shown. It is evident from Figs. 3.5(a) and 3.5(b), that, at low coverage, water molecules adsorb on the edges at well-defined binding sites and with particular orientations (as described above). We note that for both edges and the Al-surface, at monolayer coverage, the adsorbed water layer exhibits little order and no long-range patterns are discernable [Figs. 3.5(c) and 3.5(d)]. A detailed analysis of the nature of the monolayer on the Al-surface, and its possible relevance to ice nucleation under atmospheric conditions, was given in a previous communication [29]. To briefly summarize, although some hexagonal structures do form on the Al-surface [Fig. 3.5(d)], the relevant lattice parameter is considerably larger than that required to match ice Ih. Structure factors calculated for the surface layer employing reciprocal lattice vectors appropriate for hexagonal ice [59, 60] are close to zero. The expected strain on the ice embryo due to this mismatch would depress the nucleation temperature to well below the atmospherically relevant range [29, 61]. Therefore, our 48  3.3. Results and Discussion  (a)  (b)  (c)  (d)  Figure 3.5: Snap shots obtained at 235 K using the TIP5P-E model for low water coverages on the unprotonated edge (a) and the protonated edge (b), and for monolayer coverages on the unprotonated edge (c) and the Al-surface (d). The water oxygens and hydrogens are blue and white, respectively.  49  3.3. Results and Discussion results do not support the theory that kaolinite serves as a good ice nucleus because the crystallographic properties of kaolinite are such that the Al-surface is particularly hospitable to water structures that closely match hexagonal ice.  3.3.4  Water Density and Hydrogen Bonding at Monolayer Coverage  Water density profiles, ρ(z), (in g/cm3 ) and hydrogen-bond numbers (as functions of z) for monolayer coverage on the Al-surface and both edges are plotted in Fig. 3.6. In this figure z = 0 corresponds to the center of the highest surface hydroxyl hydrogen for the Al-surface, and to the center of the lowest exposed oxygen atom for the edges. We note that the density is quite high near the surfaces indicating a rather compact layer, and shows some structural features associated with the atomic granularity of the surface. For each surface the average total number of hydrogen bonds per water molecule is shown together with the water-surface and water-water contributions. Generally, we see that near the surface, in the most compact region of the layer, water molecules tend to have about four hydrogen bonds, as one would expect, but the number can be a little higher or lower depending on the surface and the particular location. Also, not unexpectedly, there are fewer hydrogen bonds in the outer fringes of the layer. Comparing the average number of water-surface hydrogen bonds in the monolayer with the numbers given in Table 3.1 for a single water molecule, some differences are noticeable. These differences arise because in the monolayer interactions with neighboring water molecules strongly influence the geometrical arrangements. This tends to reduce the differences between the water models giving qualitatively similar monolayer structures, as noted above. A density and hydrogen bond analysis was also carried out at 298 K. The results obtained were similar to those discussed above, except that the structural features tended to be a little broader and less sharp than those observed at 235 K. This is consistent with the larger fluctuations expected at the higher temperature. 50  H-bond number  ρ(z)  3.3. Results and Discussion  5 4 3 2 1 3 2 1 3 2 1 5 4 3 2 1 0  1  2  3  4  5  6  7 0  1  2  3  4  5  6  7  z (Å) Figure 3.6: The water density (g/cm3 ) profile (top most panel) and the average number of hydrogen bonds per water molecule at 235 K, obtained for monolayer coverage using the SPC/E (left) and TIP5P-E (right) models. The dashed (black), dotted-dashed (red), and solid (green) lines in the density plot correspond to the unprotonated edge, protonated edge, and Al-surface, respectively. The lower three panels give the hydrogen-bond numbers for the Al-surface (bottom), protonated edge (middle), and unprotonated edge (top). In these panels the solid (black), dotted-dashed (blue), and dashed (orange) curves represent water-surface bonds, water-water bonds, and total number, respectively.  51  3.4. Summary and Conclusions  3.4  Summary and Conclusions  GCMC simulations have been employed to examine water uptake and structure on the Al-surface, the Si-surface, and two edge-like surfaces of kaolinite. Simulations were carried out for a range of chemical potentials (vapor pressures) at the atmospherically relevant temperatures, 298 and 235 K. Two water models, the SPC/E and TIP5P-E, were considered. At low temperature (10 K), these models prefer different orientations on some surfaces, but the single-molecule binding energies are in good agreement, and, for the Alsurface, compare well with earlier DFT calculations [28]. At 298 and 235 K, both models gave similar adsorption isotherms and monolayer structures. From a modeling perspective, we found that arbitrary truncation of a kaolinite slab along the axis perpendicular to the surface of interest (the z axis in this paper) can lead to unphysical fields far from the surface. These fields can lead to spurious results that strongly depend on where the slab is truncated. We show that this problem can be overcome by simply including two slabs in the simulation cell, oriented such that the unphysical, long-range fields cancel exactly. With this simulation arrangement, the results obtained become quite insensitive to slab thickness. As previously known [33], the Si-surface is relatively hydrophobic, and at both temperatures considered we find no significant water adsorption before saturation. The Al-surface and both edges take up water and form monolayers before saturation, at atmospherically relevant pressures. The edges have by far the greatest affinity for water, adsorbing at very low chemical potentials, and up to monolayer coverage. This, together with the fact that edge-like surfaces account for a significant fraction of the total surface area of kaolinite particles [50, 51], suggests that edges might play an important role in heterogeneous surface chemistry where water is needed for hydrolysis. The mechanism of water uptake on the edges is distinctly different from that of the Al-surface. Adsorption on the edges appears to be driven mainly by strong water-surface interactions. On the edges, water molecules initially adsorb at well-defined periodic bind52  3.4. Summary and Conclusions ing sites, and the coverage grows essentially continuously with chemical potential until a monolayer is achieved. In contrast, on the Al-surface there is practically no sub-monolayer adsorption, and monolayer formation appears as a fairly sharp transition. This indicates that collective behavior associated with water-water interactions is essential to monolayer formation on the Al-surface. Further evidence for this is given by the large hysteresis observed on water desorption from the Al-surface. Very little hysteresis is observed for the edges. The monolayer structure on the Al-surface and on both edges was analyzed through density and hydrogen-bond number profiles. The density profiles ρ(z) reached a maximum of 3 − 5 g/cm3 , depending on the surface, indicating that the monolayers are rather dense. In the dense regions of the monolayers, a water molecule tends to form about four hydrogen bonds, some with the surface, and some with other water molecules of the monolayer. The details of the density and hydrogen-bond number profiles do show some variation with the water model employed, but, overall, the monolayer structures are very similar. Finally, we note that for all surfaces considered, at 298 K, we observe only monolayer coverage before saturation (bulk condensation of the vapor). This contrasts with recent experimental results for kaolinite particles [41], where much higher coverages are reported for relative humidities well below 100%. Our simulations suggest that the experimental observations cannot be explained with atomistically smooth surfaces, and that surface roughness likely plays a very important role in water adsorption by kaolinite particles (Chapters 4 and 5).  53  Bibliography [1] Bauer, S.E.; Balkanski, Y.; Schulz, M; Hauglustaine, D. A.; Dentener, F. J. Geophys. Res.-Atm. 2004, 109 (D2). [2] Lunt, D. J.; Valdes, P. J. J. Geophys. Res.-Atm. 2002, 107 (D23). [3] Tegen, I.; Harrison, S. P.; Kohfeld, K.; Prentice, I. C.; Coe, M.; Heimann, M. J. Geophys. Res.-Atm. 2002, 107 (D21). [4] Tegen, I.; Miller, R. J. Geophys. Res.-Atm. 1998, 103 (D20), 25975. [5] Woodward, S. J. Geophys. Res.-Atm. 2001, 106 (D16), 18155. [6] Glaccum, R. A.; Prospero, J. M. Marine Geology 1980, 37 (3-4), 295. [7] Cziczo, D. J.; Murphy, D. M.; Hudson, P. K.; Thomson, D. S. J. Geophys. Res.-Atm. 2004 109 (D4). [8] DeMott, P. J.; Cziczo, D. J.; Prenni, A. J.; Murphy, D. M.; Kreidenweis, S. M.; Thomson, D. S.; Borys, R.; Rogers, D. C. Proc. Nat. Acad. Sci. 2003, 100 (25), 14655. [9] Eastwood, M. L.; Cremel, S.; Gehrke, C.; Girard, E.; Bertram, A. K. J. Geophys. Res.-Atm. 2008, 113. [10] Dymarska, M.; Murray, B. J.; Sun, L.; Eastwood, M. L.; Knopf, D. A.; Bertram, A. K. J. Geophys. Res.-Atm. 2006, 111 D4. [11] Zuberi, B.; Bertram, A. K.; Cassa, C. A.; Molina, L. T.; Molina, M. J. Geophys. Res. Let. 2002, 29 (10) [12] Twohy, C. H.; Kreidenweis, S. M.; Eidhammer, T.; Browell, E. V.; Heymsfield, A. J.; Bansemer, A. R.; Anderson, B. E.; Chen, G.; Ismail, S.; DeMott, P. J.; Van den Heever, S. C. Geophys. Res. Let. 2009, 36. [13] Twohy, C. H.; Poellot, M. R. Atm. Chem. Phys. 2005, 5, 2289. [14] Forster, P.; Ramaswamy, V.; Artaxo, P.; Berntsen, T.; Betts, R.; Fahey, D. W.; Haywood, J.; Lean, J.; Lowe, D. C.; Myhre, G.; Nganga, J.; Prinn, R.; Raga, G.; Schulz, M.; Dorland, R. V. Changes in Atmospheric Constituents and in Radiative Forcing., in Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change 2007, edited by S. Solomon, D. Qin, M. Manning, Z. Chen, M. Marquis, 54  Chapter 3. Bibliography K.B. Averyt, M. Tignor, and H.L. Miller, pp. 129-234, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA. [15] Cantrell, W.; Heymsfield A. J. Bull. Amer. Meteor. Soc. 2005, 86(6), 795. [16] Usher, C. R.; Michel, A. E.; Grassian, V. H. Chem. Rev. 2003, 103 (12), 4883. [17] Mogili, P. K.; Kleiber, P. D.; Young, M. A.; Grassian, V. H. Atm. Env. 2006, 40 (38), 7401. [18] Mashburn, C. D.; Frinak, E. K.; Tolbert, M. A. J. Geophys. Res.-Atm. 2006, 111 (D15). [19] Goodman, A. L.; Bernard, E. T.; Grassian, V. H. J. Phys. Chem. A 2001, 105 (26), 6443. [20] Krueger, B. J.; Ross, J. L.; Grassian, V. H. Langmuir 2005, 21 (19), 8793. [21] Hanisch, F.; Crowley, J. N. Physical Chemistry Chemical Physics 2001, 3 (12), 2474. [22] Laskin, A.; Wietsma, T. W.; Krueger, B. J.; Grassian, V. H. J. Geophys. Res.-Atm. 2005, 110 (D10). [23] Santschi, C.; Rossi, M. J. J. Phys. Chem. A 2006, 110 (21), 6789. [24] Tunega, D.; Gerzabek, M. H.; Lischka, H. J. Phys. Chem. B 2004, 108, 5930. [25] Tunega, D.; Benco, L.; Haberhauer, G.; Gerzabek, M. H.; Lischka, H. J. Phys. Chem. B 2002, 106, 11515. [26] Tunega, D.; Haberhauer, G.; Gerzabek, M. H.; Lischka, H. Langmuir 2002, 18, 139. [27] Hu, X. L.; Michaelides, A. Surf. Sci. 2007, 601, 5378. [28] Hu, X. L.; Michaelides, A. Surf. Sci. 2008, 602, 960. [29] Croteau, T.; Bertram, A. K.; Patey, G. N. J. Phys. Chem. A 2008, 12, 10708. [30] Pruppacher, H. R.; Klett, J. D. Microphysics of Clouds and Precipitation, 2nd rev. ed.; Kluwer Academic: Dordrecht, The Netherlands, 1997; pp 329-331. [31] Warne, M. R.; Cosgrove, T. Phys. Chem. Chem. Phys. 2000, 2, 3663. [32] Vasconcelos, I. F.; Bunker, B. A.; Cygan, R. T. J. Phys. Chem. C 2007,111, 6753. [33] Delville, A. J. Phys. Chem. 1995, 99, 2033. [34] Wassermann, B.; Reif, J.; Matthias, E. Phys. Rev. B 1994, 50, 2593. [35] Shevkunov, S. V. Colloidal Journal 2005, 67, 497. [36] Shevkunov, S. V. Colloidal Journal 2007, 69, 360. 55  Chapter 3. Bibliography [37] Taylor, J. H.; Hale, B. N. Phys. Rev. B 1993, 47, 9732. [38] Cerd´a, J.; Michaelides, A.; Bocquet, M.-L.; Feibleman, P. J.; Mitsui, T.; Rose, M.; Fomin, E.; Salmeron, M. Phys. Rev. Lett., 2004, 93, 116101-1. [39] Michaelides, A.; Morgenstern, K. Nature Mater, 2007, 6, 597. [40] Nutt, D. R.; Stone, A. J. Langmuir 2004, 20, 8720. [41] Schuttlefield, J. D.; Cox, D.; Grassian, V. H. J. Geophys. Res. 2007, 112, D21303. [42] Cygan, R. T.; Liang, J.-J.; Kalinichev, A. G. J. Phys. Chem. B 2004, 108, 108, 1255. [43] Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981; pp 331. [44] Bish, D. L. Clays Clay Miner. 1993, 41, 783. [45] Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. [46] Lisal, M.; Nezbeda I.; Smith, W. R. J. Phys. Chem. B 2004, 108, 7412. [47] Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910. [48] van der Spoel, D.; van Maaren, P. J.; Berendsen, H. J. C. J. Chem. Phys. 1998, 108, 10220. [49] Vega, C.; Sanz, E.; Abascal, J. L. J. Chem. Phys. 2005, 122, 114507. [50] Brady, P. V.; Cygan, R. T.; Nagy, K. L. J. Colloid and Interface Science 1996, 183, 356. [51] Zbik, M.; Smart, R. S. C. Clays and Clay Minerals 1998, 46(2), 153. [52] Wang, Y.-S.; Siu, W.-K. Can. Geotech. J. 2006, 43, 587. [53] Shelley, J. C.; Patey, G. N. Mol. Phys. 1996, 88, 385. [54] Yeh, I.-C.; Berkowitz, M. L. J. Chem. Phys. 1999, 111, 3155. [1] Allen, M. P.; Tildesley, D. J. Computer simulation of liquids, (Oxford, 1987). [56] Shelley, J. C.; Patey, G. N. J. Chem. Phys. 1995, 102, 7656. [57] Luzar, A.; Chandler, D. Nature 1996, 53, 379; Phys. Rev. Lett., 1996, 76, 928. [58] Luzar, A. J. Chem. Phys. 2000, 113, 10663. [59] Kittel, C. Introduction to Solid State Physics, Wiley & Sons, Inc., New York, 2nd edition, 1953. [60] Kroes, G. J. Surf. Sci. 1992, 275, 365. 56  Chapter 3. Bibliography [61] Turnbull, D.; Vonnegut, B. Ind. Eng. Chem. 1952, 44, 1292.  57  Chapter 4 Water Adsorption on Kaolinite Surfaces Containing Trenches∗ 4.1  Introduction  Every year between 800 − 2000 Tg of mineral dust particles are emitted into the atmosphere, which is a major fraction of the total global aerosol emission [1–6]. The origin of most aerosolized mineral dust is attributed to large dust storms occurring over the desert and arid regions of the west coast of North Africa, the Middle East, Central and South Asia, and China [7]. Aerosolized mineral dust is mainly composed of illite, chlorite, montmorillonite, kaolinite, quartz and calcite, with kaolinite comprising approximately 5−10% [8, 9]. Monitoring of dust events has shown that these particles can be transported over large distances with atmospheric lifetimes on the order of days [1, 10]. In the atmosphere mineral dust particles can take up water and act as cloud condensation nuclei or ice nuclei [11–23]. Through these mechanisms mineral dust particles can influence climate by changing the properties of clouds [24]. A better understanding of the interactions of water with mineral dust particles should lead to a better understanding of the cloud condensation and ice nucleation properties of these particles. The chemistry of the atmosphere can also be influenced by mineral dust particles as they provide a medium on which reactions can occur. These heterogeneous reactions can modify the chemical balance of gas-phase atmospheric species, such as N2 O5 , nitric acid and ozone, and the chemical composition of individual particles. Recently, it was shown that relative humidity (RH) (defined as P/Po ×100, where P is the vapor pressure of water A version of this chapter has been accepted for publication. T. Croteau, A. K. Bertram, G. N. Patey, Water Adsorption on Kaolinite Surfaces Containing Trenches. ∗  58  4.1. Introduction and Po is the saturation vapor pressure of liquid water) influences reactions involving these trace atmospheric gases. For example, Mogili et al. [25] found a significant increase in the uptake and rate of heterogeneous hydrolysis of N2 O5 on quartz as the relative humidity increased. Similarly, enhanced uptake of nitric acid with increasing relative humidity was reported on Na-montmorillonite [26], oxide [27], carbonate [28–31], and dust samples [32– 34]. To the contrary, studying the effect of adsorbed water on the reaction of ozone with alumina and hematite surfaces, Mogili et al. [25] found a strong decrease in ozone uptake as the relative humidity was increased. A better understanding and quantification of the interactions of water with mineral surfaces should also lead to a better understanding of these different effects. Molecular simulations are well suited for studying water adsorption on mineral surfaces. This technique was used recently in our investigation of water adsorption on atomistically smooth kaolinite particles [35, 36]. In agreement with ab initio studies, we found that the Si- and Al- surfaces are hydrophobic and hydrophilic, respectively [35–41]. In addition, we found that some edges (lateral faces) are very hydrophilic, taking up water at very low relative humidity. However, the maximum coverage observed on any surface studied on kaolinite did not reach more than slightly overgrown monolayers below saturation (i.e., RH < 100%). This is in clear contrast with the recent study of Schuttlefield et al. [42], who reported, using ATR-FTIR spectroscopy coupled with quartz crystal microbalance (QCM) measurements, that large amounts of water (multilayers) are adsorbed on kaolinite and other clay minerals at RH values below 100%. The above experimental and simulation results support the idea that surface defects or irregularities such as steps, cracks, trenches, etc., must play a major role in the adsorption of water occurring on clay minerals. Related to the above discussion, ice freezing experiments indicate that ice nucleation on mineral dust surfaces also occurs at specific types of surface defects, so called active ˚2 in sites [43–47]. For example, Marcolli et al. [43] showed that active sites (about 10A  59  4.2. Model and Method size) were needed to describe freezing experiments on mineral dust. Using grand canonical Monte Carlo (GCMC) simulations, Cailliez et al. [48] investigated the effect of surface defects on the hydration thermodynamic behavior of silicalite-1 zeolite. They found a change in the hydration behavior from “hydrophobic” to “hydrophilic” as the number of attractive surface defects was increased. Shevkunov [49] also studied the effect of surface defects on the adsorption of water on a AgI crystal using a Monte Carlo method. He observed that crystal defects of size 15 − 20˚ A most effectively stimulate the formation of a macroscopic condensed phase. Following our previous study on atomistically smooth kaolinite surfaces, we now investigate water adsorption on kaolinite surfaces containing defects using GCMC calculations. We examine water adsorption and structure at 298 K. We discuss binding energies, and compare our simulations with recent laboratory studies of water uptake on kaolinite surfaces. One of the goals of this study is to determine if defects can lead to larger than monolayer coverages, and potentially explain the discrepancy between laboratory measurements and simulations on flat surfaces. The particular defects that we consider are nanoscale trenches with widths ranging from 14.78 ˚ A to 73.91 ˚ A and a depth of 23.19 ˚ A. The motivation for this geometry and these length scales is discussed below. The remainder of the paper is organized in three parts. The models and simulation issues are considered in Section II, results are presented and discussed in Section III, and our conclusions are summarized in Section IV.  4.2  Model and Method  The models employed to represent the clay lattice (CLAYFF) and water molecules (SPC/E) were used in earlier studies and have been described elsewhere [35, 50–52]. The total configurational energy of our system, U, is given by the sum of Coulombic and Lennard-Jones  60  4.2. Model and Method (LJ) interactions e2 U= 4πǫo  i=j  qi qj + rij  Do,ij i=j  Ro,ij rij  12  Ro,ij −2 rij  6  ,  (4.1)  where the sums are over all lattice and water interaction sites, and rij is the distance between sites i and j. The first term represents the Coulombic interactions, where qi is the point charge on site i, e is an elementary charge, and ǫo is the dielectric permittivity of free space. The second term represents the LJ contribution and the energy and distance parameters for the different site-site interactions, Do,ij and Ro,ij , respectively, are given in Ref. [50]. For the cross interactions, the Lorentz-Berthelot rules, Do,12 =  Do,11 Do,22  and Ro,12 = (Ro,11 + Ro,22 )/2, are used. As mentioned above, the defects we consider are so-called trenches, as shown in Figs. 4.1 and 4.2. The motivation for this particular defect type comes from scanning electron microscope (SEM) images of kaolinite particles [53–59]. Such images show that the basal planes of kaolinite are relatively smooth (at least on the length-scale of the resolution of an SEM) but the surfaces perpendicular to the basal planes are rough and irregular, with features that often look like trenches or cracks that should have sides (i.e., walls) that correspond to the Al- and Si-surfaces, similar to our design (Figs. 4.1 and 4.2). The trench bottoms and the surface areas between the trenches consist of unprotonated edges (100 plane). In the SEM images [53–59] the dimensions of the trenches and cracks are larger than the dimensions used in our simulations, but it seems reasonable to assume that similar structures are present at smaller length-scales below the resolution of SEM. Due to computational limitations, trenches of the dimensions observed by SEM cannot be simulated. However, trenches on the size of several nanometres are within the limits of current computational power. A rigid unit cell corresponding to the most stable kaolinite [Al2 Si2 O5 (OH)4 ] configuration, with parameters a = 5.1535 ˚ A, b = 8.9419 ˚ A, and c = 7.3906 ˚ A, was used to construct the trenches [60]. As seen in Fig. 4.1, the trenches have finite width (x-axis) 61  4.2. Model and Method  EMPTY SPACE  1111111111111111 0000000000000000 0000000000000000 1111111111111111 0000000000000000 1111111111111111 00000000000000000 11111111111111111 0000000000000000 1111111111111111 00000000000000000 11111111111111111 00000000000000000 11111111111111111 0000000000000000 1111111111111111 00000000000000000 11111111111111111 00000000000000000 11111111111111111 0000000000000000 1111111111111111 (100) 00000000000000000 11111111111111111 00000000000000000 11111111111111111 0000000000000000 1111111111111111 00000000000000000 11111111111111111 TRENCH 11111111111111111 WATER 00000000000000000 0000000000000000 1111111111111111 plane 00000000000000000 11111111111111111 00000000000000000 11111111111111111 0000000000000000 1111111111111111 00000000000000000 11111111111111111 00000000000000000 11111111111111111 0000000000000000 1111111111111111 00000000000000000 11111111111111111 000000000000000001111111111111111 11111111111111111 0000000000000000 00000000000000000 11111111111111111  z  y x  EMPTY SPACE Figure 4.1: The simulation cell.  62  4.2. Model and Method  Figure 4.2: The trenches considered at 298 K. Note that the systems are replicated along the x axis to match the dimension of the largest case. From top to bottom, the trenches are 1, 2, 3 and 4. Simulations were performed for trench 1 without (1A) and with (1B) the inter-trench spacing filled (clay atoms inside the black rectangle). The black arrows indicate the trenches.  63  4.2. Model and Method and depth (z-axis) but are made infinite in the y dimension by periodically replicating the system. The trench walls are composed of Al- and Si- surfaces and the horizontal surfaces (xy plane) correspond to unprotonated edge surfaces [35, 36]. These surfaces are constructed by replicating the unit cell as many times as needed to attain the desired dimensions along a Miller index that will produce an Al-surface, a Si-surface, or an unprotonated edge. All trenches we consider are based on Fig. 4.1, and only vary in the x dimension of the unprotonated edge segments (bottom of the trench and top surface). We constructed four main trenches labeled 1, 2, 3, and 4, all having different dimensions based on the above unit cell parameters (see Fig. 4.2 and Table 4.1). The first step in building a trench is to construct an Al-surface with x and y dimensions of 30.921 ˚ A and 35.7676 ˚ A, respectively. The surface is then rotated 90o counterclockwise to bring the basal surfaces along the z axis (the Si-surface forms the inner trench wall). Then, an unprotonated edge surface with the same y dimension as the Al-surface and varying x dimension, depending on the trench (see Fig. 4.2), is placed to the right of the Si-surface so that it forms the floor of the trench. Then, another Al-surface is placed at the end of the unprotonated edge surface to form the second trench wall. Finally, immediately to the right of the top of the second trench wall is placed another unprotonated edge surface of x dimension of 14.78 ˚ A (trenches 1, 3, and 4) or 29.56 ˚ A (trench 2). A complete list of all trench parameters is given in Table 4.1. There are two types of trenches with parameters corresponding to those of trench 1. One which has the inter-trench spacing empty (1A) and the other (1B) filled with clay atoms (see black outline in Fig. 4.2). The dimensions of trench 1 differ from trench 2 only in the length of the top unprotonated edge surface. Trenches 1 and 2, 3, and 4, differ in width. The trenches all have a depth of 23.19 ˚ A. The double slab systems (Fig. 4.1) contain 4896 (trench 1A, inter-trench spacing empty), 7344 (trench 1B, inter-trench spacing filled), 5712 (trench 2), 6528 (trench 3), and 8160 (trench 4) clay atoms. A minimum of ∼ 1300 water molecules and up to ∼ 7000 were found in our systems for the 64  4.2. Model and Method range of chemical potentials studied and different trench types. As discussed in Ref. [36], two opposing clay slabs were used in order to cancel the longrange field due to arbitrary truncation in the z dimension (Fig. 4.1). Unless otherwise specified, 50 ˚ A separates the closest two unprotonated edge surfaces. Fig. 4.3 shows the energy, averaged over orientations of a single water molecule at 298 K, as a function of distance from the bottoms of trenches 1A and 1B. These average energies were collected over 1 × 105 Monte Carlo rotations for each position. The water-lattice interaction energy remains very attractive as long as the water molecule is located inside the trench. When the water molecule reaches the top of the trench (∼ 30 ˚ A) the interaction energy increases rapidly. In both cases, the energy approaches zero about 10-15 ˚ A above the top of the trench. Note that as the water molecule moves away from the trench walls (green curve in Fig. 4.3) the interaction energy decreases. Thus, as the trench becomes wider the interaction energy near the center of the trench decreases. Quantitative differences are observed in the strength of the interaction for trenches 1A and 1B. Trench 1B leads to more attractive interactions inside and near the top of the trench. The effect that this has on the adsorption of water molecules is discussed below. The configurational energies were obtained using the usual 3D Ewald method [61] leaving an empty gap of 107 ˚ A in the z dimension to avoid any unphysical interaction with the other periodic images [62, 63]. The trenches were repeated periodically in the x and y dimensions. A total of 6858 wave vectors were employed in the Ewald sums ˚−1 , and a reciprocal space cut-off of 2.0 A ˚−1 . together with the parameter α = 0.175 A The LJ and real space interactions were cut off at 18.0 ˚ A. The GCMC method was employed, using the procedure previously described [35, 36]. Insertions, deletions and displacements (each move corresponding to one step) of water molecules were attempted with equal probability at fixed chemical potential. Very long simulations were necessary to attain equilibrium (on the order of 109 or more steps) due to the large number of water molecules present. The number of clay atoms also significantly  65  4.3. Results and Discussion Table 4.1: Dimensions of the simulation cells in terms of the lattice parameters, a = 5.1535 ˚ A, b = 8.9419 ˚ A, and c = 7.3906 ˚ A. The y dimension and trench depth are 4b ˚ ˚ (35.77 A) and 4.5a (23.19 A), respectively, for all systems. Trench trench width x dimension (˚ A) (˚ A) 1A, 1B 2c (14.78) 6c (44.34) 2 2c (14.78) 8c (59.12) 3 6c (44.34) 10c (73.91) 4 10c (73.91) 14c (103.47) affected the computing time required (on the order of several weeks or even a few months). Based on our previous studies [35, 36], it was judged sufficient (for qualitative purposes) to have a single layer of clay atoms around the trench structure, thus minimizing the computing cost. Nevertheless, a comparison is presented for trench 1 with and without the inter-trench spacing filled (see Fig. 4.2 and discussion below). Following equilibration, averages were calculated including the energy and water content. Adsorption isotherms at 298 K were obtained for all trenches.  4.3 4.3.1  Results and Discussion Water Adsorption and Structure at 298 K  Adsorption isotherms obtained at 298 K are presented in Fig. 4.4 as functions of relative humidity. As the vapor pressures involved in the simulations are low (below 0.1 atm), a satisfactory estimate of the pressure is obtained using the ideal gas relationship, P = eβµtr /βΛ3 , where Λ is the thermal de Broglie wavelength and µtr is the translational contribution to the chemical potential. P/Po is set to one at the bulk condensation of SPC/E water at 298 K (µ = −58.5 kJ/mol) [64]. In Fig. 4.4, the y axis is plotted as the number of monolayers in order to facilitate direct comparison with experiment [42]. The number of water molecules in a monolayer is defined as the total surface area (including horizontal surfaces and walls of the trenches) divided by 9˚ A2 , the area occupied  66  4.3. Results and Discussion  0  E (kJ/mol)  0 -20  -20  -40  -40  -60  -60  -80  -80  -100  -100  -120  -120  -140  -140  -160  -160  -180  -180  -200  -200  -220  0  20  40  60  80  100  20  40  60  80  100  -220  z (Å) Figure 4.3: The interaction energy averaged over the orientations of a single water molecule as a function of distance from the bottom of trench 1A (left) and 1B (right). The red and blue curves correspond to a water molecule moving up along both trench walls, while the green curve corresponds to a water molecule moving up the middle of the trench.  67  4.3. Results and Discussion by a water molecule. Therefore, the estimated number of monolayers is given by 9˚ A2 · Navg /Atot , where Navg and Atot are the average number of water molecules and total surface area, respectively. Note that using this definition a monolayer has 1.1 × 1015 water molecules/cm2 . From our previous results, the maximum coverage obtained on atomistically smooth kaolinite surfaces at 298 K corresponds to slightly overgrown monolayers [36]. This is in clear contrast with the current results which suggest multilayer coverages in all cases (Fig. 4.4), if one were to assume that the water molecules are uniformly distributed. As seen in Figs. 4.5 and 4.6, trenches which are 14.78 ˚ A wide and 23.19 ˚ A deep (trenches 1A, 1B and 2) are always full of water even at very low vapor pressures (0.0003% RH). Due to this filling, the apparent coverage at 1% and higher is greater than 2 monolayers. This clearly shows that trenches can lead to enhanced water uptake, and apparent multilayer adsorption. Effect of Slab Separation Tests were made to check the effect of slab separation in our two-slab geometry (Fig. 4.1) by doubling the separation distance from 50 ˚ A to 100 ˚ A for all trenches considered at selected relative humidities. The results obtained in the relevant RH regime show only minor qualitative differences. For example, the number of water molecules adsorbed for trench 1B at RH values ranging from 24% to 82% is almost the same at both slab separations, as can be seen from the results (purple stars) included in Fig. 4.4. Similar minor changes were observed for the other trenches in this RH range. The effect of slab separation will be a little more important near saturation, affecting the pressure at which “bulk” condensation (slit filling) occurs. As we are currently not focusing on the exact location of the filling transition, but rather on the adsorption behavior over a wide range of relative humidities below saturation, the results obtained with a slab separation of 50 ˚ A are satisfactory.  68  4.3. Results and Discussion  6  Monolayers  5  4  3  2  1  0  0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1  P/Po Figure 4.4: Adsorption isotherms obtained at 298 K for trenches 1A (red triangles up), 1B (orange stars), 2 (black circles), 3 (green squares), and 4 (maroon plus signs) with a slit separation of 50 ˚ A, as functions of RH. Results for an atomistically smooth unprotonated edge surface [36] (blue triangles down) are included for comparison. The results plotted as purple stars were obtained for trench 1B at a slit separation of 100 ˚ A.  69  4.3. Results and Discussion Effect of Clay Atoms in the Inter-trench Space Comparing the results for trenches 1A and 1B, we see that filling the inter-trench spacing with clay atoms (trench 1B) leads to somewhat enhanced water uptake, by a factor of ∼ 1.2 at all RH values (see Fig. 4.4). Thus, at the same relative humidities more water molecules are adsorbed for trench 1B. At 0.0003% RH, trench 1A has a dry upper unprotonated edge surface, whereas it is partially covered with water molecules for trench 1B (Figs. 4.5a and b). Both cases have small mounds at the step edges. At higher relative humidities, greater mounds are present at the step edges which form a thick inter-connecting coverage at 82% RH for trench 1B (Fig. 4.5b). Even though larger coverages are observed when the inter-trench spacing is filled, the general adsorption behaviour is qualitatively similar for both models. The single clay layer models are more feasible computationally, and give a reasonable qualitative description of the adsorption. Quantitatively, one can take the view that these models provide lower limits on the amount of water adsorbed. Effect of Inter-trench Separation Length Doubling the length of the top unprotonated edge separating the periodic trench images in the x dimension from 14.78 ˚ A (trench 1A) to 29.56 ˚ A (trench 2) leads to a small decrease in coverage by a factor of ∼ 1.1. Although trench 2 takes up more water molecules on average, the apparent monolayer coverage is slightly smaller. This is mainly due to the extra contribution to the total surface area coming from the longer top unprotonated edge of trench 2, which is covered only by a single monolayer. As noted above, both trenches 1A and 2 are completely filled in the very low RH regime (0.0003% RH), with small mounds at the step edges and dry upper unprotonated edge surface (Figs. 4.5a and 4.6). As the relative humidity is increased to 1.0%, the step-edge mounds increase in size and water molecules begin to partially cover the top unprotonated edge (middle images in Figs. 4.5a and 4.6). Finally, at 82% RH, large mounds have grown at the step edges, and a full monolayer covers the top unprotonated edge surface, consistent with  70  4.3. Results and Discussion  0.0003%  0.0003%  1.0%  1.0%  82%  82%  (a)  (b)  Figure 4.5: Snapshots of trench 1A (left) and 1B (right) as functions of RH at 298 K. The corresponding chemical potentials are µ = −90 kJ/mol (top), −70 kJ/mol (middle), and −59 kJ/mol (bottom). The water oxygen atoms are blue, lattice oxygen atoms are red, hydrogen atoms are white, aluminum atoms are grey, and silicon atoms are brown.  71  4.3. Results and Discussion our previous study [36]. Overall, the effect of having a smaller top unprotonated edge separation between periodic trench images is minimal, only affecting adsorption very near bulk condensation. Effect of Trench Width Increasing the trench width by a factor of 3 from 14.78 ˚ A (trench 1A, 1B and 2) to 44.34 ˚ A (trench 3) significantly influences the water coverage. The wider trench 3 is not full in the low relative pressure regime (below 11% RH). For example, at ∼ 0.0003% RH the trench is only partially filled with thick water aggregates located at both step edges and has a dry unprotonated edge at the bottom (see Fig. 4.7). At higher RH, more water molecules are taken up in the trench leading to larger coverages. At 11% RH, which is in the range of atmospheric relevance, trench 3 fills completely and large mounds develop at the step edges (Fig. 4.7, middle image). At 82% RH, the large mounds (∼ 16 ˚ A thick) formed at the step edges grow to connect the periodic trenches. Trench 3 leads to larger apparent coverages than the narrower trenches 1 and 2, ranging from ∼ 1.2 (0.0003% RH) to ∼ 5.6 monolayers (82% RH). The widest trench we consider (trench 4) is 73.93 ˚ A in width. Snapshots for this trench are shown in Fig. 4.8, and we note that in this case the trench is completely filled at 82% RH. However, at 11%, and particularly at 55% RH, large mounds have grown at the step edges. At 82% RH, the apparent coverage is ∼ 6.3 monolayers.  4.3.2  Binding Energies  As is evident from the discussion above, the trenches we consider are strongly hydrophilic, and can take up large quantities of water. This is due to the fact that inside a trench, as shown in Fig. 4.3 for a single water molecule, the water-lattice interactions are very strong. As an example, Table 4.2 gives the average water-water and water-lattice energies for trench 2 at various RH. Standard deviations based on ten “independent” estimates  72  4.3. Results and Discussion  0.0003%  1.0%  82%  Figure 4.6: Snapshots of trench 2 as functions of RH at 298 K. The corresponding chemical potentials are µ = −90 kJ/mol (top), −70 kJ/mol (middle), and −59 kJ/mol (bottom). The atoms are denoted as in Fig. 4.5.  73  4.3. Results and Discussion  0.0003%  11%  82%  Figure 4.7: Snapshots of trench 3 as functions of RH at 298 K. The corresponding chemical potentials are µ = −90 kJ/mol (top), −64 kJ/mol (middle), and −59 kJ/mol (bottom). The atoms are denoted as in Fig. 4.5.  74  4.3. Results and Discussion  11%  55%  82%  Figure 4.8: Snapshots of trench 4 as functions of RH at 298 K. The corresponding chemical potentials are µ = −64 kJ/mol (top), −62 kJ/mol (middle), and −59 kJ/mol (bottom). The atoms are denoted as in Fig. 4.5.  75  4.3. Results and Discussion Table 4.2: Average configurational energies at 298 K for trench 2. The numbers in brackets are estimated standard deviations. Relative µ Uwater-water Uwater-surface Utotal Humidity (kJ/mol) (kJ/mol) (kJ/mol) (kJ/mol) 82% -59 -26.68(0.27) -65.36(0.38) -92.04(0.52) 55% -60 -25.64(0.19) -71.28(0.75) -96.92(0.57) 24% -62 -26.18(0.10) -75.29(0.10) -101.47(0.11) 11% -64 -26.05(0.10) -76.23(0.12) -102.29(0.09) 1% -70 -23.69(0.22) -80.92(0.78) -104.61(0.67) 0.0003% -90 -21.34(0.13) -96.03(0.13) -117.37(0.14) are included. Note that all water molecules found in these systems, sitting on the top unprotonated edge surface or inside the trench, were included in the average. Below saturation, the water-lattice energy of interaction is very attractive ranging from ∼ −65 kJ/mol (P/Po ≈ 0.82) to ∼ −95 kJ/mol (P/Po ≈ 3 × 10−6). On the other hand, the water-water interactions make a smaller contribution to the total energy and are of the order of -20 to -30 kJ/mol. As a reference, bulk water using SPC/E with 3D periodicity has an average internal energy of ∼ −47 kJ/mol. The relatively small water-water energy for our structured surfaces is largely due to the finite size of these systems. The strong water-lattice interaction energies in the trenches really emphasizes the importance of the granularity of the clay particles for apparent multilayer adsorption to occur. They also suggest that deeper trenches of comparable widths would also fill completely at the RH values considered here.  4.3.3  Comparison of Calculated Adsorption Isotherms with Experiment  In their experimental study using ATR-FTIR spectroscopy coupled with quartz crystal microbalance (QCM) measurements, Schuttlefield et al. [42] reported that at 296 K multiple monolayers of water are adsorbed. In fact, based on their results at 50% and 80% RH, approximately 6-25 and 10-40 monolayers of water, respectively, are uptaken depending on the source of the kaolinite sample. They expressed the adsorbed mass 76  4.3. Results and Discussion of water as multilayer coverage based on the estimated total surface area of their clay samples. This was obtained by multiplying the mass of the clay, as obtained with the QCM, by the corresponding specific BET surface area. Finally, the mass of water was converted into monolayers by approximating a monolayer as 1×1015 water molecules/cm2 . (Note that we assume 1.1 × 1015 water molecules/cm2 which gives slightly lower coverage estimates.) It is clear the single monolayer coverage previously obtained for atomistically smooth kaolinite surfaces [35, 36] is inconsistent with these experiments. Our current simulations (Fig. 4.4) for the largest trench considered gives an apparent coverage of up to ∼ 6.3 monolayers, which is closer to experiment than the single monolayers obtained on atomistically smooth surfaces. Due to computational limitations, we have considered only trenches that are 23.2 ˚ A deep. If we assume that deeper trenches of the same width would show similar behavior (which is not unlikely given the strength of the water-lattice interactions in a trench), we can extrapolate to obtain apparent coverages for deeper trenches. For rectangular trenches as in Fig. 4.1, the apparent monolayer coverage (number of monolayers) is given by the simple expression, A2 /(2 · zt + xt + xs ), where ρt and ρs are the volume and surface (ρt · xt · zt + ρs · xs )·9˚ water densities inside the trench and on the top unprotonated edge, respectively, xt and zt are the trench width and depth, and xs is the length of the top unprotonated edge in the x dimension. This expression has two obvious limits: if the trench depth zt is much A2 /2; if the trench width xt larger than xt and xs , the expression simplifies to ρt · xt · 9˚ A2 . Using the first simplification is much larger than zt and xs , it reduces to ρt · zt · 9˚ (much deeper trench), and estimating ρt from our simulation results, we can obtain upper limits for the different trench widths considered here. This is not inconsistent with SEM imaging [53–59] as the trench widths and inter-trench spacings are often much smaller than the trench depths. The trench density ρt is obtained by dividing the average number of water molecules found below the top of a trench by the corresponding trench volume. For trenches 1 and 2, 3, and 4 the estimated upper limits are approximately 3.6, 10.4, and  77  4.4. Summary and Conclusions 14.8 monolayers, respectively, which are closer again to the reported experimental results [42]. These estimates based on our simulations serve to illustrate the substantial influence that surface defects can have on water uptake on clay particles. Clearly, defects must be considered in attempting to understand their water uptake, and the associated atmospheric implications. It is important to note that we have studied only one specific type of defect over a limited range of dimensions. More simulations on other possible types of defects are needed, as well as experiments to better characterize the types of defects that occur on real mineral dust particles.  4.4  Summary and Conclusions  Water adsorption simulations were performed at 298 K on four distinct kaolinite trenches using the GCMC method. The water molecules were represented by the SPC/E model. The interaction of the water molecules with the kaolinite lattice was found to be very attractive inside the trenches. This leads to apparent multilayer coverages in all our systems for the whole range of atmospherically relevant relative humidities. This is in clear contrast with our previous results on atomistically smooth kaolinite surfaces, where the maximum coverage observed at 298 K was a slightly overgrown monolayer. For trenches 1 and 2 (∼ 14 ˚ A wide) the interior of the trench remains filled with water even at very low relative humidities. In agreement with Ref. [36], the top unprotonated edge surface picks up water gradually as the relative humidity is increased finally forming a monolayer before saturation. Mounds of water molecules are located at the step edges, and these grow larger as the relative humidity is increased. From a modeling perspective, truncation of the kaolinite lattice mainly affects the strength of the water-lattice interaction inside a trench. If more lattice atoms are included, the water-lattice interaction tends to be stronger, which in turn tends to increase the coverage at a given RH value. Nevertheless, qualitatively, the general adsorption trends do not show a strong dependence on 78  4.4. Summary and Conclusions lattice truncation. This is mainly because inside a trench the water-lattice interactions are very strong even for the most severe (also the most economical) lattice truncation, where the trenches are surrounded by a single kaolinite layer. The larger trenches considered (trenches 3 and 4) either partially or completely fill in the atmospherically relevant RH range. In these wider trenches the water-lattice interaction decreases slowly away from the trench wall, and this reduces the adsorption at lower RH values. Nonetheless, these systems also develop thick mounds of water at the trench edges and fill completely below saturation. This yields apparent multilayer coverages, up to 6.3 for trench 4. Extrapolation of our results to deeper trenches gives coverage estimates that lie in, or are close to, the ranges reported experimentally [42]. If our simulations are applicable to real mineral dust particles, they show that large amounts of water can accumulate in trenches, and likely in other similar defects, under atmospherically relevant RH values. These regions of high water content would seem to be ideal locations for heterogeneous chemistry such as N2 O5 hydrolysis to occur. We would expect the environment in these “pools” of water to be significantly different from flat surfaces, where only monolayers form. From the perspective of classical nucleation theory, these water pools also suggest a favorable location for ice nucleation. Ice nucleation could occur deep within a trench, where the ice-water(liquid) interface would possibly pose less of a barrier to nucleation, than the ice-air interface that is expected to contribute to the barrier on flat surfaces. This interesting possibility will be the focus of a future study (see Chapter 5).  79  Bibliography [1] Bauer, S. E.; Balkanski, Y.; Schulz, M.; Hauglustaine, D. A.; Dentener, F. J. Geophys. Res.-Atm. 2004, 109 (D02304). [2] Lunt, D. J.; Valdes, P. J. J. Geophys. Res.-Atm. 2002, 107 (D23). [3] Tegen, I.; Harrison, S. P.; Kohfeld, K. E.; Prentice, I. C.; Coe, M.; Heimann, M. J. Geophys. Res.-Atm. 2002, 107 (D21), 4576. [4] Tegen, I.; Miller, R. J. Geophys. Res.-Atm. 1998, 103 (D20), 25975. [5] Woodward, S. J. Geophys. Res.-Atm. 2001, 106 (D16), 18155. [6] Zender, C. S.; Bian; H.; Newman, D. J. Geophys. Res. 2003, 108 (D14), 4416. [7] Prospero, J. M.; Ginoux, P.; Torres, O. ; Nicholson, S. E.; Gill, T. E. Rev. Geophys. 2002, 40 (1), 1002. [8] Usher, C. R.; Michel, A. E.; Grassian, V. H. Chem. Rev. 2003, 103 (12), 4883. [9] Glaccum, R. A.; Prospero, J. M. Marine Geology 1980, 37 (3-4), 295. [10] Husar, R. B.; et al. J. Geophys. Res. 2001, 106, 18317. [11] Cziczo, D. J.; Murphy, D. M.; Hudson, P. K.; Thomson, D. S. J. Geophys. Res.-Atm. 2004 109 (D4). [12] DeMott, P. J.; Cziczo, D. J.; Prenni, A. J.; Murphy, D. M.; Kreidenweis, S. M.; Thomson, D. S.; Borys, R.; Rogers, D. C. Proc. Nat. Acad. Sci. 2003, 100 (25), 14655. [13] DeMott, P. J.; Sassen, K.; Poellot, M. R.; Baumgardner, D.; Rogers, D. C.; Brooks, S. D.; Prenni, A. J.; Kreidenweis, S. M. Geophys. Res. Let. 2003, 30, 1732. [14] Eastwood, M. L.; Cremel, S.; Gehrke, C.; Girard, E.; Bertram, A. K. J. Geophys. Res.-Atm. 2008, 113, D22203. [15] Eastwood, M. L.; Cremel, S.; Wheeler, M.; Murray, B. J.; Girard, E.; Bertram, A. K. Geophys. Res. Let. 2009, 36, L02811. [16] Dymarska, M.; Murray, B. J.; Sun, L.; Eastwood, M. L.; Knopf, D. A.; Bertram, A. K. J. Geophys. Res.-Atm. 2006, 111 D4.  80  Chapter 4. Bibliography [17] Twohy, C. H.; Kreidenweis, S. M.; Eidhammer, T.; Browell, E. V.; Heymsfield, A. J.; Bansemer, A. R.; Anderson, B. E.; Chen, G.; Ismail, S.; DeMott, P. J.; Van den Heever, S. C. Geophys. Res. Let. 2009, 36. [18] Twohy, C. H.; Poellot, M. R. Atm. Chem. Phys. 2005, 5, 2289. [19] Herich, H.; Tritscher, T.; Wiacek, A.; Gysel, M.; Weingartner, E.; Lohmann, U.; Baltensperger, U.; Cziczo, D. J. Phys. Chem. Chem. Phys. 2009, 11, 7804. [20] Welti, A.; Lnd, F.; Stetzer, O.; Lohmann, U. Atm. Chem. Phys. 2009, 9, 6705. [21] Zimmermann, F.; Weinbruch, S.; Schtz, L.; Hofmann, H.; Ebert, M.; Kandler, K.; Worringen, A. J. Geophys. Res.-Atm. 2008, 113, D23204. [22] Connolly, P. J.; Mhler, O.; Field, P. R.; Saathoff, H.; Burgess, R.; Choularton, T.; Gallagher, M. Atm. Chem. Phys. 2009, 9, 2805. [23] Ansmann, A.; Tesche, M.; Althausen, D.; Mller, D.; Seifert, P.; Freudenthaler, V.; Heese, B.; Wiegner, M.; Pisani, G.; Knippertz, P.; Dubovik, O. J. Geophys. Res.Atm. 2008, 113, D04210. [24] Forster, P.; Ramaswamy, V.; Artaxo, P.; Berntsen, T.; Betts, R.; Fahey, D. W.; Haywood, J.; Lean, J.; Lowe, D. C.; Myhre, G.; Nganga, J.; Prinn, R.; Raga, G.; Schulz, M.; Dorland, R. V. Changes in Atmospheric Constituents and in Radiative Forcing., in Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change 2007, edited by S. Solomon, D. Qin, M. Manning, Z. Chen, M. Marquis, K.B. Averyt, M. Tignor, and H.L. Miller, pp. 129-234, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA. [25] Mogili, P. K.; Kleiber, P. D.; Young, M. A.; Grassian, V. H. J. Phys. Chem. A 2006, 110, 13799. [26] Mashburn, C. D.; Frinak, E. K.; Tolbert, M. A. J. Geophys. Res.-Atm. 2006, 111 (D15). [27] Goodman, A. L.; Bernard, E. T.; Grassian, V. H. J. Phys. Chem. A 2001, 105, 6443. [28] Al-Hosney, H. A.; Grassian, V. H. Phys. Chem. Chem. Phys. 2005, 7, 1266. [29] Goodman, A. L.; Underwood, G. M.; Grassian, V. H. J. Geophys. Res. 2000, 105, 29053. [30] Krueger, B. J.; Ross, J. L.; Grassian, V. H. Langmuir 2005, 21, 8793. [31] Santschi, C.; Rossi, M. J. J. Phys. Chem. A 2006, 110 (21), 6789. [32] Hanisch, F.; Crowley, J. N. Phys. Chem. Chem. Phys. 2001, 3 (12), 2474. [33] Laskin, A.; Wietsma, T. W.; Krueger, B. J.; Grassian, V. H. J. Geophys. Res. 2005, 110, D10208. 81  Chapter 4. Bibliography [34] Underwood, G. M.; Li, P.; Al-Abadleh, H.; Grassian, V. H. J. Phys. Chem. A 2001, 105, 6609. [35] Croteau, T.; Bertram, A. K.; Patey, G. N. J. Phys. Chem. A 2008, 112(43), 10708. [36] Croteau, T.; Bertram, A. K.; Patey, G. N. J. Phys. Chem. A 2009, 113(27), 7826. [37] Tunega, D.; Gerzabek, M. H.; Lischka, H. J. Phys. Chem. B 2004, 108, 5930. [38] Tunega, D.; Benco, L.; Haberhauer, G.; Gerzabek, M. H.; Lischka, H. J. Phys. Chem. B 2002, 106, 11515. [39] Tunega, D.; Haberhauer, G.; Gerzabek, M. H.; Lischka, H. Langmuir 2002, 18, 139. [40] Hu, X. L.; Michaelides, A. Surf. Sci. 2007, 601, 5378. [41] Hu, X. L.; Michaelides, A. Surf. Sci. 2008, 602, 960. [42] Schuttlefield, J. D.; Cox, D.; Grassian, V. H. J. Geophys. Res. 2007, 112 (D21303). [43] Marcolli, C.; Gedamke, S.; Peter, T.; Zobrist B. Atmos. Chem. Phys. 2007, 7, 5081. [44] Archuleta, C. M.; DeMott, P. J.; Kreidenweis, S. M. Atmos. Chem. Phys. 2005, 5, 2617. [45] Hung, H.-M.; Malinowski, A.; Martin S. T. J. Phys. Chem. A 2003, 107, 1296. [46] Welti, A.; Luond, F.; Stetzer, O; Lohmann, U. Atmos. Chem. Phys. Discuss. 2009, 9, 6929. [47] Pruppacher, H. R.; Klett, J. D. Microphysics of Clouds and Precipitation, 2nd rev. ed.; Kluwer Academic: Dordrecht, The Netherlands, 1997; pp 329-331. [48] Cailliez, F; Stirnemann, G.; Boutin, A.; Demachy, I.; Fuchs, A. H. J. Phys. Chem. C 2008, 112, 10435. [49] Shevkunov, S. V. Colloidal Journal 2006, 68, 370. [50] Cygan, R. T.; Liang, J.-J.; Kalinichev, A. G. J. Phys. Chem. B 2004, 108, 1255. [51] Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981; pp 331. [52] Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. [53] Papoulis, D.; Tsolis-Katagas, P. Clay Min. 2008, 43, 631. [54] Valaskova, M.; Rieder, M.; Matejka, V.; Capkova, P.; Sliva, A. Appl. Clay. Sci. 2007, 35, 108. [55] Aparicio, P.; Perez-Bernal, J. L.; Galan, E.; Bello, M. A. Clay Min. 2004, 39, 75. 82  Chapter 4. Bibliography [56] Ekosse, G. Appl. Clay Sci. 2000, 16, 301. [57] Keller, W. D. Clays Clay Miner. 1976, 24, 107. [58] Keller, W. D. Clays Clay Miner. 1978, 26, 1. [59] Bauluz, B.; Mayayo, M. J.; Yuste, A.; Gonzalez Lopez, J. M. Clay Min. 2008, 43, 459. [60] Bish, D. L. Clays Clay Miner. 1993, 41, 783. [61] Allen, M. P.; Tildesley, D. J. Computer simulation of liquids, (Oxford, 1987). [62] Shelley, J. C.; Patey, G. N. Mol. Phys. 1996, 88, 385. [63] Yeh, I.-C.; Berkowitz, M. L. J. Chem. Phys. 1999, 111, 3155. [64] Liu, J. C.; Monson, P. A.; Langmuir 2005, 21, 10219.  83  Chapter 5 Structural Analysis of Water Adsorbed on Kaolinite Particles Containing Trenches at 235 K and Possible Connections with Ice Nucleation in the Atmosphere∗ 5.1  Introduction  Although heterogeneous ice nucleation occurring in the atmosphere has been studied for many years its effect on climate remains poorly understood. A better description of the heterogeneous nucleation process leading to the formation of ice clouds is certainly necessary to help predicting precipitation patterns as well as in estimating the amount of solar radiation reaching the surface of the Earth [1–6]. An important and abundant type of aerosol particles that influence ice cloud formation in the atmosphere are mineral dust particles [7]. A list of the common clays minerals found in aerosolized dust samples include quartz, illite, muscovite, chlorite, kaolinite, and calcite [8]. Our focus here is on kaolinite particles. These particles are known to be effective ice nuclei (IN) based on laboratory studies [9–11] and they make up a significant amount of the aerosolized mass of mineral dust [8] (about 5 − 10 %). The reason why kaolinite particles are good IN remains unclear. One possible hypothesis is that there exists a good crystallographic match between the kaolinite and hexagonal ice surfaces [12]. The kaolinite surface was thought to arrange the adsorbed A version of this chapter will be submitted for publication. T. Croteau, A. K. Bertram, G. N. Patey, Structural Analysis of Water Adsorbed on Kaolinite Particles Containing Trenches at 235 K and Possible Connections with Ice Nucleation in the Atmosphere ∗  84  5.1. Introduction water, according to specific sites, in a structure that has a lattice constant similar to ice. These structures would then be good templates for the formation of ice embryos. But we have recently shown that this is likely not the case [13]. In fact, the lattice constants of the ring structures formed on kaolinite surfaces do not match those of hexagonal ice. This results in significant strain of the ice embryo which does not favour the growth of ice. These results are supported by ab initio studies [14, 15]. The simulation studies [13–15] support the idea that surface defects or irregularities such as steps, cracks, trenches, etc. must play a major role in ice nucleation by kaolinite particles, and likely other types of clay particles. Closely related, ice freezing experiments also indicate that ice nucleation on mineral dust surfaces occurs at specific types of surface defects, so called active sites [12, 16–19]. For example, Marcolli et al. [16] showed that active sites (about 10 ˚ A2 in size) were needed to describe freezing experiments on Arizona test dust, a type of mineral dust. In a recent study [20] we investigated water adsorption in kaolinite defects at 298 K. Specifically we considered nanoscale trenches with widths ranging from 14.78 ˚ A to 73.91 ˚ A and a depth of 23.19 ˚ A. The motivation for this geometry and these length scales is discussed below. This work showed that at 298 K water completely filled these trenches at subsaturation RH values. In the following, we continue with these studies looking at the structure of water in the same trenches at 235 K. One of the goals of this study is to determine if water in these trenches has order, which is imposed by a local field and/or confinement. Another goal is to assess whether these types of defects might be important in ice nucleation in the atmosphere. These simulations could also point the way for future experiments and simulations, that may lead to a better understanding of heterogeneous ice nucleation on kaolinite, and other clay surfaces. Related to our work, over the last decade or so, several simulations on ice formation in water have been reported (for examples, see Refs. [21–35]). In some cases, ice formation was achieved by confining the water or by applying a field. This has the effect of  85  5.2. Model and Method reducing the number of available disordered hydrogen-bond network structures compared to an unconstrained, three dimensional system, thus enhancing ice nucleation. As an example, in their molecular dynamics simulations, Svishchev and Kusalik [21] were able to crystallize liquid water into some amorphous forms of ice, as well as into cubic ice I, depending on the density, by applying a field at a temperature of 250 K. By confining the water in a hydrophobic slit pore less than one nanometre wide Koga et al. [22–24] were able to produce a liquid to an amorphous bilayer transition. In other similar cases, ice was observed in nanotubes [25, 26], pores [27] and slits [26, 28]. In the following, we analyze the water structures formed using density and order parameter profiles, molecular imaging and hydrogen bond numbers. The results are based on grand canonical Monte Carlo simulations performed at 235 K using both the SPC/E and TIP5P-E water models. The remainder of the paper is organized in three parts. The models and simulation issues are considered briefly in Section II, results are presented and discussed in Section III, and our conclusions are summarized in Section IV.  5.2  Model and Method  Grand canonical Monte Carlo simulations of water adsorption were performed on kaolinite surfaces containing trenches at 235 K. The trench systems used are based on the schematic displayed in Fig. 4.1. The motivation for this particular defect type comes from scanning electron microscope (SEM) images of kaolinite particles. Such images show that the basal planes of kaolinite are relatively smooth (at least on the length-scale of the resolution of an SEM) but the surfaces perpendicular to the basal planes are rough and irregular, with features that often look like trenches or cracks that should have sides (i.e., walls) that correspond to the Al- and Si-surfaces. The trench bottoms and the surface areas between the trenches consist of unprotonated edges (100 plane). In the SEM images the dimensions of the trenches and cracks are larger than the dimensions used in our simulations, but it seems reasonable to assume that similar structures are present at smaller length-scales 86  5.3. Results and Discussion below the resolution of SEM. The same four trenches as in our previous study [20], labeled 1, 2, 3 and 4, were used (see Fig. 4.2). The complete description of the simulation setup and construction of these trenches has already been discussed in Ref. [20] (Chapter 4). Let us recall that the trenches in Fig. 4.2 are made of unprotonated edges (bottom of the trench and top surface), Al- and Si-surfaces (trench walls). They have a trench depth of 23.19 ˚ A and a y axis of 35.77 ˚ A although it is made infinite by periodically replicating the system in the xy plane. The trench width varies from 14.78 ˚ A (trenches 1A, 1B, and 2) to 44.34 ˚ A (trench 3) and 73.91 ˚ A (trench 4). A complete list of the parameters is given in Table 4.1. Differing from Ref. [20], we set the distance separating the top and bottom slabs in our double slab systems (Fig. 4.1) to 100 ˚ A instead of 50 ˚ A to ensure that the water molecules are not interacting with the other opposing surface (which, in turn, is used to cancel the effect of lattice truncation). This allows the formation of large mounds above the top unprotonated edge surface, without significant interactions with the other opposing surface. Also, we used two different well-known water models, namely the SPC/E [36] and TIP5P-E [37], to ensure as far as possible, that our results and conclusions don’t depend on the model employed. The force field parameters for the surface are those of CLAYFF [38], a model developed for clay particles. This model was also used in our previous studies of water adsorption using kaolinite [13, 20].  5.3 5.3.1  Results and Discussion Water Adsorption at 235 K  The results we report were obtained at chemical potentials of −54 kJ/mol and −49.5 kJ/mol for all trenches studied using the SPC/E and TIP5P-E models, respectively. Since chemical potentials for the bulk transitions of SPC/E and TIP5P-E water do not appear to have been obtained at 235 K, we use the estimates given by our simulations of atomistically  87  5.3. Results and Discussion flat kaolinite surfaces [13] to define relative humidities (RH). These are estimated as ∼ 60 % and ∼ 77 % RH for SPC/E and TIP5P-E, respectively. Note that in all likelyhood these values are upper limits on the RH, because in our simulations [13] the surface-water interactions are expected to slightly favor the condensed phase. Images of the different coverages found in our systems, at the above RH, are displayed in Figs. 5.1 and 5.2. As was observed in Ref. [20], the coverage found in trench 1B is larger than for its counterpart, trench 1A, which has the inter-trench spacing empty (see Fig. 5.1). But as will be seen later, this has no qualitative effect on the structures formed. The same observation applies for the different water models. Although smaller coverages are found for the narrow trenches using TIP5P-E compared to SPC/E (see Fig. 5.2a and 5.2b), as will be seen below, the structural results are very similar whether we use the SPC/E or TIP5P-E models. As shown in Figs. 5.1 and 5.2 water is strongly taken up at the indicated RH. All trenches are filled except for Trench 4 (Figs. 5.2e and f). It is only a matter of convergence, which is slow at this temperature (this system is still uptaking water molecules), for the trench to fill completely. Note that the mounds covering trench 3 (Figs. 5.2c and d) are also still growing after ∼ 2 × 109 Monte Carlo steps. Based on our 298 K results [20], all trenches will fill at the specified relative humidity given sufficient time. Nevertheless, the partially filled Trench 4 is instructive for assessing the structure of water in the trenches (see below). Also note that, since the results for the narrow trenches 1A, 1B, and 2 are very similar, we show results only for trench 1B which, in turn, gives an “upper limit” for the adsorption, and also for structural order, as will be seen below.  5.3.2  Water Density Profiles  We now need to determine if the structure adopted by the water molecules in our systems is similar to any recognizable form of ice, or at if it least corresponds to some solid state. The first calculation we performed to characterize the water structure are water density  88  5.3. Results and Discussion  (a)  (b)  Figure 5.1: Snapshots of trenches 1A (a) and 1B (b) using SPC/E at 235 K and ∼ 60 % RH. The water oxygen atoms are blue, lattice oxygen atoms are red, hydrogen atoms are white, ”lone pairs” are orange, aluminum atoms are grey, and silicon atoms are brown.  89  5.3. Results and Discussion  (a)  (b)  (c)  (d)  (e)  (f )  Figure 5.2: Snapshots of trench 1B (a,b), 3 (c,d) and 4 (e,f) using SPC/E (a,c,e, ∼ 60 % RH) and TIP5P-E (b,d,f, ∼ 77 % RH) at 235 K. The atoms are denoted as in Fig. 5.1. 90  5.3. Results and Discussion profiles in all three dimensions. As can be seen in Figs. 5.3, 5.4, and 5.5, very clear structural evidence (high intensity peaks) suggesting a solid state can be extracted from these calculations, especially for the narrowest trenches. Firstly, Fig. 5.3 shows that inside trenches that are 14.78 ˚ A wide (i.e. trench 1B), the water molecules can be separated into six well-defined layers along the x-axis, one next to both the Al- and Si-surfaces, and four distinct layers in between. The high intensity of the peaks shows that the narrow trenches are very dense. In fact, trench 1B is the densest of all trenches (∼ 1.6g/cm3 ) due to its greater attractive water-surface interactions (see Ref. [20]). The peaks near the walls, where the interaction is strongest, represent denser layers. As one moves away from the walls, the intensities of the peaks decreases. The fact that the water molecules are confined between two walls certainly imposes more order, at least near the walls. This effect clearly decreases as the trench width is increased leading to less ordered structures away from the walls for trenches 3 and 4 (see Fig. 5.3). Note from Fig. 5.3, that the results are not strongly dependent on the water model we use. The density distributions calculated along the y axis show high intensity peaks throughout the trench for trench 1B (see Fig. 5.4), which suggests that the water molecules are in a solid state as this axis is not constrained by any surface. However, the absence of confinement along this axis leads to structural and layering patterns that are not as well defined. For this reason no images of the layering in this direction are shown below. The z direction, on the other hand, is constrained from the bottom by an unprotonated edge surface leading to a density distribution where well-defined patterns can be identified. Along this axis the water molecules are most ordered for the narrow trenches (see Fig. 5.5). The first peak corresponding to the water layer interacting immediately with the unprotonated edge surface has the highest intensity, and is followed by a series of sharp peaks corresponding to successive layers. Again, as the trench width is increased (trenches 3 and 4), these secondary peaks become less intense (Fig. 5.5).  91  5.3. Results and Discussion  8 6 4 2 0 3  ρ (g/cm )  5  0  2  4  6  8  10  12  14  16  6  12  18  24  30  36  42  48  4 3 2 1 0 0 5 4 3 2 1 0  0  10  20  30  40  50  60  70  80  x (Å) Figure 5.3: The water density (g/cm3 ) profile obtained at 235 K along the x axis (trench width) for trenches 1B (top panel), 3 (middle panel) and 4 (bottom panel) using the SPC/E (solid black curves) and TIP5P-E (dashed red curves) water models at ∼ 60 % and ∼ 77 % RH, respectively. Note the different scales on the x axis for all three trenches.  92  5.3. Results and Discussion  3 2  3  ρ (g/cm )  1  1.6 1.2 0.8  1.6 1.2 0.8 0  2  4  6  8  10  12  14  16  18  20  y (Å) Figure 5.4: The water density (g/cm3 ) profile obtained at 235 K along the y axis for trenches 1B, 3, and 4 (from top to bottom) using the SPC/E (solid black curves) and TIP5P-E (dashed red curves) water models at ∼ 60 % and ∼ 77 % RH, respectively.  93  5.3. Results and Discussion  5 4 3 2  3  3  ρ (g/cm )  1  2 1  3 2 1 0  0  2  4  6  8  10  12  14  16  18  20  z (Å) Figure 5.5: The water density (g/cm3 ) profile obtained at 235 K along the z axis (trench depth) for trenches 1B, 3, and 4 (from top to bottom) using the SPC/E (solid black curves) and TIP5P-E (dashed red curves) water models, at ∼ 60 % and ∼ 77 % RH, respectively.  94  5.3. Results and Discussion  5.3.3  Structural Imaging  In this section we focus only on trench 1B, which has the best defined and most ordered structure of all systems. Also, for better clarity, we only show images for the SPC/E model, which has fewer sites, since the results are very similar for both water models. As mentioned above, the narrow trenches have six well defined layers along the trench width (x coordinate confined between both trench walls). These are clearly represented in Fig. 5.6 where all six layers have been isolated for trench 1B. The peaks near the walls, where the interaction is strongest, represent denser and more ordered layers. Based on the above density profiles and images, the water structure at the surfaces appears to be more ordered. As observed in our previous study [13], the surface-induced ring structures (see Fig. 5.6d) do not resemble ice, at least not ice Ih, for which calculation of a structure factor on these layers does not reveal any specific order. Nevertheless, they do show a great deal of order, the water molecules are well aligned and the hydrogen atoms all appear to point in the same general direction (z). The layer represented in Fig. 5.6c is next to the Si-surface thus explaining why the majority of the water molecules have their hydrogens lying flat, preferring to interact with other water molecules rather than with the surface, which is considered hydrophobic [13] (only weakly interacting). On the other side of the trench, the layer next to the Al-surface shows a greater tendency to interact with the clay atoms as more hydrogen atoms are oriented towards it. We note some hexagonal patterns in Fig. 5.6d, but these match the surface and not hexagonal ice as discussed in Ref. [13], where a misfit strain of 14 % was calculated between the hexagonal patterns formed on the Al-surface and ice Ih. When looking at the separated layers along the z axis as shown in Fig. 5.7a, a very interesting observation can be made. All water molecules have their hydrogen atoms pointing upward from the first layer up. Figs. 5.7c and d show this even more clearly where the first layer is projected flat in the xy plane both from the top and bottom, showing networks of white hydrogen atoms and blue water oxygens, respectively. This  95  5.3. Results and Discussion structure is retained all the way up even outside the trench (see Fig. 5.8), although the ordering starts to decrease at this location (see below), where a qualitative difference is observed between the orientation of the water molecules above the trench and those binding to the top unprotonated edge. It is interesting to note, that different patterns are found on the different unprotonated edges depending on the location. The same effect is observed for the larger surfaces although it is not as pronounced. The structure found on the top unprotonated edge surface is similar to that observed in our study of flat kaolinite surfaces, and is, therefore, what one would expect for an infinitely wide trench.  5.3.4  Orientational Ordering  As a way of quantifying the alignment of the water molecules in a trench, we calculated a dipole order parameter profile in the z direction. For a given configuration this is defined as 1 m(z) = N(z)µw  N (z)  µi ,  (5.1)  i=1  where N(z) is the number of water molecules included, as a function of z, in a rectangular volume element of thickness 0.116 ˚ A, µw is the dipole moment of either an SPC/E or TIP5P-E water molecule, and the vector sum is over all water dipoles included in the volume element. The results reported below were obtained as averages over 2000 configurations. As a reference, if all water molecules were to be perfectly aligned (perfect ferroelectric order), then the magnitude of the average dipole order parameter would be < |m(z)| >= 1. On the other hand, a value of < |m(z)| >≈ 0 is expected if the dipoles are disordered. The calculations were done with NVT runs using configurations corresponding to ∼ 60 % and ∼ 77 % RH at 235 K for SPC/E and TIP5P-E, respectively. The average components < mx (z) >, < my (z) >, and < mz (z) > of m(z) are displayed in Fig. 5.9 for trench 1B, and show clear alignment of the water molecules inside the trench for both water models. The curves shown in Fig. 5.9 confirm our observations from simple imaging of the different layers in the z axis, namely, that the net moment 96  5.3. Results and Discussion  8 7  1  3  ρ (g/cm )  6  2  3  4  5  6  3  5 4 3 2 1 0  0  2  4  6  8  10  12  14  16  x (Å)  z y x  (a)  (b)  1  6  (c)  (d)  Figure 5.6: Images of the individual layers along the x axis for trench 1B using SPC/E at 235 K and ∼ 60 % RH: (a) all layers along the x axis, (b), (c) and (d) are flat views of the left side for the first, third, and sixth layers as seen in (a). An artificial spacing of 5 ˚ A was inserted between the layers in (a). The water oxygen atoms are blue and hydrogen atoms are white. The yellow hexagonal pattern overlayed in (d) matches the surface but not ice Ih (see Chapter 2).  97  5.3. Results and Discussion  14  6  12  5  8  3  6  2  4  z (Å)  10  4  3  2  1 4  3  2  ρ (g/cm )  1  0  0  3  z y x  (a)  1  (b)  1  (c)  (d)  Figure 5.7: Images of the individual layers along the z axis for trench 1B using SPC/E at 235 K and ∼ 60 % RH: (a) first six layers along the z axis starting from the bottom of the trench, (b) and (c) are top views of the first and third layers, respectively, and (d) is a bottom view of the first layer. An artificial spacing of 5 ˚ A was inserted between the layers in (a). The atoms are denoted as in Fig. 5.6.  98  5.3. Results and Discussion  (a)  (b)  (c)  Figure 5.8: Top (b) and bottom (c) views of the first layer (black rectangle in (a)) above the top unprotonated edge at 235 K and ∼ 60 % RH using SPC/E for trench 1B. The atoms are denoted as in Fig. 5.6. is essentially pointing upward along the z axis. However, at around 24 ˚ A above the bottom of the trench, which matches the distance where the water molecules leave the trench, the alignment of the water molecules starts to loose some of its ordering in the z component. After that, the water molecules become less ordered as you move away from the trench. The other two components (< mx (z) > and < my (z) >) are much less significant, particularly for TIP5P-E, where they oscillate about zero for both the x and y components. Also note, that the x and y components are unaffected by the change in structure initiated upon exiting the trench. The average dipole order parameter profiles displayed in Fig. 5.10 for trenches 1B, 3 and 4 emphasize the ferroelectric character of the water structures found inside the trenches for both water models. The magnitude is slightly below 1 inside trench 1B. As mentioned before, trenches 3 and 4 are less ordered. Nevertheless, the curves displayed in Fig. 5.10 show that there is still significant ordering inside the larger trenches. Note again, that most of the ordering comes from the z component. This is why the alignment of the water molecules starts to decrease at around 24 ˚ A above the bottom of the trench  99  5.3. Results and Discussion for all trenches, as discussed above. In addition to the average, and average vector components, of the dipole order parameter profiles displayed in Figs. 5.9 and 5.10, we report in Table 5.1 the average dipole order parameter numbers for all trenches studied, both for all water molecules located below and above the top unprotonated edge surface. This location serves as a separation since there is a marked change in structure at this interface between the water molecules located above the trench and above the top unprotonated edge surfaces. Both Figs. 5.9 and 5.10, as well as Table 5.1, clearly show that the water molecules inside the trench are highly ordered. The average dipole order parameter values including all the water molecules inside the trench are as high as 0.94 and 0.96 for trench 1B, using the SPC/E and TIP5P-E water models, respectively (see Table 5.1). Having the inter-trench spacing empty leads to slightly smaller values of ∼ 0.9 for trenches 1A and 2. The values for the larger trenches 3 and 4, ∼ 0.7 and ∼ 0.6, respectively, show that as the width of the trench is increased the ordering of water molecules decreases. Nevertheless, all these values demonstrate significant orientational order, and indicate that the results don’t depend on the water model. Outside the trench, the values are much smaller (between ∼ 0.2 and ∼ 0.4). These values agree with Figs. 5.9 and 5.10 which show a clear decrease of the dipole order parameter as you leave the trench, and with Fig. 5.8 which shows clear images of the different sections with differing orientation of the water molecules above the top unprotonated edge. The presence of an ordering field acting inside the trench was verified by the calculation of the dipole order parameter profile of a single water molecule along both walls (Al- and Si-surface) as well as in the middle of trench 1B. The results presented in Fig. 5.11 show, that for the majority of the water molecules located away from the walls, there is a field orienting the water upwards inside the trench (∼ 23˚ A), and it extends some 15 ˚ A above the top of the trench. Beyond that distance, the field drops significantly as indicated by the decrease in the z component of the dipole order parameter. Near both walls, the field  100  5.3. Results and Discussion  Table 5.1: Calculated dipole order parameters (obtained for the water molecules included either inside or outside the trench) and hydrogen bond numbers (with an angle condition of either 35o or 40o ). The associated standard deviations are included in brackets. Trench Order parameter Type SPC/E TIP5P-E inside outside inside outside 1A 0.907 (0.001) 0.31 (0.01) 1B 0.9463 (0.0008) 0.29 (0.01) 0.9560 (0.0003) 0.39 (0.01) 2 0.886 (0.001) 0.20 (0.01) 0.892 (0.002) 0.22 (0.01) 3 0.721 (0.002) 0.26 (0.01) 0.737 (0.003) 0.29 (0.01) 4 0.596 (0.003) 0.251 (0.008) 0.600 (0.002) 0.227 (0.004) Trench Hydrogen bond number Type SPC/E 35o 40o 1A 3.94 (0.01) 4.16 (0.01) 1B 4.01 (0.01) 4.310 (0.009) 2 3.91 (0.01) 4.14 (0.02) 3 3.724 (0.007) 3.832 (0.008) 4 3.632 (0.007) 3.727 (0.007) mainly orients the water molecules in the z direction, but it can also have a significant contribution from the x component. This has the effect of tilting the water molecules slightly sideways at some locations along the walls. Along the Al-surface, the effect of the x and z components of the field also extend some 15 ˚ A above the top of the trench.  5.3.5  Hydrogen Bonding  When calculating the number of water-water hydrogen bonds we only used the water molecules away from the walls and the same conditions as in our previous studies were used for the SPC/E model [13, 20]. However, it was found that the results were quite sensitive to the angle condition. For example, the average number of hydrogen bonds inside trench 1B was calculated as 4.01 and 4.31 using angle conditions of 35o and 40o for SPC/E, respectively (see Table 5.1). It is quite surprising that such a small change in angle condition led to such a big modification in the estimated number of hydrogen bonds. This is probably due to the high density of the system. Therefore, care must be  101  <mz(z)>  5.3. Results and Discussion  1 0.8 0.6 0.4 0.2  <my(z)>  0.4 0.2 0 -0.2 -0.4  <mx(z)>  0.4 0.2 0 -0.2 -0.4 0  4  8  12  16  20  24  28  32  z (Å) Figure 5.9: Profiles of the average dipole order parameter vector components obtained along the z axis starting from the bottom of trench 1B at 235 K using SPC/E at ∼ 60 % RH (black curves) and TIP5P-E at ∼ 77 % RH (red curves).  102  5.3. Results and Discussion  1 0.8  <|m(z)|>  0.6 0.4 0.2 1 0.8 0.6 0.4 0.2  0  4  8  12  16  20  24  28  32  z (Å) Figure 5.10: Profiles of the average magnitude of the dipole order parameter obtained along the z axis starting from the bottom of the trench at 235 K for trenches 1B (red curve), 3 (black curve), and 4 (blue curve) using SPC/E (bottom panel) and TIP5P-E (top panel) at ∼ 60 % and ∼ 77 % RH, respectively.  103  5.3. Results and Discussion  0.8 0.4 0 -0.4 -0.8  <mα(z)>  0.8  0.4 0  -0.4 -0.8 0.8 0.4 0 -0.4 -0.8 0  10  20  30  40  50  60  70  z (Å) Figure 5.11: Profiles of the average dipole order parameter vector components obtained along the z axis for 1 water molecule near the Si-surface (bottom panel), middle of the trench (middle panel) and near the Al-surface (top panel) for trench 1B at 235 K using the SPC/E water model. The black, red and blue curves represent the x, y and z components of the order parameter, respectively.  104  5.4. Relationship between Observations and Ice Nucleation in the Atmosphere taken when applying this calculation to a dense system. Nevertheless, the values shown in Table 5.1 indicate that the narrow trenches 1 and 2 are the most hydrogen bonded of the systems studied with values above 3.9. Trench 1B, the densest system, has the largest value. Increasing the trench width leads to a decrease in the number of hydrogen bonds even if the systems are still quite dense. However, these larger systems (trenches 3 and 4) are denser near the walls (sides and bottom of the trench). Since we only used the water molecules away from the walls to calculate the average hydrogen bond number, it is no surprise to obtain smaller values for the larger systems, as they are less dense in that region compared to the narrow trenches. In fact, trench 4 is not completely filled explaining why it has a lower value. As a comparison, the regions used to calculate the hydrogen bond numbers for trenches 1B and 3 have densities of ∼ 1.6g/cm3 and ∼ 1.1g/cm3, respectively.  5.4  Possible Relationship between Observations and Ice Nucleation in the Atmosphere  The above results provide hints that the trench-like defects studied may be important for ice nucleation, although we cannot be certain whether or not these surfaces will be favorable for bulk ice nucleation and growth. First, the hydrogen bond numbers are close to four for the narrow trenches, which is closer to ice than to liquid water. Second, the highly ordered water structures determined by our density profiles suggest that a dense solid state is present in the trenches. One might expect these environments be more favorable to bulk ice nucleation and growth than a liquid water environment, or atomically flat kaolinite surfaces, that have a large misfit with bulk ice. Molecular dynamics calculations that investigate the possible growth of bulk ice in these environments might prove interesting, since this technique has had some success in observing the growth of ice structures [21–24, 26, 28, 29]. Also, the dipole order parameter calculations indicate the presence of a proton ordered  105  5.4. Relationship between Observations and Ice Nucleation in the Atmosphere ferroelectric structure in the trenches. Usually, ice phases are proton disordered and there is no ferroelectric order. However, some proton ordered bulk ice structures are known to exist and can be ferroelectric. For example, ice II is proton ordered and has no proton disordered analogue, ice IX is the proton ordered version of ice III (ice III is itself partially ordered), and ice VIII is the proton ordered version of ice VII [39]. Experiments also show that proton rearrangement in normal ice Ih can be enhanced by doping the crystals with alkali hydroxides [40–45], by applying an electric field [46–49] or, by using a surface which promotes the growth of an ordered ice phase [50, 51]. Perhaps the structured water in the trenches would be favorable for nucleating and growing proton ordered, ferroelectric ice, which could than be converted into the thermodynamically favorable hexagonal ice Ih further from the surface. It would be very interesting if a series of experiments could be performed, to determine if trenches are really an hospitable environment for the formation of ferroelectric types of ice. It can also be speculated that adding mobile ions to our trench systems would greatly influence the results, as these would tend to counteract the ordering field. This would be consistent with recent laboratory studies [52], that show a strong effect of ions on the ice nucleation properties of kaolinite particles, as well as for other types of minerals. Simulations and experiments would also be needed in this case to establish the possible effects and implications of the presence of ions.  106  Bibliography [1] Kaufman, Y. J.; Tanre, D.; Boucher, O. Nature 2002, 419, 215. [2] Murray, B. J.; Knopf, D. A.; Bertram, A. K. Nature 2002, 434, 202. [3] Sastry, S. Nature 2005, 438, 746. [4] Shaw, R. A.; Durant, A. J.; Mi, Y. J. Phys. Chem. B 2005, 109, 9865. [5] Sassen, K. Nature 2005, 434, 456. [6] Forster, P.; Ramaswamy, V.; Artaxo, P.; Berntsen, T.; Betts, R.; Fahey, D. W.; Haywood, J.; Lean, J.; Lowe, D. C.; Myhre, G.; Nganga, J.; Prinn, R.; Raga, G.; Schulz, M.; Dorland, R. V. Changes in Atmospheric Constituents and in Radiative Forcing., in Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change 2007, edited by S. Solomon, D. Qin, M. Manning, Z. Chen, M. Marquis, K.B. Averyt, M. Tignor, and H.L. Miller, pp. 129-234, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA. [7] Cantrell, W.; Heymsfield A. J. Bull. Amer. Meteor. Soc. 2005, 86(6), 795. DeMott, P. J.; Cziczo, D. J.; Prenni, A. J.; Murphy, D. M.; Kreidenweis, S. M.; Thomson, D. S.; Borys, R.; Rogers, D. C. Proc. Nat. Acad. Sci. 2003b, 100, 25, 14655. [8] Glaccum, R. A.; Prospero, J. M. Mar. Geol. 1980, 37(3-4), 295. [9] Zimmerman, F.; Ebert, M.; Worringen, A.; Schutz, L.; Weinbruch, S. Atmos. Environ. 2007, 41, 8219. [10] Zuberi, B.; Bertram, A.K.; Cassa, C. A.; Molina, L. T.; Molina, M. J. Geophys. Res. Lett. 2002, 29, 142. [11] Dymarska, M; Murray, B. J.; Sun, L; Eastwood, M.; Knopf, D. A.; Bertram, A. K. J. Geophys. Res. 2006, 111, D04204, doi:10.1029/2005JD006627. [12] Pruppacher, H. R.; Klett, J. D. Microphysics of Clouds and Precipitation, 2nd rev. ed.; Kluwer Academic: Dordrecht, The Netherlands, 1997; pp 329-331. [13] Croteau, T.; Bertram, A. K.; Patey, G. N. J. Phys. Chem. A 2008, 112 (43), 10708. Croteau, T.; Bertram, A. K.; Patey, G. N. J. Phys. Chem. A 2009, 113 (27), 7826. [14] Hu, X. L; Michaelides, A. Surf. Sci. 2007, 601, 5378. [15] Hu, X. L; Michaelides, A. Surf. Sci. 2008, 602, 960. 107  Chapter 5. Bibliography [16] Marcolli, C.; Gedamke, S.; Peter, T.; Zobrist B. Atmos. Chem. Phys. 2007, 7, 5081. [17] Archuleta, C. M.; DeMott, P. J.; Kreidenweis, S. M. Atmos. Chem. Phys. 2005, 5, 2617. [18] Hung, H.-M.; Malinowski, A.; Martin S. T. J. Phys. Chem. A 2003, 107, 1296. [19] Welti, A.; Luond, F.; Stetzer, O; Lohmann, U. Atmos. Chem. Phys. Discuss. 2009, 9, 6929. [20] Croteau, T.; Bertram, A. K.; Patey, G. N. J. Phys. Chem. A (accepted) 2009, Included here as Chapter 4. [21] Svishchev, I. M.; Kusalik, P. G. J. Am. Chem. Soc. 1996, 118, 649. [22] Koga, K.; Tanaka, H.; Zeng, X. C. Nature 2000, 408, 564. [23] Slovak, J.; Koga, K.; Tanaka, H.; Zeng, X. C. Phys. Rev. E 1999, 60(5), 5833. [24] Slovak, J.; Tanaka, H.; Koga, K.; Zeng, X. C. Pysica A 2001, 292, 87. [25] Sliwinska-Bartkowiak, M.; Jazdzewska, M.; Huang, L. L.; Gubbins, K. E. Phys. Chem. Chem. Phys. 2008, 10, 4909. [26] Tanaka, H; Koga, K. Bull. Chem. Soc. Jpn. 2006, 79(11), 1621. [27] Seyed-Yazdi, J.; Farman, H.; Dore, J. C.; Webber, J. Beau W.; Findenegg, G. H. J. Phys.: Condens. Matter 2008, 20, 205108. [28] Hirunsit, P.; Balbuena, P. B. J. Phys. Chem. C 2007, 111, 1709. [29] Matsumoto, M.; Saito, S.; Ohmine, I. Nature 2002, 416, 409. [30] Vrbka, L.; Jungwirth, P. J. Phys. Chem. B 2006, 110, 18126. [31] Radhakrishnan, R.; Trout B. L. Phys. Rev. Lett. 2003, 90(15), 158301. [32] Radhakrishnan, R.; Trout B. L. J. Am. Chem. Soc. 2003, 125, 7743. [33] Brukhno, A. V.; Anwar, J.; Davidchack, R.; Handel, R. J. Phys.: Condens. Matter 2008, 20, 494243. [34] Quigley, D.; Rodger, P. M. J. Chem. Phys. 2008, 128, 154518. [35] Bauerecker, S.; Ulbig, P.; Buch, V.; Vrbka, L.; Jungwirth, P. J. Phys. Chem. C 2008, 112, 7631. [36] Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. [37] Lisal, M.; Nezbeda I.; Smith, W. R. J. Phys. Chem. B 2004, 108, 7412. [38] Cygan, R. T.; Liang, J.-J.; Kalinichev, A. G. J. Phys. Chem. B 2004, 108, 1255. 108  Chapter 5. Bibliography [39] Malenkov, G. J. Phys.: Condens. Matter 2009, 21, 283101. [40] Worz, O.; Cole, R. H. J. Chem. Phys. 1969, 51, 1546. [41] Kawada, S.; J. Niinuma, J. J. Phys. Soc. Jpn. 1977, 43, 715. [42] Johari, P.; Whalley, E. J. Chem. Phys. 1981, 75, 1333.5 [43] Tajima, Y.; Matsuo, T.; Suga, H. J. Phys. Chem. Solids 1984, 35, 1135. [44] Yamamuro, O.; Kuratomi, N.; Matsuo, T.; Suga, H. Solid State Commun. 1990, 73, 317. [45] Tyagi, M.; Murthy, S. S. N. J. Phys. Chem. A 2002, 106, 5072. [46] Dengel, O.; Eckener, U.; Plitz, H.; Riehl, N. Phys. Lett. 1964, 9, 291. [47] Johari, G. P.; Jones, S. J. J. Chem. Phys. 1975, 62, 4213. [48] Jackson, M.; Whitworth, R. W. J. Chem. Phys. 1995, 103, 7647. [49] Jackson, M.; Whitworth, R. W. J. Phys. Chem. B 1997, 101, 6177. [50] Su, X.; Lianos, L.; Shen, Y. R.; Somorjai, G. Phys. Rev. Lett. 1998, 80, 1533. [51] Iedema, M. J.; Dresser, M. J.; Doering, D. L.; Rowland, J. B.; Hess, W. P.; Tsekourus, A. A.; Cowin, J. P. J. Phys. Chem. B 1998, 102, 9203. [52] Eastwood, M. L.; Cremel, S.; Gehrke, C.; Girard, E.; Bertram, A. K. J. Geophys. Res.-Atm. 2008, 113.  109  Chapter 6 Conclusions and Perspective The simulations presented in this thesis are all related to the adsorption and structure of water on kaolinite surfaces with and without trenches. In Chapters 2 and 3, it was demonstrated that only slightly overgrown monolayers can be adsorbed on atomistically flat kaolinite surfaces, in the atmospherically relevant range of relative humidity at 235 and 298 K. The structures adopted by the water molecules on the different surfaces do not support the commonly held idea, that ice nucleation is occurring due to a good match with the surface [1]. In fact, as verified by ab initio studies [2, 3], the hexagonal strucuture formed on the pseudo-hexagonal Al-surface does not match ice Ih. Given that the results obtained using flat surfaces were unable to explain the large water adsorption reported in some experimental studies [4], we turned our attention towards another important factor in water adsorption, and possible ice nucleation, the presence of active sites in the form of surface defects. The defects simulated are trenches as presented in Chapters 4 and 5. In Chapter 4, the use of trenches of different sizes revealed that large amounts (apparent multi-layers) of water could be adsorbed on kaolinite particles. These results allowed us to make a direct comparison with the recent experimental studies, that showed multilayer coverages on kaolinite and other clay particles [4]. Our last study consists of an analysis of the water structure formed in trenches at 235 K. Structural evidence shows that the water adsorbed inside the trenches is a dense ferroelectric solid, with a geometrical layered structure in all three dimensions, and alignment of the water molecules. However, it does not closely resemble any known bulk ice structure. This research is at the junction of physical and atmospheric chemistry. The simulations performed give us a better understanding of the interactions of water with clay surfaces, which is a key factor in understanding the formation of ice on aerosol particles found in 110  Chapter 6. Conclusions and Perspective the atmosphere. The use of realistic models for both the clay surface and water molecules gives us significant confidence in the results obtained. Also, the fact that two different water models give very similar results is an indication of the reliability of the simulations. The results obtained on flat surfaces agreed with other simulation and ab initio calculations [2, 3, 5–8], which further support the theoretical foundation of our approach, and gives us confidence in the comparison of our results with experiments. The simulation studies reported in the literature were mostly performed for closed systems. The importance of using an open system in the understanding of the naturally occurring water coverage on kaolinite, which is related to the strength of the water-surface interaction and the relative humidity, is emphasized by our results. A novelty of this thesis is the inclusion of edges and trenches in the description of water adsorption on kaolinite. As it turns out, edges may play a significant role in heterogeneous chemistry, where water is needed, because these surfaces have an enormous water affinity. On the other hand, the study of trenches revealed the presence of large water “pools”. This is a very important result as ice cannot form without the presence of a significant amount of water. A very important conclusion of this thesis is that active sites such as cracks or trenches are a necessary condition for the presence of large quantities of water, and possibly for the formation of ice. As seen in Chapters 4 and 5, multilayer coverage comparable to experimental results are only found in trenches, not on flat kaolinite surfaces. However, this observation needs to be confirmed by performing water adsorption simulations on other types of particles (i.e., other clay minerals) containing active sites in the form of defects. Then, the results could be compared to a set of water adsorption or ice nucleation experiments performed on the same particles. The importance of surface roughness should also be investigated experimentally by choosing a set of particles which are smooth (as free as possible of defects) and others which are not. Atomic force microscopy (AFM) could be used to determine the surface roughness.  111  Chapter 6. Conclusions and Perspective An initial objective of this thesis was to explain the difference between the ice nucleation ability of various aerosols. The hypothesis of a good surface match between the water molecules and the kaolinite surface was first ruled out, as the hexagonal arrangements found on the Al-surface of kaolinite did not have the same size as an ice structure. Also, the strong interactions between the first water layer and the surface make it “hydrophobic” and prohibit any further growth. The next idea was that active sites must play a role in the ice nucleation ability of particles. This second argument was supported by our simulations of trenches where large amounts of water were found. However, as mentioned before, the formation of a recognizable bulk ice structure was not achieved, but the structured water found in the trenches may be important in hexagonal ice nucleation and growth. A great deal of research still has to come before we can finally understand precisely the effect of aerosols on climate. Researchers will now have to study, amongst other things, the growth of ice embryo under different environments (i.e., different types of cracks, socalled active sites) and conditions on different aerosol particles, both experimentally and with the help of computer simulations. Our simulations point to the idea that the initial ice embryo arising from defects might be ferroelectric. It would be very helpful if a set of experiments could be designed to try and determine if some kind of ferroelectric ice is observed at the initial stages of nucleation and growth. Also, more computer simulations aimed at the growth of bulk ice in different environments and under different conditions, are necessary. The inclusion of mobile charged particles is also of interest, as these may well neutralize the polarizing effect of the field, which aligns the water molecules in our trenches. This could help explain why the nucleation ability of aerosol particles decreases after they are treated with salts or acids as observed in experiments [9].  112  Bibliography [1] Pruppacher, H. R.; Klett, J. D. Microphysics of Clouds and Precipitation, edited by Kluwer Academic, Dordrecht, 2nd rev. ed., 1997. [2] Hu, X. L.; Michaelides, A. Surf. Sci. 2007, 601, 5378. [3] Hu, X. L.; Michaelides, A. Surf. Sci. 2008, 602, 960. [4] Schuttlefield, J. D.; Cox, D.; Grassian, V. H. J. Geophys. Res. 2007, 112 (D21303). [5] Delville, A. J. Phys. Chem. 1995, 99, 2033. [6] Tunega, D.; Gerzabek, M. H.; Lischka, H. J. Phys. Chem. B 2004, 108, 5930. [7] Tunega, D.; Benco, L.; Haberhauer, G.; Gerzabek, M. H.; Lischka, H. J. Phys. Chem. B 2002, 106, 11515. [8] Tunega, D.; Haberhauer, G.; Gerzabek, M. H.; Lischka, H. Langmuir 2002, 18, 139. [9] Eastwood, M. L.; Cremel, S.; Gehrke, C.; Girard, E.; Bertram, A. K. J. Geophys. Res.-Atm. 2008, 113.  113  Part I Appendices  114  Appendix A Simulation Models Simulation of a molecular system containing water and a surface such as kaolinite requires two basics things: models for the water molecules and the kaolinite surface, which include the Hamiltonian specifying the interactions between these two components, and a recipe for evolving the system towards an equilibrium state consistent with the Hamiltonian.  A.1  Molecular Models  In this thesis the water models used are the three-site SPC/E [1] and five-site TIP5PE [2]. The kaolinite lattice is described with the Clay Force Field (CLAYFF) recently developed by Cygan et al. [3]. A representation of both water models is given in Fig. A.1. In addition, the nonbonded parameters used for all the atoms included in our system, namely, the clay surface and water molecules, are given in Table A.1. Note here, that both the water molecules and kaolinite lattice are represented as interacting sites and that the models are rigid. The reader can refer to Chapters 2 to 5 for a more complete description of the models. The total configurational energy for our models is given by  Uconf =  1 4πǫo  i=j  qi qj + rij  Do,ij i=j  Ro,ij rij  12  −2  Ro,ij rij  6  ,  (A.1)  where the first and second terms are the Coulombic and Lennard-Jones (LJ) potentials, qi is the charge on site i, ǫo is the permittivity of free space, Ro,ij and Do,ij are the distance and energy parameters for the different LJ site-site interactions. The summations are over all interaction sites in the system including water molecules and the kaolinite 115  A.1. Molecular Models  species SPC/E water TIP5P-E water hydroxyl hydrogen hydroxyl oxygen bridging oxygen tetrahedral silicon octahedral aluminum  symbol wh,wo wh,wo,lp ho oh ob st ao  Do (kcal/mol) 0.1554 0.1780 0.1554 0.1554 1.8405 x 10−6 1.3298 x 10−6  Ro charge (˚ A) (e) 3.166 qwo = −0.8476, qwh = 0.4238 3.097 qlp = −0.2410, qwh = 0.2410 0.425 3.5532 -0.950 3.5532 -1.050 3.7064 2.100 4.7943 1.575  Table A.1: Nonbonded parameters for kaolinite ([Al2 Si2 O5 (OH)4 ]) using the CLAYFF Force Field [3] and water models. The subscripts wo, wh, and lp represent water oxygen, water hydrogen, and lone pair, respectively. lattice. For the cross interactions, the Lorentz-Berthelot rules, Do,ij =  Do,ii Do,jj and  Ro,ij = (Ro,ii + Ro,jj )/2, are used. The LJ potential is composed of two terms, an attractive and repulsive terms. The attraction comes from London dispersion forces. They are due to the induced dipoleinduced dipole and higher induced multipole forces. For a spherically symmetric particle, the induced dipole forces give rise to a potential uatt (r) = C6 /r 6 (LJ potential), and the higher multipoles give terms of all higher even powers in r. The inclusion of higher order terms depends on the level of accuracy required. The repulsive part of the LJ potential accounts for Pauli repulsion at short ranges due to overlapping electron orbitals. The LJ potential uses a mathematical form that represents well the repulsive part as an inverse power law urep (r) ∝ 1/r 12. The Coulombic term in Eq. A.1 describes the electrostatic interactions occurring between particles that have charges centered on the nuclei. This term gives rise to long-range electrostatic attractive and repulsive interactions. Now that molecular models have been described for all parts of the system, we can write explicitly an expression for the Hamiltonian as the sum of a kinetic and configurational energies (H = Ukin + Uconf ). Since the surface is kept rigid the clay atoms will only contribute to the configurational term. The total kinetic term is thus derived by the regular expression for the kinetic energy of a rigid body [4, 5] representing either SPC/E or TIP5P-E water molecules 116  A.1. Molecular Models  NW  |pi|2 + 2mW  Ukin = i=1 {i}  {i}  {i}  2  {i}  pω − pφ cos(θi ) cos(φi) − pθ sin(θi )sin(φi ) + 2I1 sin2 (θi )  {i}  2  {i}  pω − pφ cos(θi ) sin(φi ) − pθ sin(θi )cos(φi) 2I2 sin2 (θi )  (A.2)  {i} 2  pφ + 2I3  ,  where ω, θ, and φ are the Euler angles associated with the orientation of SPC/E or TIP5PE water molecules, I1 , I2 , I3 are the principal moments of inertia of a water molecule, {i}  p{ω,θ,φ} and pi are the conjugate momenta associated with the Euler angles and the center of mass momentum of water molecule i, respectively.  −0.241e o  1.0 A  o  0.9572 A o  o  109.47  3.166 A  o  104.52  o  109.47 0.70 A o  −0.8476e o  3.097 A +0.241e  +0.4238e  (a)  (b)  Figure A.1: Schematics of the SPC/E (left) and TIP5P-E (right) water models.  117  A.1. Molecular Models  A.1.1  Long-Range Forces  One of the limitations of computer simulations is the size of the system that can be simulated. The model systems must be small enough to be tractable. Large macroscopic systems are currently beyond the scope of computational power. When dealing with small systems of finite size, care must be taken to avoid any spurious surface effects. This limitation can be mitigated by the use of a technique called periodic boundary conditions. In this technique, a system of finite size is made infinite by periodically replicating the simulation box throughout space. Every molecule in any box has an image in all the other boxes. In practice, if the center of mass of a molecule leaves the box its periodic image will enter the same box through the opposite face. This method has the effect of illiminating boundaries in the central box, and hence surface effects. In this thesis, test were made to make sure that the size of the central box did not affect the results. Another limitation comes from the calculation of the potential energy. For a single Monte Carlo step, this calculation is of the order of the number of molecules in the system. Clearly, it is impossible to include the interactions of the infinite array of periodic images. For short-ranged interactions such as the van der Walls interactions (or Lennard-Jones potential) the solution is to invoke the minimum image convention, with a spherical cutoff at half the box length. In this scheme, when calculating the intermolecular interaction of any molecule, we place it at the centre of a box with identical dimensions to the central box. This molecule now interacts with the nearest molecules (or periodic images) falling in this region. Usually, it is sufficient to include the interactions of the molecules solely located inside the central box. However, for long-ranged interactions such as Coulombic interactions the minimum image convention method does not work as contributions from periodic images outside the cell will be significant. In the case of point charges, contributions from the periodic images is dealt with using Ewald sums. This method was initially developed for 3-D periodic systems (EW3D) [6, 7] but many methods were later developed and improved  118  A.1. Molecular Models for quasi two-dimensionally periodic systems (EW2D), such as flat surfaces or trenches, with a finite third dimension [8–11]. Another useful approach which is commonly used for quasi 2-D systems (and which was used in this thesis) is to apply the conventional Ewald 3-D method, given that a sufficiently large empty space is included in the z-direction to avoid any interaction with the periodic images [12] (see Fig. 2.1 for an example). When including the contribution of all ions in the central box and all periodic images the Coulombic term found in Eq. A.1 becomes Ucoul  1 1 = 2 4πǫo  N  N  i  j=1  ′ n  qi qj |rij + n|  ,  (A.3)  where n = (nx Lx , ny Ly , nz Lz ) is a vector identifying all the periodic images of the central box. The prime indicates that for n = 0 (i.e., the central box) the summation does not include the interaction of an ion with itself (i.e., i = j). For all other values of n the i = j interactions are included as they correspond to the interaction of an ion with its periodic images. The factor of 1/2 accounts for the double counting of interactions between the same pair of ions ij and ji. The Ewald 3-D method (EW3D) is implemented in two steps. Firstly, it effectively neutralizes all ions by surrounding them with a charge distribution, which is a Gaussian, of opposite charge centered on the ion ρqi (r) = qi κ3 exp(−κ2 r 2 )/π 3/2 ,  (A.4)  where the width of the distribution is determined by κ, and r is the position relative to the center of the distribution. This charge distribution has the effect of screening the interactions of the point charges. The combination of point ions and Gaussian charges improves greatly the summation as it is now short ranged in real space, where the minimum image convention is used. Secondly, a cancelling distribution equal and opposite to Eq. A.4 must be added in order to remove the influence of such a charge distribution on  119  A.1. Molecular Models the potential energy of the original system of point charges. The interactions due to the cancelling distributions is solved in reciprocal space as a Fourier series, then transformed back to real space. In addition to these changes a self-energy correction term arising from the interaction of the cancelling distribution with its own point charge must be removed to complete the Ewald sum. For molecular systems, another self-energy correction term is subtracted from the potential energy to account for intra-molecular Coulombic interactions. When using conducting boundaries the final result for the total electrostatic energy is given by [6]  Ucoul  1 1 = 2 4πǫo  N  N  ∞ ′  i=1 j=1  |n|=0  + (1/πL3 ) k=0  zi zj  erf c(κ|rij + n|) |rij + n|  zi zj (4π 2 /k 2 )exp(−k 2 /4κ2 )cos(k · rij ) (A.5)  N  κ  −  4π 3/2 ǫo  −  1 4πǫo  zi2 i=1 ns 2 κzia /π 1/2 +  i  a=1  1 2  ns  ns  zia zib erf (κdab )/dab  .  a=1 b=a  The complementary error function found in the first term falls to zero with increasing (κ|rij + n|). Thus, by choosing a large enough κ the sum in real space comes down to a sum over the particles in the primary cell with n = 0 (i.e., minimum image convention with a cutoff at half the box length). In calculating the reciprocal sum of the second term with vectors k = 2πn/L2 , it is important to keep a sufficiently large number of k-vectors in reciprocal space for the sum to converge to the appropriate potential energy for a given system. The third and fourth terms account for the self interactions between periodic images and intra-molecular interactions, respectively, and must be removed. In the last term ns , correspond to the number of sites on molecule i and dab is the intramolecular separation between sites a and b.  120  Bibliography [1] Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. [2] Lisal, M.; Nezbeda I.; Smith, W. R. J. Phys. Chem. B 2004, 108, 7412. [3] Cygan, R. T.; Liang, J.-J.; Kalinichev, A. G. J. Phys. Chem. B 2004, 108, 1255. [4] Goldstein, H. Classical Mechanics (Addison-Wesley, Reading, Massachusetts, 1980), 2nd ed. [5] McQuarrie, D. A. Statistical Mechanics (Harper Collins, New York, New York, 1976). [6] Allen, M. P.; Tildesley, D. J. Computer simulation of liquids, (Oxford University Press, Oxford, United Kingdom, 1987). [7] Ewald, P. Ann. Phys. 1921, 64, 253. [8] Parry, D. E. Surf. Sci. 1975, 49, 433. [9] Heyes, D. M.; Barber, M.; Clarke, J. H. R. J. Chem. Soc., Faraday Trans. II 1977, 73, 1485. [10] de Leeuw, S. W.; Perram, J. W. Mol. Phys. 1979, 37, 1313. [11] Spohr, E. J. Chem. Phys. 1997, 107, 6342. [12] Yeh, I.-C.; Berkowitz, M. L. J. Chem. Phys. 1999, 111, 3155.  121  Appendix B Simulation Methodology B.1  Monte Carlo Simulation Method  Now that we have defined the interactions among the particles included in our systems, we need a method for calculating the equilibrium macroscopic properties which allows easy treatment of a fluctuating number of water molecules. In statistical mechanics, the ensemble that represents this situation is called the grand canonical ensemble which has a fixed temperature, volume, and chemical potential, µ, (i.e., where the number of water molecules fluctuates). In agreement with this ensemble, the grand canonical Monte Carlo method was developed [1]. In this method, the configuration space is explored by attempting trial moves (water translation, rotation, insertion, and deletion). Each move is either accepted or rejected based on criteria that depend on whether the energy of the system rises or falls.  B.1.1  Statistical Mechanics  First, let us describe the basics principles of the Monte Carlo method using the canonical ensemble for a set of N identical atoms contained in a given volume (V ) at a constant temperature (T ). As for any other ensemble, calculation of the classical canonical partition function, ZN V T , is a central task  ZN V T =  1 h3N N!  dpN drN exp[−H(rN , pN )/kB T ] ,  (B.1)  122  B.1. Monte Carlo Simulation Method where rN and pN stand for the coordinates and momenta of all N particles included in the system, and H(rN , pN ) is the Hamiltonian of the system, composed of the kinetic (K) and potential (U) energy of the system (H = K + U). The partition function is a very important quantity that can be used to calculate average thermodynamic properties  A =  dpN drN A(rN , pN )exp[−βH(rN , pN )] , dpN drN exp[−βH(rN , pN )]  (B.2)  where β = 1/kB T . In Eq. B.2, the observable A is expressed in terms of the coordinates and momenta of the N-body system. The kinetic part (K) is a simple quadratic function of the momenta and can be easily separated and integrated out. However, when computing averages of functions A(rN ) the multidimensional integral that results is practically impossible to solve analytically and a numerical technique must be used. In principle, it would be possible to estimate the configurational integral by generating a finite but large number of particle configurations (τmax ), calculating the energy of each configuration, and weighing their contribution according to the Boltzmann factor (exp[−βU(r)]) as follows  ZN V T =  VN τmax  τmax  exp[−βU(rτ )] .  (B.3)  τ =1  Still, in order to evaluate such a sum a large number of random trial configurations is necessary. However, very few trial configurations randomly sampled from the Boltzmann distribution make a significant contribution. Again, this calculation is impractical due to the extremely large number of configurations necessary to obtain the correct answer. A much more efficient way of sampling the states is needed, which would concentrate on the states that make a more significant contribution to the probability distribution. Most of the time, we are not interested in calculating the partition function but rather  123  B.1. Monte Carlo Simulation Method in solving averages of the type  A =  drN A(rN )exp[−βU(rN )] . drN exp[−βU(rN )]  (B.4)  This is where the Monte Carlo method comes into play. Metropolis et al. [2] developed an algorithm that allows an efficient sampling of ratios such as in Eq. B.4. It provides a very efficient way to explore the phase space. The Metropolis Monte Carlo algorithm creates a Markov chain which is a sequence of random configurations or states which satisfies the following necessary conditions: a) each trial belongs to a finite set of outcomes; and b) each trial only depends on the outcome of the trial that immediately precedes it. A transition probability πmn links any two states m and n of the Markov chain. This transition is the probability of going from state m to state n. The matrix elements of πmn must obey one necessary condition: when an equilibrium distribution is reached the average number of accepted moves from state m to any other state n is exactly canceled by the reverse transitions. This condition is known as detailed balance or microscopic reversibility  ρm πmn = ρn πnm ,  (B.5)  where ρm and ρn are the elements of the limiting probability distribution (ρ) corresponding to states m and n. The transition probability matrix must satisfy other conditions. It must conserve probability πnm = 1 .  (B.6)  n  Additionnally, πnm must have no memory of the previously attempted moves. In other words, the transition probability matrix must only depend on states m and n. Finally, we require that a certain number of different ways exist to go from state m to state n. When all these conditions are satisfied, a Markov chain leading to the correct statistical 124  B.1. Monte Carlo Simulation Method distribution functions is produced and is guaranteed to give proper equilibrium properties. Metropolis et al. [2] showed that a transition probability matrix linking two consecutives states m and n has the following form πmn = Cmn πmn = Cmn  ρn ≥ ρm  ρn ρm  (B.7)  ρn < ρm .  It is composed of the product of two quantities, πmn = αmn Amn , with αmn being the probability that a transition between states m and n is attempted and Amn is the probability that the transition is accepted.  B.1.2  Grand Canonical Monte Carlo Simulation  We now make use of the Metropolis method [1, 2] to determine the transition probabilities for the different possible transitions occurring in our system based on the grand canonical ensemble where the number of water molecules in the system, NW , fluctuates according to a fixed chemical potential. The grand canonical partition function for this system is expressed by ∞  ZµV T = NW  (NW !)−1 eβµW NW (2Λ3W Θ1 Θ2 Θ3 )NW =0  2π  π  2π  −L  2π  R  e−β(Uconf (NW ,{ri },{Ωi }) x  ... 0  0  0  L  0  0  NW  sin(θi )dΩi d3 ri , i=1  (B.8) where the momentum terms have been integrated out. In Eq. B.8, ΛW and Θ1 , Θ2 , and Θ3 are the translational thermal wavelength and rotational thermal angles of the water molecules, respectively ΛW = √  h , 2πmW kB T  125  B.1. Monte Carlo Simulation Method Θi = √  h , 2πIi kB T  where h, mW , Ii , kB , and T , are the Planck’s constant, the molecular mass of water, i th principal moment of inertia of a water molecule, Boltzmann’s constant, and temperature, respectively. At each step either a particle insertion, deletion, or displacement (rotation or translation) is attempted with equal probability δ. Thus, these three types of moves are selected randomly each with a probability of δ = 1/3. In addition, when a displacement is selected there is a probability of 1/2 that either a rotation or a translation is attempted. When attempting a particle insertion the water molecule is given a random orientation and a random location restricted to the interlayer between the two slabs. Every attempt described above is either accepted or rejected with a certain probability. If the move is downhill in energy it is immediately accepted. If the move is uphill in energy a random number is generated between 0 and 1, and the move is accepted upon the condition that the random number is smaller than the corresponding Boltzmann factor for the change in the appropriate property (see below) between the two states. We now define the different acceptance probabilities for the different transformations occurring in our grand canonical Monte Carlo simulation [3]. 1. Displacement of a water molecule: A water molecule is selected at random and is displaced a randomly selected distance along one of the three cartesian axis. The attempt and acceptance probabilities are given by 1 − 2δ 2NW  dxdydz V  ,  (B.9)  Amn = min 1, e−β(Un −Um )  .  (B.10)  αmn =  2. Rotation of a water molecule: A water molecule is selected at random and is rotated randomly about all three Euler axis. The attempt and acceptance probabilities are given  126  B.1. Monte Carlo Simulation Method by αmn =  1 − 2δ 2NW  sin(θ)dωdθdφ 8π 2  Amn = min 1, e−β(Un −Um )  ,  (B.11)  .  (B.12)  3. Water molecule insertion: A water molecule is inserted by placing the centre of mass at a random location, and then rotating the water molecule to a random orientation. The attempt probability is αmn =  δsin(θ)dxdydzdωdθdφ . 8π2 V  (B.13)  The acceptance probability is  Amn = min 1,  4π 2 V e−β[(Un −Um )−µW ] Λ3W Θ1 Θ2 Θ3 (NW (n) + 1)  .  (B.14)  4. Water molecule deletion: A water molecule is chosen at random and is removed from the system. The attempt and acceptance probabilities are given by  αmn =  Amn = min 1,  δ , NW  Λ3W Θ1 Θ2 Θ3 NW (n) −β[(Un −Um )+µW ] e 4π 2 V  (B.15)  .  (B.16)  In the above equations, Um − Un is the difference in potential energy between the two states, V is the total volume water molecules can occupy in our systems, NW is the number of water molecules, and µW is the chemical potential of water. Repetition of the above process satisfies microscopic reversibility and produces a Markov chain leading to the correct statistical distribution functions and proper equilibrium properties. To ensure rapid convergence the acceptance rates of a rotation or translation were set to approximatively 50% by adjusting the maximum displacement parameters. For the atomistically flat kaolinite surfaces a minimum of 5x107 Monte Carlo (MC) steps were necessary to attain equilibrium. However, in the case of trenches very 127  B.1. Monte Carlo Simulation Method long simulations were necessary before equilibrium was reached (∼ 1x109 MC steps or more). Following equilibration, averages were collected for a minimum of 2x108 MC steps and included, energy, water content, density and order parameter profiles, and hydrogen bond number. This procedure was repeated for a wide range of chemical potentials and at 235K and 298K for both SPC/E and TIP5P-E water to produce adsorption isotherms for all kaolinite surfaces and trenches studied.  128  Bibliography [1] Allen, M. P.; Tildesley, D. J. Computer simulation of liquids, (Oxford University Press, Oxford, United Kingdom, 1987). [2] Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087. [3] Lakatos, G. W. Transport and Structure in Nanoscale Channels (University of British Columbia, Vancouver, Canada, 2006).  129  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0060989/manifest

Comment

Related Items