Can Ice-Like Structures Form on Non-Ice-Like Substrates? The Example of the K-feldspar Microcline.

Feldspar minerals are the most common rock formers in Earth's crust. As such they play an important role in subjects ranging from geology to climate science. An atomistic understanding of the feldspar structure and its interaction with water is therefore desirable, not least because feldspar has been shown to dominate ice nucleation by mineral dusts in Earth's atmosphere. The complexity of the ice/feldspar interface arising from the numerous chemical motifs expressed on the surface makes it a challenging system. Here we report a comprehensive study of this challenging system with ab initio density functional theory calculations. We show that the distribution of Al atoms, which is crucial for the dissolution kinetics of tectosilicate minerals, differs significantly between the bulk environment and on the surface. Furthermore, we demonstrate that water does not form ice-like overlayers in the contact layer on the most easily cleaved (001) surface of K-feldspar. We do, however, identify contact layer structures of water that induce ice-like ordering in the second overlayer. This suggests that even substrates without an apparent match with the ice structure may still act as excellent ice nucleating agents.

Introduction

Feldspar minerals belong to the group of tectosilicates which make up nearly 75% of the crust of the Earth. Their abundance makes them relevant to various fields such as water contamination, soil, geochemistry, and atmospheric sciences.[1,2] Despite their ubiquity, gaps in our understanding of the physical chemistry of these materials remain. How water adsorbs on their surfaces, how ice might form, and aspects of feldspar dissolution behavior are prominent examples.

One particular mineral, microcline (a K-feldspar), was put into the spotlight recently. Using a droplet freezing technique, Atkinson et al. showed that K-feldspar dominates ice nucleation in Earth’s atmosphere by mineral dust.[3] While these experiments are well suited to establishing ice nucleating ability, they did not provide the atomistic insight necessary to understand why K-feldspar is such an excellent ice nucleator. This work has prompted a number of follow up studies,[4−10] and of course water mineral surfaces, although not specifically about microcline, have previously been widely studied.[2,11−16] The findings of Atkinson et al.[3] are also to some extent surprising because feldspar surfaces bear no obvious resemblance to the structure of ice. It is commonly thought that good ice nucleators should present surfaces that template ice-like structures. Lattice match as in the case of AgI[17] or hexagonal arrangements of hydroxyl groups in the case of kaolinite[18] are famous examples. Some experiments suggest, however, that other materials without an apparent templating effect can still be good ice nucleators. One example is testosterone,[19] and others are marine biogenic particles.[20] These cases not only demonstrate that lattice match is probably not a good predictor for ice nucleating ability, but also that we still lack an understanding of what it is that makes a good ice nucleator.[21]

Recent experimental studies have focused on establishing the ice nucleating ability of various materials,[4−6,8,9] and recent simulations have aimed at characterizing specific water–solid interfaces[22−28] and at understanding general trends in ice nucleating on more idealized model systems.[29−36] They were able to show that it is a subtle balance between lattice match, surface topology, and adsorption energy that impacts ice nucleation. This means that ice nucleating ability is strongly system dependent and requires investigations on a case by case basis. Therefore, by looking at a variety of examples of good ice nucleators, one might be able to deduce a general theoretical framework for what makes materials good at nucleating ice. Understanding what it is that makes particularly interesting substrates such as feldspar so good at nucleating ice is therefore desirable. Nucleation on real feldspar particles, however, is a complicated process, not least because many different surfaces will be exposed on a real feldspar particle. Furthermore, defects, trenches, and wedges might play a role in the nucleation process. No single study can address all of these issues, but here we lay the groundwork by establishing the structure of the clean and water covered microcline (001) surface, which is the surface that is cleaved most easily for feldspars.[13]

The rest of the paper is structured as follows: In section we elaborate on the computational method. We then discuss the bulk structure of feldspar in section for two reasons. First, we show that our computational setup for our DFT calculations is able to reproduce experimental findings very well, which gives us confidence in our methods. Second, understanding the bulk is a prerequisite to the surface and the surface/water interface. After covering bulk structures we focus on the surface, for which much less experimental data are available, in section . We show that the thermodynamics of the Al/Si distribution on the surface differ significantly from that in bulk and connect this with the dissolution behavior of feldspar minerals. In section we move on to the mineral water interface, where we discuss the interaction of water with the microcline (001) surface. At first, we focus on monolayer adsorption and then follow this by a discussion of multilayer ice growth. We demonstrate that the “ice-unlike” structure of feldspar does not lead to ice-like structures in the contact layer. However, the particular structure of the contact layer has the ability to induce ice-like structures in the second water layer. It is these ice structures that could make feldspar the excellent ice nucleator it is. As well as providing initial understanding of clean and water covered microcline we hope this study will stimulate the general discussion of what makes a material a good ice nucleator.

Computational Method

Density functional theory (DFT) calculations were performed with the periodic plane wave code VASP.[37,38] Core electrons were replaced by projector augmented wave (PAW) potentials,[39] and valence electrons were expanded in plane waves. All calculations were done on microcline, a triclinic potassium feldspar with chemical composition K4Al4Si12O32 (space group C1̅). Its structure is shown in Figure a and b.

Figure 1

(a) One possible arrangement of AlO4– (light blue) and SiO4 (dark blue) tetrahedra in microcline. K+ ions are shown in pink and oxygen atoms in red. A primitive unit cell contains 4 AlO4– tetrahedra, 4 K+ ions (to balance the charge), and 12 SiO4 tetrahedra. T1 sites expose bridging oxygen atoms to the (001) surface (these oxygens are highlighted with green circles), and T2 expose bridging oxygens to the (010) surface (oxygens highlighted with yellow circles). (b) Full atomistic model of the unit cell. Al atoms are depicted in blue, Si atoms in yellow, K+ ions in pink, and oxygen atoms in red. This color code is used for feldspar throughout the paper. All eight T1 and eight T2 sites are labeled. On the bottom of (b) we show the orientation of the unit cell, as well as the location of the two most stable cleavage planes in (a). (c) Typical slab model used for surface and interface calculations. Water fragments (OH and H) used to saturate undercoordinated bonds are shown in green.

Starting from the experimental unit cell and structure[40] we did a full geometry relaxation with a 600 eV cutoff for the plane wave basis set, which is 50% larger than the recommended 400 eV cutoff for ionic relaxation so as to avoid problems with Pulay stresses (see the Supporting Information (SI) for more information). This optimization was followed by an ionic relaxation with a 400 eV cutoff to obtain the equilibrium geometry for any bulk structure. All other calculations were also carried out with a 400 eV cutoff. Electronic structures were converged to within 1.0 × 10–5 eV, and ions were relaxed until the forces acting on them were below 0.005 eV/Å. K-points were sampled with a 6 × 4 × 6 Monkhorst–Pack k-point mesh[41] for bulk calculations and a 3 × 2 × 1 Monkhorst–Pack k-point mesh for all slab calculations. The (001) slab was generated by cleaving a 1 × 1 × 2 supercell and saturating all undercoordinated atoms (Si, Al, or O) by dissociating water molecules to form Al–OH and Si–OH groups where necessary. A vacuum parallel to the surface normal of at least 12 Å was used to separate slabs in adjacent cells, and a dipole correction along the surface normal was applied.[42,43]Figure c shows the slab model resulting from this procedure.

We used the Perdew–Burke–Ernzerhof (PBE)[44,45] exchange-correlation functional for all calculations, but we also used the optB88-vdW[46−48] functional to test the role of van der Waals forces. Both functionals lead to the same qualitative conclusions, as shown in the SI. The stability of bulk structures with different Al/Si order (for a precise definition see sec. ) was investigated by comparing their total energy Ebulk relative to the most stable bulk structure with energy Ebulkmin per Al atom. In the primitive unit cell, four Al atoms can be distributed among 16 sites, which gives rise to a total of 1820 possible structures. Using the experimental crystal structure of microcline,[40] all 1820 structures were generated and minimized with the classical force field ClayFF[49] using GROMACS.[50,51] The four most stable bulk structures found for every distinct Al order were then optimized with DFT as described above. Additionally, several structures were chosen by hand (some of them at random, some based on chemical intuition) and optimized to guarantee a good sampling of configurational space. In the SI we show that energies obtained with PBE and ClayFF do indeed correlate well with each other, which demonstrates that this a reasonable procedure. In total 38 different bulk configurations were considered with DFT.

The stability of different Al/Si orders on the surface was evaluated in a similar manner. For this we compared the total energy of surfaces Esurf with different Al order relative to the most stable surface structure with energy Esurfmin per Al.

Water adsorption energies were calculated as followswhere EH2O/surf denotes the total energy of the slab with water adsorbed on it, n the number of water molecules adsorbed, EH the total energy of a water molecule in vacuo, and Esurf the total energy of the slab in vacuo. With this definition, a negative adsorption energy corresponds to favorable (exothermic) adsorption. The dimensionality of the configurational space is too large to perform a comprehensive search of all adsorption minima, especially for multilayer water structures which contain up to 24 water molecules. Therefore, a guided structure search approach based on the potential energy surface (PES) of water on microcline, chemical intuition based on maximizing the number of hydrogen bonds formed, and structural motifs found in ice Ih was chosen to address this problem. The key idea is to guide the structure search in such a way that mainly structures which either have a strong interaction with the surface or strong intermolecular bonds are searched for, leaving out all the structures that are of unlikely relevance. We prescreened the configuration space with ClayFF[49] combined with the SPC water model,[52] and low-energy candidates were then further relaxed using DFT. Detailed information about our procedure is provided in the SI. In total, more than 100 000 structures were screened with ClayFF, and more than 150 structures were minimized with DFT.

Results and Discussion

Bulk

Before investigating the water/feldspar interface one needs to have an atomistic understanding of the bulk. Microcline’s bulk structure is fairly well established, but there are issues over the Al order at an atomistic level that remain unclear. We start with an introduction of the general structural features of microcline, which is followed by a discussion of Al order from two perspectives. First, the effect of Al order on the unit cell dimensions is investigated, and second, the energetics of different Al order states are discussed. Our results show that DFT is suitable to describe the structure and energetics of feldspar well by comparing our results to experiment. They also reveal, however, that a microscopic structural feature, the appearance of two distinct Al–O bonds, is not accounted for in macroscopic models.

Structure of Microcline

The structure of feldspars in general can be described as “mirrored crankshaft-chain” frameworks of polymerized Al/SiO4 tetrahedra.[53] Both Al and Si atoms are 4-fold coordinated by oxygen, forming AlO4 and SiO4 tetrahedra. K+ ions counterbalance the formal −1 charge of AlO4 tetrahedra. Figure a shows the unit cell of microcline with AlO4 and SiO4 tetrahedra highlighted. The primitive unit cell contains a total of 16 tetrahedra, of which 4 have Al and 12 have Si in the center, as well as 4 K+ ions. Aluminum can potentially be located in the center of any of the 16 tetrahedra making up the unit cell. The possible arrangements of Al among these sites is referred to as Al/Si disorder and has been the subject of numerous studies; see, e.g., ref (54).

One can distinguish two different tetrahedral sites in monoclinic feldspars, namely T1 and T2 sites. T1 sites expose oxygen atoms to the (001) surface (Figure a, green circles), while T2 sites expose oxygen atoms to the (010) surface (Figure a, yellow circles). To clarify this nomenclature, the location of T1 and T2 sites is also highlighted in Figure b. In a unit cell there are eight T1 and eight T2 sites.[56] If all Al occupy the same type of site, the structure is said to be fully ordered, if they are equally distributed between T1 and T2 sites, the structure is said to be fully disordered. We stress that with order we only refer to the distribution of Al between T1 and T2 sites throughout the paper. Besides being relevant to the physicochemical properties of feldspar minerals, Al/Si order also plays a crucial role in, for example, zeolites.[57]

Al Order

We address the Al order from two different viewpoints in this section. First, we discuss the relation between Al order and unit cell dimension. Then we discuss the energetics of different Al order states.

A connection between Al order and feldspar unit cell dimension has long been known.[58−60] For a detailed explanation of the origin of this connection, the interested reader is referred to ref (61). For present purposes it is sufficient to state that Kroll and Ribbe[55] proposed the following widely used connection between Al order/disorder and unit cell geometryIn this equation, T1 takes a value between 0 (all Al are in T2 sites) and 1 (all Al are in T1 sites). b and c* are unit cell and reciprocal unit cell vectors, respectively (measured in Å and Å−1). This or similar methods[55,62] have been useful in understanding Al disorder in numerous minerals. However, they all rely on measurements performed on macroscopic samples in which the data are averaged over multiple polycrystalline domains. Computational studies have addressed Al order issues using classical force fields and Monte Carlo methods (see for example refs (63 and 64) and references therein). Using DFT this relationship can be tested at the atomistic level, which, to the best of our knowledge, has not been done before.

For every bulk structure computed we estimated the predicted Al disorder using eq from the unit cell parameters obtained from the lattice relaxation. Figure shows the results. The x-axis represents the actual Al order in a structure and the y-axis the order estimated with eq . The solid red line shows the predictions based on the model,[55] and the black data points are estimates based on our DFT results. Overall we find that the model and our DFT estimates agree very well with each other. The linear relationship between cell dimension and T1 order is clearly well reproduced. What is interesting to note, however, is that different structures at any given T1 occupation result in a number of different estimated disorder states, as can be seen from Figure . In other words, the suggested 1:1 correspondence between estimated T1 order and actual order breaks down at the atomistic level. This might not be too surprising considering that the model itself is based on data averaged over large samples. The underlying reason is that eq assumes one characteristic Si–O bond length and one characteristic Al–O bond. However, two different Si–O bond lengths can be identified, depending on whether Al atoms are next nearest neighbors to Si or not. We provide a more detailed explanation of why this breakdown happens at the microscopic level in the SI.

Figure 2

Connection between Al order and bulk unit cell geometry. The x-axis shows the Al occupation of T1 sites. T1 = 0.00 means that all Al are in T2 sites, T1 = 1.00 that all Al are in T1 sites. We estimate the Al disorder based on DFT lattice relaxations combined with the model proposed by Kroll and Ribbe.[55] The y-axis shows the T1 occupation estimated using this model (eq ). The red line shows the behavior eq predicts. Each black circle represents a different structure which was fully optimized with DFT. Because we optimized several structures with the same actual T1 order that differed in their atomistic structures, multiple data points for each Al/Si order were obtained. The model and DFT agree very well on the trend. Also shown on the top are the most stable bulk structures found for T1 = 0.00 and T1 = 1.00.

After having considered the effect of Al order on the unit cell geometry, we now discuss its impact on the energy. Figure plots Al occupation of T1 sites as a function of energy relative to the most stable bulk structure. On the x-axis the Al distribution in T1 sites is plotted, 0.0 meaning that all Al are located in T2 sites and 1.0 meaning that all Al are located in T1 sites. Our calculations show that in microcline Al prefers to occupy T1 sites exclusively. This again is in very good agreement with experiment,[55] which shows that microcline is highly ordered, with Al occupying almost exclusively T1 sites. The second most stable order state is for all Al to occupy T2 sites, which is 30 meV per Al atom higher in energy. This is considerably more than the thermal energy available at room temperature, which agrees well with the fact that microcline is found experimentally to be highly ordered, with Al concentrated in T1 sites. Besides energetics many other factors, such as conditions under which the mineral formed (temperature, pressure, pH), will contribute to the final composition of a feldspar mineral. None of these are taken into account here; however, our results show that it is thermodynamically favorable to highly order Al into T1 sites.

Figure 3

Relation between Al order and energetics of microcline (bulk and surface). The x-axis shows the Al occupation of T1 sites. T1 = 0.00 means that all Al are in T2 sites and T1 = 1.00 that all Al are in T1 sites. The energy per Al atom of different bulk (left y-axis, black dots) and surface structures (right y-axis, red squares) relative to the most stable structure is plotted. For clarity, only the most stable structure identified for each T1 occupation is shown. We find that the most stable structures are the ones with T1 = 1.00 for bulk microcline and T1 = 0.00 for the surface. The topmost parts of a T1 = 0.00 and a T1 = 1.00 surface are also shown. While the surface corresponding to T1 = 0.00 only exposes Si–OH groups, the surface at T1 = 1.00 has 50% Si–OH and 50% Al–OH groups.

Each of the most stable bulk structures at a given coverage obeys Löwenstein’s rule (Al–O–Al linkages are unfavorable);[65] no such Al–O–Al linkages appear in any of the most stable structures. Furthermore, Dempsey’s rule, which states that the number of Al–O–Si–O–Al linkages shall be minimized[66] at a given Al/Si ratio, also holds, with the number of such linkages being four for each of the most stable structures.

To summarize our bulk results: we established that our method agrees with experimental findings. We can reproduce both the empirical relation between Al order and lattice dimensions as well as the preferential ordering of Al into T1 sites. Furthermore, we find that there are two distinct Si–O bond lengths instead of only one, depending on local environment around that Si–O bond. Having an atomistic understanding of the bulk, one can now move on to the surface.

The (001) Surface

Here we focus on the (001) surface because this is the most easily cleaved surface[13] and therefore most likely to play a significant role in dissolution and ice formation. We will show that the energetics of Al order on the surface differ significantly from the bulk. While Al in bulk prefer T1 sites, they appear to be more stable in T2 sites on the surface.

To investigate the surface structure of microcline (001), slab models from the fully relaxed bulk structures with different Al/Si order were built, as discussed in section . The (001) cleavage plane is shown in Figure a. Cleaving the bulk structure along this plane leads to the breaking of Si–O and Al–O bonds. The resulting undercoordinated atoms were saturated with dissociated water molecules, and the resulting slab model is shown in Figure c.

Upon relaxation, none of the surfaces underwent significant reconstructions. This is consistent with previous studies on for example mordenite (001).[12] Our main goal here is to understand how surface stability depends on Al order. We show the results in Figure , where the x-axis corresponds to Al T1 site occupation and the y-axis on the right of the graph is the stability of a given surface relative to the most stable surface, Esurfmin. Note that the Al order is directly linked to the number of Al–OH groups on the surface. Structures with T1 = 0.00 have only Si–OH groups exposed on the surface, whereas structures with T1 = 1.00 have 50% Al–OH and 50% Si–OH groups on their surface.

Our main finding is that, in contrast to the bulk, surfaces with all Al in T2 sites (T1 = 0.00, no Al–OH groups) are most stable. The fact that surface composition and bulk composition differ suggests that the chemical difference between them, the presence or absence of hydroxyl groups, is the driving force for this. Our result simply implies that SiOH groups are more stable than AlOH groups. The likely reason behind this is the stronger electrostatic attraction between a Si4+ and an OH– compared to an Al3+ and an OH–. The result agrees well with experimental evidence showing that Al groups dissolve prior to Si groups. Previous computational studies[1,2,14,16,67,68] consistently found that, based on kinetic arguments (reaction barriers), Al should dissolve faster than Si. Our results add to this in that they show that it is also thermodynamically more favorable to enrich the surface in Si–OH groups.

Water/Feldspar Interface

Having established the structure of the (001) surface of microcline, we now look at its interaction with water. In this section we will look at the water/feldspar interface ranging from water monomers to multilayer ice structures on top of microcline. Even though monolayers identified do not resemble the structure of ice, we find second overlayers that are indeed ice-like. This shows how even substrates with a non-ice-like surface morphology can stabilize ice-like multilayers.

The feldspar surface is rather complex, and it exhibits a variety of chemical motifs with which water can interact. Specifically, there are hydroxyl groups attached to either Si or Al, which can form hydrogen bonds with water, either as hydrogen bond donors or acceptors. In addition there are K+ ions which can form coordinative bonds with H2O. Lastly and of less importance, bridging oxygen atoms (Al–O–Si or Si–O–Si) can act as hydrogen bond acceptors.

To understand how water interacts with microcline (001), the adsorption of water in a broad range of adsorption structures over a wide range of coverages was examined. Two particular surfaces were considered, the most stable surface (T1 = 0.0) and a surface with Al highly ordered in T1 sites, T1 = 1.0. The former will be exposed in weathered feldspar samples, whereas the latter will be present in freshly cleaved samples. We discuss the differences between Al- and Si-rich surfaces whenever necessary, but our main focus will be on the freshly cleaved surface. One reason for focusing on this (T1 = 1.0) surface rather than the most stable (T1 = 0.0) surface is that there are major chemical processes happening during the weathering process. It is for example well established that K+ ions will partially be replaced with H3O+ ions on the surface.[54] Furthermore, because weathering happens in non-neutral conditions, hydroxyl groups and bridging oxygens will be partially protonated/deprotonated depending on the pH. Having a realistic representation of a weathered microcline surface is therefore beyond the reach of this work. Looking at the surface that will be exposed after fresh cleavage on the other hand involves considerably fewer approximations concerning surface structure. In this study we therefore focus on a clean, nonweathered feldspar surface. All issues arising from the fact that real feldspar surfaces in nature will be imperfect are beyond the scope of this work.

The particular T1 = 1.0 surface chosen was one of the most stable T1 = 1.0 surfaces with Al distributed between different T1 sites so as to maximize the number of different chemical motifs on the surface. The slab considered had adjacent SiOH–SiOH, AlOH–AlOH, and AlOH–SiOH groups. Figure a shows results for water adsorption on this surface and a T1 = 0.0 surface as a function of coverage.

Figure 4

(a) Adsorption energy of water as a function of coverage. We define 1.0 monolayer to consist of 8 water molecules because there are 8 strong interaction sites for water available in a primitive unit cell (6 hydroxyl groups and 2 potassium ions). The dashed line represents the cohesive energy of bulk ice Ih from PBE (−0.66 eV/water). Blue data points correspond to adsorption energies on a T1 = 1.0 surface and red data points to energies on a T1 = 0.0 surface. Coverages up to 1.5 ML form a single contact layer with microcline. Energies of multilayers (coverage of 2.5 and 3.0 ML) are also shown. Overlayer structures 1.0 ML + basal and 1.0 ML + prism both have a coverage of 2.5 ML, and overlayer structures basal + basal and prism + prism both have a coverage of 3.0 ML; they are therefore difficult to distinguish on the basis of coverage. The inset shows the energetics of different multilayer structures more clearly, with (i), (ii), (iii), and (iv) being 1.0 ML + basal, 1.0 ML + prism, prism + prism, and basal + basal, respectively. (b) to (e) show the most stable adsorption structure found for a water monomer (0.125 ML) and monolayers with coverage 0.5, 1.0 and 1.5 ML, respectively. For clarity, we only show the unit cell in (b), and (c)–(e) have the same unit cell. We also add red lines connecting neighboring water molecules to highlight the network H2O forms.

Water Monomers and Clusters

The most stable monomer adsorption site is in between one AlOH and one SiOH group, as shown in Figure b. In this configuration, water interacts with two hydroxyl groups, receiving one hydrogen bond (from Si–OH) and donating another (to Al–OH). The remaining hydrogen atom of the water molecules forms an additional hydrogen bond with a bridging oxygen atom meaning that the water monomer is involved in three hydrogen bonds. The adsorption energy of a water molecule in this position is −0.75 eV. On the most stable T1 = 0.0 surface, the adsorption energy is significantly weaker, only −0.59 eV. Because Al–OH groups are more reactive than Si–OH groups (hence the difference in surface energies) they bind water more strongly. Even though the energetics of water is quite different on both surfaces, structural features of an adsorbed water molecule are quite similar. The most stable adsorption site for a water monomer on the T1 = 0.0 surface is the same as the one shown in Figure b. On the T1 = 1.0 surface the O–O hydrogen bond lengths involving the water molecule and the surface are 2.65 Å (hydrogen bond accepted from Si–OH), 2.71 Å (hydrogen bond donated to Al–OH), and 3.13 Å (hydrogen bond donated to bridging O) for the water monomer on the T1 = 1.0 surface. On the T1 = 0.0 surface, these hydrogen bond lengths are 2.71 Å (hydrogen bond accepted from Si–OH), 2.84 Å (hydrogen bond donated to Si–OH), and 3.26 Å (hydrogen bond donated to bridging O). The general elongation of hydrogen bonds on the T1 = 0.0 surface agrees with the weaker adsorption energy.

In this adsorption site, the water monomer remains intact. We checked for stable dissociated states but did not find any that came within 1.5 eV of the intact monomer. The fact that water monomers do not dissociate on this surface is not surprising since our surface model is fully hydroxylated.[69]

Upon increasing the coverage water dimer adsorption was examined, again considering a broad range of structures. Metastable hydrogen bonded dimers were identified. However, their adsorption energy was at least 0.05 eV/H2O less than the adsorption energy of two separate adsorbed monomers. Similar behavior was seen for higher coverages: isolated water monomers were more stable than hydrogen bonded clusters. This means that water clustering on a freshly cleaved microcline (001) surface is not favored. However, on the most stable surface (T1 = 0.0) this behavior at low coverages is qualitatively different. Because of the weaker adsorption energy, clustering of water is favorable on the T1 = 0.0 surface. Upon increasing coverage, the water adsorption energies become more similar on the two surfaces (see Figure a). This is not surprising because at higher coverages the relative importance of surface water interactions versus water–water interactions decreases. As with the monomer, the structural features of the adsorbed water cluster do not differ much, and changing the surface mainly impacts the adsorption energies.

Monolayer Coverage

To understand how ice can grow on surfaces, it is crucial to look at structures that might form on the substrate at high coverages. At these high coverages, the differences in terms of water adsorption between the two different surfaces considered become very small, and so our focus in the discussion will be on the freshly cleaved (001) surface. We looked at a wide range of ice structures including traditional ones (basal and prismatic faces) and nontraditional ones adapting to the potential energy surface of the substrate.

Figure d shows the most stable contact layer structure with an adsorption energy of −0.56 eV/H2O found at a coverage of 1.0 ML. This structure is composed of square-like arrangements of water molecules. Three out of four water molecules in a square form hydrogen bonds with surface hydroxyl groups, and the fourth water molecule is coordinated to a K+ ion. The average intermolecular water O–O bond length in this structure is 2.87 Å, and the O–O–O angles range from 77° to 103° (compared to 109° in ice I). This contact layer is rather flat, the average height difference between nearest neighbor water molecules being 0.3 Å (as opposed to 0.9 Å in a basal bilayer), and does not resemble any traditional ice-like structure.

Overlayers that would resemble the structure of an ice I basal or prism layer have a coverage of 1.5 ML. In terms of energetics, overlayers at 1.0 and 1.5 ML are very similar, the most stable structure having an adsorption energy of −0.60 eV/H2O. The most stable contact layer at such a coverage is shown in Figure e. Irrespective of our starting point (basal or prism face) for the structure search, low energy structures converge to structures similar to the one shown in Figure e. In these structures the hexagonal arrangement of water molecules persists (as can be seen from the top view). While ice bilayer structures have water molecules located at two well-defined heights, water molecules in these monolayers are at a broad distribution of heights, as exemplified by the side view in Figure e. In the most stable 1.5 ML ice structure on microcline, intermolecular O–O bond lengths range from 2.59 Å to 2.91 Å and intermolecular O–O–O angles from 95° to 128°. Low energy monolayer coverage structures at 1.5 ML therefore do not resemble either basal or prism structures. The broad range of bond lengths and angles is a result of the water overlayer optimizing its interaction with the surface, which comes at the cost of distorting hydrogen bonds.

Previous work has shown that, for example, on metal surfaces ice bilayers do not form in general.[70−73] Even surfaces with a near perfect lattice match to ice do not necessarily tend to form ice bilayers.[74] Interestingly on other surfaces, such as for example kaolinite, more regular ice-like structures can form in the contact layer.[18] In this case this is made possible by the particular arrangement of hydroxyl groups on the kaolinite surface. They are arranged in an almost uniform hexagonal pattern and match the ice bilayer structure rather well. A hexagonal layer of ice can therefore bind very strongly to the kaolinite surface.[18] It should be noted that the hexagonal structure on kaolinite has no dangling hydrogen bonds and is therefore not amenable to further water adsorption. On microcline, the situation is very different because the surface does not match the bilayer structure at all. Not only are different interaction sites exposed at different heights (the height difference between K+ ions and the oxygen atom of a OH group is 0.6 Å at the surface) but also the spacing between different hydroxyl groups itself varies significantly (O–O distances range from 4.7 Å to 6.1 Å). As a result of this very complex surface, we do not observe contact layers resembling a regular ice I bilayer structure.

Multilayer Coverage

The ability of a substrate to stabilize multiple ice layers appears to be a trivial prerequisite for any substrate to be a good ice nucleator. Interestingly, the contact layer on kaolinite was found to be hydrophobic.[18] Based solely on this, kaolinite should not drive multilayer ice growth, which would make it a bad ice nucleator. Classical MD simulations addressing ice nucleation on kaolinite revealed however that it nevertheless promotes ice growth.[75,76] This clearly shows that looking exclusively at the most stable structures is insufficient. At least of equal importance is to understand the full configurational space of water on a substrate and any low energy metastable structures that might form. In the case of kaolinite it was for example a prismatic face of ice forming on the surface that drove ice nucleation and not the 0.04 eV/H2O more stable bilayer-like structure obtained from first-principles calculations.[18]

To explore this system we performed an extensive search of possible ice structures using our guided structure search method. We specifically looked at four different ice multilayer structures on microcline: (i) 1.0 ML in the contact layer plus an ice I basal face in the second layer (1.0 ML + basal) (coverage of 2.5 ML); (ii) 1.0 ML in the contact layer plus an ice I prism face in the second layer (1.0 ML + prism) (coverage of 2.5 ML); (iii) two layers of the ice I prism face (prism+prism) (coverage of 3.0 ML); (iv) two layers of the ice I basal face (basal+basal) (coverage of 3.0 ML). The energies of the different multilayer structures are plotted in Figure a as a function of coverage and, in the inset, for the four structures just described. Figure shows the most stable structure for each of these cases from a side and top view.

Figure 5

Multilayer ice growth on microcline (001). We show the most stable structure found for a double basal layer (a), a double prism layer (b), a combination of a contact layer with 1.0 ML coverage and a basal layer (c), and a combination of a contact layer with 1.0 ML coverage and a prism layer (d). The contact layer in each case is shown in dark blue and the second layer in light blue. Visual guidelines connecting nearest neighbor water molecules are also shown. (b) and (c) are the most stable structures of all multilayer structures and very closely mimic the prism face (b) or basal face (c) of ice I.

Before discussing the details of different multilayers, we compare monolayer and multilayer ice structures. By adding an additional water layer on top of the contact layer, the adsorption energy (Eads = −0.60 eV/H2O) either increases slightly (compared to the 1.0 ML coverage) or stays the same (compared to the 1.5 ML coverage).

The inset shows that the energetics of the different types of overlayer are very similar and range from Eads = −0.60 eV/H2O (for structures (ii) and (iii)) to Eads = −0.58 eV/H2O (for structure (iv)), (i) being in the middle with an adsorption energy of Eads = −0.59 eV/H2O. These tiny energy differences of a few meV/water are beyond the accuracy of DFT (at the Generalized Gradient Approximation level) for hydrogen bonded systems,[77,78] and we therefore do not attempt to say which structure is most stable. Even though the various multilayer ice structures have similar energies, their structures differ substantially. Some comments on this are useful. First, the second ice layer in the prism+prism structure very closely resembles the prism face of ice I, as can be seen from the red guidelines in Figure b. The average O–O distance of the second ice overlayer in this structure is 2.82 Å, and the average O–O–O angle is 113°. Note that the first overlayer in this case does not adopt a prism face-like structure but instead acts as a mediator between the substrate structure of microcline and the prism face formed in the second overlayer. The ice configuration shown in Figure c also quite closely resembles an ice-like structure, in this case the basal face of ice I (average O–O distance: 2.72 Å, average O–O–O angle: 107°). Recent unbiased molecular dynamics simulations on model systems revealed that the formation of an ice-like overlayer seems to enhance ice nucleating ability.[35] Both of these structures are therefore viable candidates for how ice might form on a clean and defect-free microcline (001) surface. In the other two scenarios, shown in Figure a and 5d, the second layer gets significantly more distorted upon relaxation. In both structures, the in-plane symmetry resembles the hexagonal structure of ice, and the side view reveals however that the buckling is very different compared to a basal and prism face. This makes both of these structures unlikely candidates for explaining what makes feldspar a good ice nucleator.

We now try to understand the structures obtained. From the calculations at 1.5 ML we know that neither basal-like nor prism-like contact layers can form (see Figure e for a reminder). It makes sense that this irregular first layer induces strong distortions in a basal-like overlayer on top of it, simply to make both layers fit onto each other. Not surprisingly, therefore, a multilayer consisting of two basal-like ice layers, as in Figure a, will distort significantly upon relaxation. In contrast, the contact layers found for a coverage of 1.0 ML were much more regular (see Figure d for one particular example). As a consequence, forming a regular basal-like overlayer is much more feasible in this case. This is why the second water layer is able to maintain its basal-like structure in the scenario shown in Figure c. On the other hand, forming a prism-face-like structure on this regular contact layer with 1.0 ML coverage is not favorable. Thus, the second layer in Figure d, which started off as prism-like face, deforms significantly.

To summarize, our results show that 1.0 and 1.5 ML contact layers on microcline can template the formation of an ice-like structure in the second layer. This is interesting because neither the substrate nor the contact layer adopts an ice-like structure. Nevertheless, basal and prismatic ice-like configurations can be formed on top of them. These ice-like arrangements of water molecules become possible in the second overlayer because the first contact layer, even though overall not being very ice-like, still exhibits a local structure akin to that of ice (O–O bond lengths for example are comparable to ice in the first contact layer). In contrast, the underlying substrate, as in the case of feldspar (001), can be very different from an ice-like structure. In other words, the non-ice-like contact layer can be thought of as having a mediating function between the ice-unlike substrate and an ice structure. This finding might help to stimulate thinking in relation to ice nucleation because it greatly enlarges the compound space of potentially efficient ice nucleators. Our finding, based on ab initio calculations, shows that the surface structure of a good ice nucleator need not necessarily to match the structure of ice.

Discussion and Conclusions

Density functional theory combined with a guided structure search approach was used to investigate the feldspar/water interface. At first, we looked at bulk microcline and focused in particular on Al order. We reproduced not only the empirical relation between Al order and lattice dimension but also the fact that microcline has Al highly ordered in its T1 site. Furthermore, our microscopic understanding allowed us to address an issue that was previously overlooked, namely, that there are two rather than one characteristic Si–O bond lengths. As a consequence the 1:1 relation suggested by the empirical model breaks down at the atomistic level.

After establishing a clear picture of the bulk structure, we moved to the (001) surface of microcline. The thermodynamically most favorable order state on the surface was different from the bulk. Instead of Al preferring T1 sites as in the bulk, at the surface they are more stable in T2 sites. This was explained by the greater stability of SiOH groups compared to AlOH groups. Our numbers provide clear support for the well-known experimental fact that Al dissolves prior to Si.[54]

When it came to water adsorption we focused on investigating the configurational space of water on microcline. We found that at low coverages single water molecules are able to bind strongly to the surface. As more water is adsorbed, the adsorption energy per water molecule decreases up to a coverage of 1.0 ML. Contact layers with a coverage of 1.0 and 1.5 ML were found to be ideal candidates to stabilize ice-like structures. Multilayer ice configurations closely resembling the prism face or the basal face of ice I were identified. This was somewhat surprising because neither the microcline surface nor the contact layer resembles an ice-like structure. Nevertheless, ice-like structures can be grown on top of them. Besides feldspar, other substrates such as testosterone[19] or marine biogenic particles[20] are also good ice nucleators without having an apparent match with ice. Our findings here could potentially apply to these systems as well; just because the substrate itself does not resemble the structure of ice it does not automatically follow that ice-like structures cannot be stable on them.

It is interesting to briefly consider how this work could be taken forward. The accuracy of PBE for hydrogen bonded systems is around 50 meV/H2O or better depending on the system.[77,78] The most obvious shortcoming of PBE is that it neglects van der Waals dispersion forces. However, as reported in the SI we have performed a series of tests for feldspar bulk and water covered surfaces which show that the inclusion of dispersion does not alter any of the key conclusions reached in this particular system. The main effect of dispersion is to reduce the energy spacing between stable and metastable structures, in a manner similar to what has been observed for bulk ice structures.[77,78] Given the limitation of PBE, and any exchange-correlation functional for that matter, we have been careful not to attempt to discriminate between low energy structures of similar stability. Rather we have focused on obtaining a general understanding of the low energy structures of ice that are accessible on the surface. Besides approximations coming from the choice of exchange-correlation functional, our study focused exclusively on a perfect (001) surface of microcline at 0 K. Real feldspar particles in our atmosphere will expose a large variety of different surfaces, all of which might contribute to the outstanding ice nucleating ability of this material. Furthermore, the effect of temperature, defects, trenches, wedges, and impurities on its surface remains unclear. Addressing these issues is important for future studies aiming to explain what makes feldspar an excellent ice nucleator.

Even though we did not directly look at ice nucleation on microcline, the DFT data can be useful for guiding and validating future force field studies which might probe the ice nucleating ability of feldspar explicitly. The fact that we now have a good idea about the microcline (001) surface as well as the finding that water does not appear to dissociate on this surface suggests that classical force field simulations are a reasonable choice to investigate this interface further. This is particularly promising because it is becoming feasible to study ice nucleation with unbiased molecular dynamics simulations.[29−31,33−35,76] Furthermore, recent force field studies suggested that the formation of ice-like water overlayers on top of substrates can promote ice nucleation.[35] We therefore envisage future studies, in which a large number of materials could be screened with the method introduced here to identify water structures at complex substrates, for which the ice nucleating ability could subsequently be probed with molecular dynamics simulations. This could make it feasible in the future to rationally design highly efficient ice nucleating agents.