Atomistic View of Mercury Cycling in Polar Snowpacks: Probing the Role of Hg2+ Adsorption Using Ab Initio Calculations

Abstract

Photochemical oxidation of atmospheric elemental mercury (Hg0) promotes reactive oxidized Hg (HgII) adsorption on particles and deposition to the polar snowpack. The deposited Hg either returns to the atmosphere via photochemical reduction or remains in the snowpack depending on the strength of adsorption. In this study, we performed ab initio calculations to understand the atomic-level cause of the fate of adsorbed Hg by determining the adsorption affinity for Hg2+, the simplest form of HgII, of barite, halite, muscovite, illite, and ice-Ih as potential adsorbents. The adsorption affinity was estimated by calculating the energy required to dissociate adsorbed Hg2+ from the adsorbents. The results reveal that Hg2+ is stable on the surfaces of the selected adsorbents, except barite, but is prone to photodissociation under solar ultraviolet radiation. This mild adsorption is expected to contribute to the bidirectional exchange of Hg between the atmosphere and the polar snowpack. Thus, this theoretical approach can provide complementary perspectives on polar Hg dynamics beyond the limitations of field and laboratory experiments. Further studies on more complicated and realistic adsorption models with different HgII species and adsorbent surfaces having diverse defective structures are required to better comprehend air–snow Hg cycling in the polar regions.

1. Introduction

Mercury (Hg) is of significant concern as a global pollutant, and its behavior and fate in the environment have been intensively studied for several decades [1,2,3,4,5,6,7,8,9,10,11]. The snowpack in polar regions serves as a natural sink in the global biogeochemical cycles of trace elements as long as it remains frozen [2]. Unidirectional deposition prevails among most trace elements that reach the polar snowpack by atmospheric transport, but a series of photochemical redox reactions induces bidirectional exchanges of Hg between the atmosphere and the surface snowpack [3,4,5,6]. Gaseous elemental Hg (GEM, Hg⁰) is the predominant species of atmospheric Hg and can be transported over long distances due to its relatively poor reactivity [4,5,6]. When Hg⁰ is oxidized to reactive gaseous mercury (RGM, Hg^II), however, it can be readily incorporated with particles to form particulate Hg (P_Hg), whose deposition causes a substantial increase (decrease) in Hg concentration in the surface snow (near-surface air) [6,8,11]. However, it has been proposed that the deposited Hg mostly returns to the atmosphere via a photochemically induced reduction of P_Hg to Hg⁰ [3,10,12]. Depending on the balance between deposition and reemission, the snowpack alternates its role in Hg cycling between a sink and a source from diurnal to seasonal time scales [11].

The lifetime of P_Hg deposited in surface snow is estimated to be 4–24 h under solar ultraviolet (UV) radiation, but it is observed that a minor fraction of P_Hg remains in the polar snowpack for a prolonged period in resistance to photochemical reemission and is sequestered from air–snow Hg cycling [3,6,10,12]. Since reducible Hg would mostly return to the atmosphere, P_Hg resistant to reemission represents the major Hg form in the deep snowpacks. This fraction can be perennially stored in glaciers or potentially becomes bioavailable to aquatic biota if glaciers melt into the sea [4,6]. The production of P_Hg has not been fully understood, but the adsorption of RGM on snow surfaces might be the most likely mechanism when considering ice as the most abundant adsorbent [1,2,8]. However, Hg concentrations have been observed to correlate with other constituents in glaciers, raising the supposition that certain impurities in ice more strongly adsorb Hg to form P_Hg. Jitaru et al. proposed strong adsorption of RGM on airborne mineral dust particles to explain synchronous increases in Hg and dust concentrations in Antarctic Dome C ice cores during the glacial periods [2]. Such a relationship was uncertain in a recent Antarctic snowpack. The correlation between sea salt Na⁺ and Hg concentrations was rather significant, suggesting the possibility that sea salt aerosols contribute to the formation of the resistant P_Hg [1]. Halides have also been proposed as potential stabilizers of P_Hg in other studies ([1] and other references therein). Laboratory and field experiments have found that the incident angle, and wavelength, intensity of solar radiation, and ambient temperature could affect the photochemical reemission flux [8,9,10,13], but explanations for P_Hg production and the detailed reaction pathways of Hg reemission have been sparse or descriptive in nature. The main difficulties in experimental studies may lie in reproducing the natural air–snow interaction or in changing controlling factors in the laboratory or the field.

Beyond the experimental limitations, advances in the ab initio calculations (i.e., the quantum chemical calculations [14,15,16,17,18]) have enabled exploration of the interactions between Hg and particles all the way down to the atomic level with perspectives on the atomistic origins of P_Hg production. This approach has distinct advantages in that Hg speciation and particles of interest can be flexibly designed and that the calculation results can be mutually compared. For example, a previous study calculated the adsorption energies of Hg⁰ on ice-Ih surfaces with proton defects and showed that defective surfaces more strongly adsorb Hg⁰ than do surfaces without defects [19]. RGM (e.g., Hg²⁺, HgO, HgCl₂, and HgBr₂) are generally regarded to be particle reactive and to facilitate stronger adsorption on particle surfaces, but no direct evidence for the strong adsorption of RGM has been shown. More importantly, there has not been systematic research on which types of particles provide a more effective adsorption surface in the polar snowpack.

The ultimate goal of this study is to gain insight into air–snow Hg cycling at the atomic level by estimating the adsorption affinities for reactive Hg species of potential adsorbents in the polar snowpack. As the first step, we investigate the adsorption of Hg²⁺, the simplest form of RGM, on barite, halite, illite, muscovite, and ice-Ih as potential adsorbents present in the polar snowpack. We compare the relative adsorption affinities by calculating and comparing the binding energies (E_Bind) of Hg²⁺ on those adsorbent surfaces. Since the E_Bind calculated in this study corresponds to the energy required to dissociate Hg²⁺ from the adsorbed surface, it can be presumed that stronger adsorption results in greater E_Bind. Based on the calculated E_Bind, we discuss the potential role of the selected adsorbents in Hg cycling.

2. Methods 2.1. Selected Adsorbents in the Polar Snowpack

We prepared cleaved surfaces of barite, halite, illite, muscovite, and ice-Ih as adsorbents for Hg²⁺. Barite (BaSO₄) is chosen as a representative of volcanic aerosols, and its structure is composed of ionic bonds between two divalent ions (Ba²⁺ and SO₄²⁻) [20,21]. Halite (NaCl) is the major component of sea salt aerosols and has ionic bonds between two monovalent ions (Na⁺ and Cl⁻) [1]. Micaceous phyllosilicates (illite and muscovite) are found as terrestrial aerosols in polar regions and have layered structures [22,23]. The cleaved surface structures of those selected adsorbents are presented in Figure 1, and their structural parameters, including the number of atoms, the thickness of the cleaved surface, and the size of the structural model, are listed in Table 1. Diverse surface structures were prepared for the micaceous phyllosilicates (illite and muscovite) and ice-Ih, in which the cleavage planes are interconnected by relatively weak van der Waals interactions and hydrogen bonds because the surface reactivity considerably changes with the atomic configurations in the uppermost layer (details will be discussed later) [24,25,26,27,28,29]. The potential adsorption sites on the cleaved surfaces are presented in Figure 2.

2.1.1. Cleaved Surface Structure of Barite

The crystal structure of barite, BaSO₄, was adopted from a previous study and cleaved to produce a (001) surface [30]. The cleavage planes of (001) and (210) are known for both natural and synthetic barites [36,37,38,39], but only the (001) surface was examined in this study. In the bulk crystal structure, the negative sulfate ions (SO₄²⁻) are interconnected with the positive Ba ions (Ba²⁺). Unpaired O and Ba atoms on the cleaved surface can serve as reactive adsorption sites for charged adsorbates. The potential adsorption sites are presented in Figure 2a.

2.1.2. Cleaved Surface Structure of Halite

We referred to the crystal structure of halite, NaCl, in the literature [31]. Only the (001) cleaved surface was considered for halite due to its highly symmetrical crystal structure. Unpaired and negatively charged Cl atoms (Cl⁻) on the cleavage plane are potential adsorption sites for positively charged adsorbates. The potential adsorption sites near the Cl atom, which are symmetrically inequivalent, are presented in Figure 2b.

2.1.3. Cleaved Surface Structures of Illite and Muscovite

The bulk crystal structures of illite and muscovite, (K,H₃O)(Al,Mg,Fe)₂(Si,Al)₄O₁₀[(OH)₂,(H₂O)] and KAl₂(AlSi₃O₁₀)(F,OH)₂, respectively, were obtained from the literature [32,33], and they are similar to each other despite their compositional differences (i.e., Si, Mg, Fe, and hydroxide are more abundant in illite than in muscovite) [40]. Their bulk crystal structures are composed of alternating layers of intercalated K atoms and so-called T–O–T layers (i.e., layer composed of a single central AlO₆ octahedral sheet joining two inward SiO₄ tetrahedral sheets on both sides) [32,33,40,41,42,43]. The T–O–T layers are interconnected by intercalated K⁺ cations through relatively weak Coulombic interactions; hence, cleavage along this interlayer plane (001) is the most natural [32,33,40,41,42,43]. Since the surface reactivity of illite is primarily governed by the ordering of the uppermost K atoms (i.e., the first nearest-neighboring atoms) rather than by the Mg and Fe atoms at the central octahedral sites, Mg- and Fe-free illite was used in this study. One-fourth of the Si atoms in the SiO₄ tetrahedral sheets were replaced with Al for both the illite and muscovite structures to balance the charges in the finite-sized surface model structures, although approximately one-sixth of the Si atoms are usually substituted by Al atoms in natural illite containing Mg and Fe [40]. The (001) cleaved surface structures of illite and muscovite were prepared to determine the possibility of Hg²⁺ adsorption on the illite and muscovite surfaces. While the ordering of the K atoms on the cleavage planes has not been directly revealed yet [24,25], previous studies suggested that the intercalated K atoms may be statistically shared on the cleavage planes directly above the T–O–T layers [25,26,44,45]. Accordingly, several types of (001) cleaved surfaces with different orderings of the uppermost K atoms were prepared for both illite and muscovite. The potential adsorption sites near the ditrigonal siloxane cavities are shown in Figure 2c,d.

2.1.4. Cleaved Surface Structure of Ice-Ih

The ice-Ih, which is the so-called hexagonal ice because of its atomic configuration, has been known as the primary phase in natural snow and glacier ice [46]. The crystal structure of ice-Ih was chosen from a previous study to produce cleaved surfaces [34]. Several cleavage planes exist for the ice-Ih (i.e., the (0001) basal plane, the (10 1¯ 0) primary prism plane, the (11 2¯ 0) secondary prism plane, and the (20-21) pyramidal plane), but compared to the other cleavage planes, the basal and primary prism planes are easily accessible in experiments [47]. The basal plane and the primary prism plane refer to the top and bottom faces and the six equivalent side faces of the hexagonal column (i.e., going anticlockwise (01 1¯ 0), ( 1¯ 100), ( 1¯ 010), (0 1¯ 10), (1 1¯ 00), and (10 1¯ 0) directions) [47]. We selected the (0001) basal plane for the adsorption affinity calculation. Previous studies have suggested that the (0001) cleaved surface is usually decorated with quasi-randomly distributed dangling H and O atoms (i.e., the lone pairs of hydrogen bonds) originating from either the termination of the top half-bilayer or the orientational disorder of H₂O molecules in the hydrogen-bonded network (e.g., the Bjerrum L-type defect) [27,28,29,35,48,49,50,51,52,53]. It has been revealed that the dangling H and O atoms on the cleaved ice-Ih surface are significant for the surface energy and the chemical reactivity [27,28,29,48]. Those dangling H and O atoms and the defects can act as reactive adsorption sites for charged adsorbates by forming hydrogen-bond-like electrostatic bonds. We prepared diverse cleavage surfaces with different orderings of H and O atoms to determine how the dangling atoms and the defects affect the adsorption affinity of ice-Ih. Dangling O atoms can be produced by termination of the half-bilayer [51], but the (000 1¯ ) cleaved surface was used instead since it was already decorated with dangling O atoms without structural modification. Therefore, both the (0001) and (0001) cleaved surfaces were used to probe the contributions of dangling H and O atoms to the adsorption affinity of ice-Ih. The surface reactivity is also affected by proton defects on the cleaved plane associated with proton hopping (i.e., the Grotthuss mechanism) [54]. For the defect structure, several types of (0001) cleaved surfaces with distinct proton defects (i.e., unoccupied atomic sites for dangling H atoms) were prepared. The details of the cleaved surface structures of ice-Ih are presented in Figure 1e, and the potential adsorption sites are indicated in Figure 2e.

2.1.5. Constructing the Cleaved Surface Structures

Structure optimization calculations were applied to the adopted bulk crystal structures from References [30,31,32,33,34] to obtain more precise crystal structures. Each optimized crystal structure was extended to a supercell structure by repeating the unit cell, the symmetrically simplest repetitive unit in the crystal structure, in a certain direction. These extended supercell structures were cleaved along naturally occurring directions. Sufficiently large vacuum spaces of more than ≈40 Å were included above the cleaved surfaces to avoid unexpected interatomic interactions between the periodically repeated layered cleaved surface structures (i.e., because of the periodic boundary condition) [17,55]. The structure optimization calculations were carried out again for these cleaved surface structures because the interlayer distances and atomic structures near the cleavage planes tend to slightly differ from those in their initial bulk structures.

2.2. Binding Energy Calculation

E_Bind of Hg²⁺ to an adsorbent can be calculated using the following equation [56]:

E_Bind = (E_CD − E_AB) − (E_C − E_A) − (E_D − E_B)

E_AB and E_CD refer to the calculated total energies for the adsorbent with Hg²⁺ in the initial (i.e., Hg²⁺ is far enough from the cleaved surface not to interact with it) and final (i.e., Hg²⁺ is stabilized at a potential adsorption site on the cleaved surface) states of adsorption, respectively. Then, (E_CD − E_AB) is the total energy difference caused by Hg²⁺ adsorption on the cleaved surface of the adsorbent. (E_C − E_A) indicates the energy difference associated with configurational changes in the cleaved surfaces during Hg²⁺ adsorption, in which E_A and E_C refer to the calculated energies for the cleaved surface structures before and after Hg²⁺ adsorption, respectively. The energy difference in Hg²⁺ before (E_B) and after (E_D) adsorption, which originates from the theoretical limitation of using a finite number of basis sets to describe the electron density distribution in the finite-sized structural model, was relatively negligible but was considered in the E_Bind calculation [14,55]. Then, E_Bind corresponds to the energy required to remove Hg²⁺ from the cleaved surface of an adsorbent. Accordingly, a stronger adsorption affinity for Hg²⁺ of an adsorbent will result in a more negative E_Bind. A negative E_Bind implies that Hg²⁺ is stabilized by physiochemical adsorption on an adsorbent. In contrast, a non-negative E_Bind indicates that Hg²⁺ is unstable on a potential adsorbent. The repulsive interatomic interaction between Hg²⁺ and a cleaved surface would result in a non-interacting system in which Hg²⁺ is separated far enough from the cleaved surface that the theoretical maximum E_Bind would be zero. On the other hand, unexpected positive binding energy (E_Bind > 0) can occur when the adsorbate is trapped in a metastable state rather than remaining in the globally minimized non-interacting system (E_Bind = 0). Usually, the E_Bind of a specific adsorbate can be simplified to E_CD − (E_C + E_D) if the adsorbate is neutrally charged and the energy difference due to the charge redistribution of the adsorbate is negligible (E_AB ≈ E_A + E_B). However, this precise E_Bind calculation, in which the binding energy is estimated by the direct comparison between the initial and final states of Hg²⁺ adsorption, can be applied to the charged system because it can ignore the effect of charge redistribution in the finite-sized structural model, imposing a significant energy difference (i.e., 15–20 eV in the size of the structural model used in this study).

2.3. Structure Optimization and Electronic Structure Calculations

Single-point energy calculations and structure optimization calculations for the selected adsorbents were carried out using the ab initio calculations [14]. Single-point energy calculations, which provided the total energy of the surface structure model, were carried out to determine the binding energy of Hg²⁺ on the cleaved surfaces. The interatomic interactions were represented with the on-the-fly (OTF) ultrasoft pseudopotential [14]. The on-site electron–electron interactions were determined using the generalized gradient approximation (GGA)-based Perdew–Burke–Ernzerhof (PBE) exchange–correlation functional [57]. The long-range order dispersion interactions, such as the van der Waals interaction, significantly influence the molecular crystal structures and cleaved surfaces with hydrogen bond-like long range interactions [58,59,60,61]. Accordingly, the hybrid semi-empirical dispersion correction term based on the Tkatchenko–Scheffler method was applied to the single-point energy calculations and the structure optimization calculations [58,59]. Since the structure optimization calculations were more time consuming, they were performed with slightly rough convergence criteria for calculation efficiency. The structure optimization calculations used the two-point steepest descent (TPSD) algorithm to search for the stable atomic configurations [62]. The convergence criteria of the structure optimization were set to 2 × 10⁻⁵ eV/atom, 0.05 eV/Å, 0.1 GPa, and 0.002 Å for the total energy, maximum forces, maximum stress, and maximum displacement, respectively. The electronic structures and the interatomic interactions were calculated with the following parameters. The plane-wave-cutoff energy, used to define the precision of the electronic wave functions, was set to 500 eV. The convergence criterion for the self-consistent field (SCF) calculations was set to 2 × 10⁻⁶ eV/atom. Because a more significant number of k points could enhance the calculation accuracy, the single Γ point in the k point grid was used for calculation efficiency. The single-point energy calculations were carried out more accurately than the structure optimization calculations to enhance the accuracy of the calculated total energy values. A higher plane-wave-cutoff energy of 600 eV was used, and a smaller convergence criterion for the SCF calculations of 1 × 10⁻⁶ eV/atom was used. For heavy metal elements including Hg, it is required to consider the relativistic effects when the detailed electronic structures associated with spin-orbital coupling are to be probed [63,64,65]. However, introducing a relativistic treatment will cause only slight changes in the calculation results for the simple adsorption. In this study, the commonly used ultrasoft pseudopotential based on the non-relativistic Schrodinger approach was applied to increase calculation efficiency [66].

3. Results and Discussion

Figure 3 shows the atomic configurations of Hg²⁺ adsorbed on potential adsorption sites on the selected adsorbents. The electrostatic charge of each atom around the bounded Hg (Hg²⁺) is labeled to represent the charge redistribution resulting from the self-consistent minimization calculations of the electron densities. The charges of the bounded Hg were 0.46–1.12 lower than the point charge of 2. The calculated E_Bind values of Hg²⁺ at the potential adsorption sites are given in Table 2 and compared between the selected adsorbents in Figure 4.

3.1. Hg²⁺ on Barite

When each Hg²⁺ was placed at the four potential adsorption sites on the (001) cleaved surface of barite (Figure 2a), the E_Bind values were calculated to be non-negative in the range of 1.0154–1.8952 eV (Table 2) and plotted in the unstable regime in Figure 4 for all cases. These results (E_Bind > 0) suggest that Hg²⁺ could not be stabilized on the barite surface. Although the result in Figure 3a presents covalent-like bonds between the Hg²⁺ and neighboring O atoms, this atomic configuration was likely in a metastable state as indicated by the non-negative E_Bind. Therefore, the adsorption affinity for Hg²⁺ of barite could be concluded to be poor. Because of the improbable adsorption of highly reactive Hg²⁺ on the barite surface, Hg might only coexist with barite as fluid inclusions or by substitution into the crystal lattice ([67,68] and references therein). For example, barite has been used as a weighting agent of drilling fluid in petroleum engineering due to its chemical stability, and it compatibly contains Hg as an impurity [67,68]. Barite may not have a significant role in Hg cycling in the polar snowpack because a chemical reaction between them is unlikely. If correlations are found among Ba, SO₄, and Hg in the polar snowpack, preexisting Hg in barite should be carefully taken into consideration.

3.2. Hg²⁺ on Halite

The two potential adsorption sites on the (001) cleaved surface of halite, directly above the Na and Cl atoms (Figure 2b), resulted in distinctly different binding energies. The E_Bind near Na was 0.5771 eV (Table 2) and plotted in the unstable regime in Figure 4 (E_Bind > 0) due to the locally applied repulsive interactions arising from the positively charged Na. In contrast, Hg²⁺ appeared to be stabilized in the form of HgCl₂ near the negatively charged Cl atoms (Figure 3b), with a negative E_Bind of −0.7535 eV (Table 2). These results suggest that the adsorption affinity for positively charged Hg²⁺ of the halite surface is dependent on the neighboring atoms, and thus, halite can stabilize Hg²⁺ on its surface by facilitating the formation of HgCl₂. The favorable adsorption of Hg²⁺ on halite is consistent with the correlation observed between the Hg and the sea salt concentrations in an Antarctic Plateau snowpack [1]. In this case, the supposition that the sea salts have enhanced Hg deposition in the polar regions seems to be reasonable at the atomic level [1].

3.3. Hg²⁺ on Illite and Muscovite

The considered adsorption surfaces of illite and muscovite are shown in Figure 2c,d with the different orderings of the uppermost K atoms. The calculated E_Bind at potential adsorption sites of illite and muscovite ranged from −1.3098 to −0.1390 eV and from −2.0415 to −0.8109 eV, respectively. For both illite and muscovite, reactive adsorption sites were observed directly above the empty ditrigonal siloxane cavities (i.e., sites lacking the uppermost K atoms; type-A1, A2, and B in Table 2). At these sites, Hg²⁺ formed covalent-like bonds with the adjacent O atoms (Figure 3c,d). Unexpectedly, the E_Bind was lower when the uppermost K sites were fully unoccupied (type-B in Table 2) rather than partially occupied (type-A1 and A2), as presented in Table 2 and Figure 4. This result occurred because the structural distortions of the empty siloxane cavities (i.e., the O atoms at the edge of the ditrigonal siloxane cavities move slightly toward the center, as shown in Figure 2c,d) could mitigate the structural changes in the cleavage planes during Hg²⁺ adsorption and the corresponding energy differences (see E_C − E_A and E_Bind in Table 2). Despite the presence of the positively charged uppermost K atoms, Hg²⁺ could be stable on the cleaved surfaces even with the fully occupied K sites. This result seems to be caused by the K atoms compensating for the electron-deficient Hg²⁺ during the charge redistribution (the charge of the stabilized Hg²⁺ was reduced to 0.02 and 0.01 for illite and muscovite, respectively), resulting in a reduction in the repulsive interatomic interactions. Overall, illite and muscovite are capable of stabilizing Hg²⁺ on their surfaces and even on the cleaved surface fully decorated with positively charged K atoms. Adsorption of heavy metals to clays and associated kinetics have been intensively investigated because of their significance toward the disposal of persistent metallic pollutants in aqueous environmental systems ([69,70] and other references therein). For instance, the adsorption affinity for Hg^II of micaceous phyllosilicates, such as montmorillonite, muscovite, and illite, in the presence of coexisting substances in the aqueous phase (e.g., fulvic and other acids and dissolved organic matter) has been explored [71,72,73,74,75,76], revealing that the clay minerals can effectively adsorb Hg^II in solution. Given that coexisting substances were not considered in this study, which is more relevant to Hg cycling between the polar atmosphere and snowpack [22,23], our results further suggest that direct adsorption of Hg^II on the clay minerals is probable.

3.4. Hg²⁺ on Ice-Ih

The (0001) cleaved surfaces of ice-Ih consisting of configurationally different uppermost bilayers are shown in Figure 2e and Figure A1. The calculated E_Bind for Hg²⁺ on ice-Ih ranged from −4.0785 to 0.0442 eV (Table 2). The E_Bind on the surface with fully ordered dangling H atoms (type-F in Table 2) had a non-negative value, implying that Hg²⁺ is slightly unstable near the positively charged dangling H atoms (Figure A1a). In contrast, the E_Bind values of Hg²⁺ on the (000 1¯ ) cleaved surface with fully ordered dangling O atoms (type-R in Table 2) are revealed to be in the stable regime in Figure 4, suggesting Hg²⁺ adsorption near the negatively charged dangling O atoms. In Figure 3e, for example, Hg²⁺ formed multiple covalent-like bonds with the surrounding dangling O atoms. Hg²⁺ was more stable when captured by multiple dangling O atoms (i.e., Figure 3e; R(01) in Figure 4) than when placed directly above a single dangling O atom (i.e., Figure A1f; R(02) in Figure 4). The proton defects on the (0001) surface (type-A1, A2, and B in Table 2) were revealed to be potential reactive adsorption sites (E_Bind < 0), at which Hg²⁺ formed covalent-like bonds with nearby unpaired O atoms (Figure A1b–e). Like those of illite and muscovite, structural distortions around the proton defects were shown to affect the surface reactivity of ice-Ih: The greater the number of proton defects present, the greater the structural distortion of the uppermost bilayer, mitigating the structural distortion associated with Hg²⁺ adsorption. Consequently, Hg²⁺ entirely decorated with proton defects became less stable than that partially decorated with proton defects. If all cases are taken together, the results show that ice-Ih is capable of stabilizing Hg²⁺ on its surface despite the wide range of E_Bind determined from the surface atomic configurations.

To understand chemical processes in polar environments or interstellar media, unusual chemical reactions on ice-Ih surfaces have been studied with sophisticated experimental techniques [50,77,78,79,80,81,82]. Theoretical approaches using the ab initio calculations have also been applied to the adsorption of specific adsorbents on ice-Ih surfaces, such as acids (HOCl, HCOOH, and CH₃COOH), halide ions (F⁻, Cl⁻, and Br⁻), alkali metals (Na and Na⁺), and heavy metals (Hg⁰) [19,49,83,84]. In these studies, the dangling H and O atoms and the proton defects were revealed to be the reactive adsorption sites on the (0001) cleaved surfaces because of their electronic instability. In particular, defective ice-Ih was shown to be able to adsorb Hg⁰ [19], although the binding energy was calculated to be less negative (from −0.31 to −0.14 eV; [19]) than E_Bind for Hg²⁺ obtained in this study (from −4.08 to −0.36 eV). The lower E_Bind for Hg²⁺ reflects relatively strong interatomic interactions between Hg²⁺ and the ice-Ih surfaces. Thus, it can be concluded that due to its highly reactive characteristics, Hg²⁺ is considerably involved in Hg cycling in polar snowpack and fractional deposition.

3.5. Relative Adsorption Affinities to the Selected Adsorbents

Comparison of the E_Bind values provides an opportunity to address the relative adsorption affinity between Hg²⁺ and the selected adsorbents (Figure 4 and Table 2). While most of the selected adsorbents were revealed to be able to stabilize Hg²⁺ on their surfaces (E_Bind < 0), ice-Ih seems to be the most effective adsorbent. However, caution is needed when concluding that natural ice is superior to the other surfaces for Hg²⁺ adsorption because the defective surface structures of ice-Ih, particularly the proton defects formed by proton hopping [19,77,78,82,85] and the dangling O atoms formed by half-bilayer termination [51], can not only enhance adsorption affinity but also increase structural instability [51,78,85]. In the natural environment, the occurrence of such reactive but energetically unstable surfaces can therefore be less probable. Thus, the adsorption affinity for Hg²⁺ of natural ice, likely decorated with fewer structural defects, might be smaller than the ice-Ih that was decorated full of structural defects estimated here, implying fewer interatomic interactions. Likewise, the adsorption affinity of illite and muscovite can be smaller than their estimated maxima depending on the natural abundance of unoccupied uppermost K sites. In summary, Figure 4 provides a comparison of the range of E_Bind values between the adsorbents and shows that all the adsorbents, except barite, can stabilize Hg²⁺ on their surfaces; however, the relative effectiveness of adsorption is determined by the abundance of reactive adsorption sites on the adsorbent surface.

3.6. Implication for Hg Cycling between the Polar Atmosphere and Snowpack

Adsorption and dissociation can change the form of Hg between reactive gaseous Hg^II and P_Hg, and hence have been considered the mechanisms that likely promote bidirectional Hg exchange between the atmosphere and the surface snow [3,4,10,12]. In this regard, the calculated adsorption affinities of adsorbents provide quantitative support for the conceptual understanding of the role of adsorption and dissociation in Hg cycling, thus expanding the atomistic-level discussion to the natural implications. Our results support that snow and marine and terrestrial aerosols contribute to the formation of P_Hg by capturing atmospheric Hg²⁺ on their surfaces in polar regions. The lifetime (e-folding decay time) of P_Hg in the surface snow after deposition, estimated via a laboratory experiment [3,4,10,12], was 4–24 h under UV radiation in an effective wavelength range from 280–400 nm. The energy of the photons corresponding to the wavelength range is within 3.10–4.43 eV and comparable to the calculated E_Bind (Table 2 and Figure 4). This finding suggests the slight chance that adsorbed Hg²⁺ may be photodissociated by UV irradiation to participate in air–snow Hg cycling. Conventional ab initio calculations are known to sometimes underestimate the actual E_Bind [86], but our results offer a rather reliable estimation of E_Bind since the semi-empirical dispersion correction was applied to account for the non-local van der Waals interactions between the approaching Hg²⁺ and the cleaved surfaces [59,86]. The halite results suggest the formation of HgCl₂ molecules by Hg²⁺ adsorption on the cleaved surface (Figure 3b). HgCl₂ is found as one of reactive gaseous Hg^II species over the polar snowpack [3,4,5,13], but direct dissociation of HgCl₂ from a pristine halite surface seems to be unlikely because of the high energy required to dissociate two bound Cl atoms from the surface (approximately −4.26 eV each). Nevertheless, we cannot exclude the possibility of HgCl₂ liberation from atomic and molecular defects where the structural instability is high. In this case, Hg cycling involves HgCl₂, and halite represents a source of chlorine for the near-surface atmosphere.

Despite intense reemission of deposited Hg, a fraction of Hg is observed to be sequestered from Hg cycling and to remain perennially in the snowpack [1,8,10]. It is still unclear what controls the degree of Hg sequestration. Although our results do not provide a decisive clue, at least the P_Hg discussed in this study seems to rarely contribute to Hg sequestration, as indicated by the relatively smaller E_Bind than expected. In addition, E_Bind values of illite, muscovite, and halite were less negative than those of ice-Ih, and hence, their stabilizing effects are unlikely to be superior to those of ice-Ih. This interpretation is contrary to the expectations of previous studies in which clay minerals and sea salts were anticipated to act as stabilizing agents for Hg deposited in the polar snowpack. Therefore, the enhancement in Hg sequestration and the consequent increase in Hg concentration in the snowpack should be discussed further, taking into consideration both other adsorbates and adsorbents with greater E_Bind and external causes, such as heavy snowfall burying deposited Hg below the sunlit layer [13].

3.7. Limitation and Future Study

Despite the efforts in this study, there are limitations to be considered in discussing the implications for natural phenomena, especially regarding cleaved surface structures, Hg compounds, and temperature.

For the selected adsorbents, several types of cleaved surfaces were explored to cope with the diverse defective structures. However, the structural defects of the barite and halite surfaces were ignored for calculation efficiency. In addition, although we prepared surface structures of illite, muscovite, and ice-Ih with different structural defects, more complicated structural defects and disorders than those considered in this study would also exist. For example, ice-Ih has been known to have ionic defects such as hydronium (H₃O⁺) and hydroxide (OH⁻) induced from structural disorders as follows [77,78,79,80,81,85]: proton disorder associated with molecular disorientation of H₂O molecules at a finite temperature (>72 K) within the constraints of the Bernal–Fowler–Pauling ice rules [27,28,29,87], quasi-randomly distributed dangling H and O atoms formed by the full- or half-bilayer terminations [51,82], molecular point defects due to sublimation of H₂O molecules [53,88,89], proton defects caused by the Grotthuss mechanism [77,78,82], and a quasi-liquid layer arising from thermally induced molecular disorder [50,77,79,82,90,91,92].

Air–snow Hg cycling would include diverse Hg compounds, such as adsorbates [6,8,10], which were not addressed in this study. Reactive gaseous Hg compounds (e.g., HgO, HgCl₂, and HgBr₂) have been expected to readily adsorb on particulates and be deposited on the surface snow [3,4,5,13]. Therefore, the adsorption and dissociation of these Hg compounds should be analyzed in future studies to comprehensively demonstrate Hg cycling. Coexisting substances need further consideration since they directly or indirectly influence the behavior of Hg with competitive or cooperative interactions. The presence of reducing agents in snow (such as benzophenone) can cause Hg^II reduction to Hg⁰, which is likely to be reemitted from the adsorbed particle surface [10].

The conventional ab initio calculations used in this study provide only the energy of the electronic contribution at 0 K [93]. Because adsorption of Hg compounds on ice surfaces in polar regions occurs at finite temperatures ( > 200 K), thermally induced structural disorder [50,77,79,82,90,91,92], hopping of excess protons [77,78,82], and vibrational energy and entropy [9] should be explored in future studies to obtain a comprehensive understanding of the consequent changes in the surface reactivity of ice and the Hg adsorption rate.

4. Conclusions

Adsorption has been postulated to play a significant role in the air–snowpack exchange of Hg in polar regions but has hardly been investigated in laboratory and field experiments. The present study was designed to address the role of adsorption at the atomic scale by using the ab initio calculation and to discuss its implication for Hg cycling. As a first step of the calculational approach, relatively simple adsorption processes of Hg²⁺ on barite, halite, illite, muscovite, and ice-Ih were explored to evaluate the adsorption affinity of the selected adsorbents. The results revealed that Hg²⁺ can be stabilized on the surfaces of halite, illite, muscovite, and ice-Ih (but not barite), forming P_Hg and facilitating Hg deposition. The calculated binding energies of Hg²⁺ were comparable to the photon energy of UV irradiation (280–400 nm wavelength), suggesting the possibility of direct photodissociation of adsorbed Hg²⁺ to participate in air–snow Hg cycling. Overall, our results highlight that mild adsorption of Hg²⁺ on ice and aerosol particles, and the subsequent photodissociation, can contribute to dynamic Hg cycling between the near-surface atmosphere and the surface snow in polar regions. Nevertheless, the current results do not seem to support the previous supposition that sea salts and mineral dust particles adsorb reactive gaseous Hg more than ice particles and promote Hg sequestration in the polar snowpack; rather, ice-Ih had a stronger adsorption affinity than halite, illite, and muscovite. This unexpectedly significant role of ice-Ih in air–snow Hg cycling needs to be further investigated in future studies, including those addressing diverse reactive gaseous Hg species as adsorbates and more varied surface structures of adsorbents. We claim that this systematic theoretical approach can provide a unique perspective on the behavior of Hg across the interface of the atmosphere and polar snowpack at the atomistic scale.

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

Crystal	Barite	Halite	Illite	Muscovite	Ice-Ih
System	Orthorhombic	Cubic	Monoclinic	Monoclinic	Monoclinic
Space group	(62) PNMA	(225) FM3¯ M	(12) C2/M	(15) C2/C	(9) CC
Lattice parameters of the unit cell
a (Å)	8.8842	5.6400	5.2021	5.2108	7.8219
b (Å)	5.4559	5.6400	8.9797	9.0399	8.6299
c (Å)	7.1569	5.6400	10.2260	20.0210	8.6299
α (°)	90.00	90.00	90.00	90.00	90.00
β (°)	90.00	90.00	95.76	101.57	121.55
γ (°)	90.00	90.00	90.00	90.00	90.00
Supercell sizes used to construct surface structure models
N × N × N	2 × 2 × 2	3 × 3 × 3	2 × 2 × 1	2 × 2 × 1	3 × 3 × 3
Direction of cleavage plane
h k l	(001)	(001)	(001)	(001)	(0001)
References
	[30]	[31]	[32]	[33]	[34,35]
Cleaved surface structures
a (Å)	17.9322	16.1735	10.4042	10.4233	13.1921
b (Å)	10.8498	16.1735	17.9594	18.1650	22.8061
c (Å)	50.1536	48.0867	48.4581	43.5444	52.6997
Area of surface (Å2)	194.5606	261.5820	186.8532	189.3395	300.8610
Thickness of layer (Å)	10.1882	8.2059	7.9319	7.8718	13.0837
Vacuum thickness (Å)	39.9655	39.8808	40.5262	35.6725	39.6160
Number of atoms (depending on atomic configurations)
	140	196	160–168	160–168	414–432

Energy (eV)	ECD − EAB 1	EC − EA 2	ED − EB 3	EBind4
Barite
Site 01	1.9653	0.0699	0.0002	1.8952
Site 02	1.7905	0.1160	0.0001	1.6744
Site 03	1.5100	0.4945	0.0002	1.0154
Site 04	1.7683	0.1315	0.0002	1.6367
Site 05	1.8801	0.0286	0.0003	1.8513
Halite
Site 01	0.3853	1.1386	0.0001	−0.7535
Site 02	0.6281	0.0511	−0.0001	0.5771
Illite
Type F site 01	0.8477	0.9866	0.0001	−0.1390
Type A1 site 01	0.0675	1.3775	−0.0001	−1.3098
Type A2 site 01	0.2617	1.3778	−0.0001	−1.1160
Type B site 01	0.0844	0.8579	−0.0002	−0.7732
Type B site 02	0.2820	0.8387	0.0001	−0.5568
Muscovite
Type F site 01	0.1645	0.9753	0.0002	−0.8109
Type A1 site 01	−0.5770	1.4648	−0.0004	−2.0415
Type A2 site 01	−0.2466	1.4187	−0.0001	−1.6653
Type B site 01	−0.5237	0.8904	−0.0003	−1.4138
Type B site 02	−0.2129	0.8790	−0.0001	−1.0918
Ice-Ih5
Type F site 01	0.0440	−0.0001	−0.0001	0.0442
Type R site 01	−2.2870	1.1908	0.0000	−3.4778
Type R site 02	−1.8707	0.7601	0.0001	−2.6309
Type A1 site 01	−3.0650	1.0133	0.0002	−4.0785
Type A1 site 02	−0.2325	0.1247	−0.0001	−0.3571
Type A2 site 01	−2.7403	0.9509	0.0003	−3.6916
Type B site 01	−2.7651	0.6175	−0.0002	−3.3824

1 Energy difference between the final and the initial states of Hg2+ adsorption. 2 Total energy change due to the configuration change of the cleaved surface structure caused by adsorption. 3 Small energy change induced by the movement of Hg2+ during adsorption, owing to the methodological imperfectness of the ab initio calculation using a finite number of basis sets to describe the electronic orbitals in a finite-size surface structure model. 4 Calculated binding energy of Hg2+; a more negative binding energy means a stronger interatomic interaction between Hg2+ and the adsorbent surface. 5 Type F, A1, A2, and B structures indicate the (0001) cleaved surfaces, and the type R structure indicates the (0001¯ ) cleaved surface of ice-Ih.

Author Contributions

Y.S.Y., Y.H., and S.D.H. conceived the idea; Y.S.Y. designed the computational approaches and performed the ab initio calculations; S.K.L. provided the computational resources; Y.S.Y. and Y.H. wrote the manuscript.

Funding

This research was supported by research grants (PE19040 and PE19200) from the Korea Polar Research Institute (KOPRI) and the NRF, Korea (2017R1A2A1A17069511).

Acknowledgments

We are grateful to three anonymous reviewers for their comments.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Figure A1.Perpendicular views of the Hg2+ (yellow dashed circle) adsorbed on different types of ice-Ih surfaces: (a) type F site 01, (b) type A1 site 01, (c) type A1 site 01, (d) type A2 site 02, (e) type B site 01, and (f) type R site 02. Only the atoms in the uppermost layers are depicted.

Figure A1.Perpendicular views of the Hg2+ (yellow dashed circle) adsorbed on different types of ice-Ih surfaces: (a) type F site 01, (b) type A1 site 01, (c) type A1 site 01, (d) type A2 site 02, (e) type B site 01, and (f) type R site 02. Only the atoms in the uppermost layers are depicted.

References

1. Han, Y.; Huh, Y.; Hur, S.D.; Hong, S.; Chung, J.W.; Motoyama, H. Net deposition of mercury to the Antarctic Plateau enhanced by sea salt. Sci. Total Environ. 2017, 583, 81–87.

2. Jitaru, P.; Gabrielli, P.; Marteel, A.; Plane, J.M.C.; Planchon, F.A.M.; Gauchard, P.A.; Ferrari, C.P.; Boutron, C.F.; Adams, F.C.; Hong, S.; et al. Atmospheric depletion of mercury over Antarctica during glacial periods. Nat. Geosci. 2009, 2, 505.

3. Dommergue, A.; Bahlmann, E.; Ebinghaus, R.; Ferrari, C.; Boutron, C. Laboratory simulation of Hg0 emissions from a snowpack. Anal. Bioanal. Chem. 2007, 388, 319–327.

4. AMAP. AMAP Assessment 2011: Mercury in the Arctic Arctic Monitoring and Assessment Programme (AMAP); AMAP: Oslo, Norway, 2011; p. 193.

5. Lindberg, S.E.; Brooks, S.; Lin, C.J.; Scott, K.J.; Landis, M.S.; Stevens, R.K.; Goodsite, M.; Richter, A. Dynamic oxidation of gaseous mercury in the arctic troposphere at polar sunrise. Environ. Sci. Technol. 2002, 36, 1245–1256.

6. Steffen, A.; Douglas, T.; Amyot, M.; Ariya, P.; Aspmo, K.; Berg, T.; Bottenheim, J.; Brooks, S.; Cobbett, F.; Dastoor, A.; et al. A synthesis of atmospheric mercury depletion event chemistry linking atmosphere, snow and water. Atmos. Chem. Phys. Discuss. 2007, 7, 10837–10931.

7. Bartels-Rausch, T.; Krysztofiak, G.; Bernhard, A.; Schläppi, M.; Schwikowski, M.; Ammann, M. Photoinduced reduction of divalent mercury in ice by organic matter. Chemosphere 2011, 82, 199–203.

8. Douglas, T.A.; Sturm, M.; Simpson, W.R.; Blum, J.D.; Alvarez-Aviles, L.; Keeler, G.J.; Perovich, D.K.; Biswas, A.; Johnson, K. Influence of snow and ice crystal formation and accumulation on mercury deposition to the arctic. Environ. Sci. Technol. 2008, 42, 1542–1551.

9. Bartels-Rausch, T.; Huthwelker, T.; Jöri, M.; Gäggeler, H.W.; Ammann, M. Interaction of gaseous elemental mercury with snow surfaces: Laboratory investigation. Environ. Res. Lett. 2008, 3, 045009.

10. Mann, E.A.; Mallory, M.L.; Ziegler, S.E.; Avery, T.S.; Tordon, R.; O’Driscoll, N.J. Photoreducible mercury loss from arctic snow is influenced by temperature and snow age. Environ. Sci. Technol. 2015, 49, 12120–12126.

11. Angot, H.; Magand, O.; Helmig, D.; Ricaud, P.; Quennehen, B.; Gallée, H.; Del Guasta, M.; Sprovieri, F.; Pirrone, N.; Savarino, J.; et al. New insights into the atmospheric mercury cycling in central Antarctica and implications on a continental scale. Atmos. Chem. Phys. 2016, 16, 8249–8264.

12. Lalonde, J.D.; Poulain, A.J.; Amyot, M. The role of mercury redox reactions in snow on snow-to-air mercury transfer. Environ. Sci. Technol. 2002, 36, 174–178.

13. King, M.D.; Simpson, W.R. Extinction of UV radiation in Arctic snow at Alert, Canada (82 °N). J. Geophys. Res. Atmos. 2001, 106, 12499–12507.

14. Clark, S.J.; Segall, M.D.; Pickard, C.J.; Hasnip, P.J.; Probert, M.J.; Refson, K.; Payne, M.C. First principles methods using CASTEP. Z. Krist. Cryst. Mater. 2005, 220, 567–570.

15. Tse, J.S. Ab Initio molecular dynamics with density functional theory. Annu. Rev. Phys. Chem. 2002, 53, 249–290.

16. Marx, D.; Hutter, J. Ab initio molecular dynamics: Theory and implementation. Mod. Methods Algorithms Quantum Chem. 2000, 1, 301–449.

17. Payne, M.C.; Teter, M.P.; Allan, D.C.; Arias, T.A.; Joannopoulos, J.D. Iterative minimization techniques for ab initio total-energy calculations: Molecular dynamics and conjugate gradients. Rev. Mod. Phys. 1992, 64, 1045–1096.

18. Pastore, G.; Smargiassi, E.; Buda, F. Theory of ab initio molecular-dynamics calculations. Phys. Rev. A 1991, 44, 6334–6347.

19. Asaduzzaman, A.M.; Schreckenbach, G. Adsorption of Na and Hg on the ice (Ih) surface: A density-functional study. J. Phys. Chem. C 2010, 114, 2941–2946.

20. Li, Y.H.; Schoonmaker, J.E. Chemical composition and mineralogy of marine sediments. In Treatise on Geochemistry; Holland, H.D., Turekian, K.K., Eds.; Pergamon: Oxford, UK, 2003; pp. 1–35.

21. Obenholzner, J.; Schroettner, H.; Delgado, H. Barite aerosol particles from volcanic plumes and fumaroles-FESEM/EDS analysis. In Proceedings of the EGS-AGU-EUG Joint Assembly, Nice, France, 6–11 April 2003; p. 8119.

22. Gaudichet, A.; De Angelis, M.; Lefevre, R.; Petit, J.; Korotkevitch, Y.; Petrov, V. Mineralogy of insoluble particles in the Vostok Antarctic ice core over the last climatic cycle (150 kyr). Geophys. Res. Lett. 1988, 15, 1471–1474.

23. Sala, M.; Dapiaggi, M.; Delmonte, B.; Marino, F.; Artioli, G.; Maggi, V.; Revel-Rolland, M.; Petit, J. Mineralogical composition of EPICA Dome C aeolian ice core dust. In Proceedings of the EGU General Assembly, Vienna, Austria, 19–24 April 2009; p. 6230.

24. Kuwahara, Y. Muscovite surface structure imaged by fluid contact mode AFM. Phys. Chem. Miner. 1999, 26, 198–205.

25. Ostendorf, F.; Schmitz, C.; Hirth, S.; Kühnle, A.; Kolodziej, J.J.; Reichling, M. How flat is an air-cleaved mica surface? Nanotechnology 2008, 19, 305705.

26. Müller, K.; Chang, C.C. Electric dipoles on clean mica surfaces. Surf. Sci. 1969, 14, 39–51.

27. Pan, D.; Liu, L.M.; Tribello, G.A.; Slater, B.; Michaelides, A.; Wang, E. Surface energy and surface proton order of ice-Ih. Phys. Rev. Lett. 2008, 101, 155703.

28. Ding, P.; Li-Min, L.; Gareth, A.T.; Ben, S.; Angelos, M.; Enge, W. Surface energy and surface proton order of the ice Ih basal and prism surfaces. J. Phys. Condens. Matter 2010, 22, 074209.

29. Buch, V.; Groenzin, H.; Li, I.; Shultz, M.J.; Tosatti, E. Proton order in the ice crystal surface. Proc. Natl. Acad. Sci. USA 2008, 105, 5969–5974.

30. Hill, R.J. A further refinement of the barite structure. Can. Mineral. 1977, 15, 522–526.

31. Walker, D.; Verma, P.K.; Cranswick, L.M.D.; Jones, R.L.; Clark, S.M.; Buhre, S. Halite-sylvite thermoelasticity. Am. Mineral. 2004, 89, 204.

32. Drits, V.A.; Zviagina, B.B.; McCarty, D.K.; Salyn, A.L. Factors responsible for crystal-chemical variations in the solid solutions from illite to aluminoceladonite and from glauconite to celadonite. Am. Mineral. 2010, 95, 348–361.

33. Catti, M.; Ferraris, G.; Hull, S.; Pavese, A. Powder neutron diffraction study of 2M1 muscovite at room pressure and at 2 GPa. Eur. J. Mineral. 1994, 6, 171–178.

34. Goto, A.; Hondoh, T.; Mae, S. The electron density distribution in ice-Ih determined by single-crystal X-ray diffractometry. J. Chem. Phys. 1990, 93, 1412–1417.

35. Bernal, J.D.; Fowler, R.H. A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions. J. Chem. Phys. 1933, 1, 515–548.

36. Fenter, P.; McBride, M.T.; Srajer, G.; Sturchio, N.C.; Bosbach, D. Structure of barite (001)- and (210)- water interfaces. J. Phys. Chem. B 2001, 105, 8112–8119.

37. Davey, R.J.; Black, S.N.; Bromley, L.A.; Cottier, D.; Dobbs, B.; Rout, J.E. Molecular design based on recognition at inorganic surfaces. Nature 1991, 353, 549–550.

38. Black, S.N.; Bromley, L.A.; Cottier, D.; Davey, R.J.; Dobbs, B.; Rout, J.E. Interactions at the organic/inorganic interface - Binding motifs for phosphonates at the surface of barite crystals. J. Chem. Soc. Faraday Trans. 1991, 87, 3409–3414.

39. Hartman, P.; Strom, C.S. Structural morphology of crystals with the barite (BaSO4) structure: A revision and extension. J. Cryst. Growth 1989, 97, 502–512.

40. Murray, H.H. Structure and composition of the clay minerals and their physical and chemical properties. In Developments in Clay Science; Murray, H.H., Ed.; Elsevier: Amsterdam, The Netherlands, 2006; Volume 2, pp. 7–31.

41. Stixrude, L.; Peacor, D.R. First-principles study of illite–smectite and implications for clay mineral systems. Nature 2002, 420, 165.

42. Teich-McGoldrick, S.L.; Greathouse, J.A.; Cygan, R.T. Molecular dynamics simulations of uranyl adsorption and structure on the basal surface of muscovite. Mol. Simul. 2014, 40, 610–617.

43. Wang, L.; Liu, R.; Hu, Y.; Sun, W. pH effects on adsorption behavior and self-aggregation of dodecylamine at muscovite/aqueous interfaces. J. Mol. Graph. Model. 2016, 67, 62–68.

44. Odelius, M.; Bernasconi, M.; Parrinello, M. Two dimensional ice adsorbed on mica surface. Phys. Rev. Lett. 1997, 78, 2855–2858.

45. Purton, J.; Allan, N.; Blundy, J. Impurity cations in the bulk and the {001} surface of muscovite: An atomistic simulation study. J. Mater. Chem. 1997, 7, 1947–1951.

46. Kenneth, G.L. The physics of snow crystals. Rep. Prog. Phys. 2005, 68, 855.

47. Olijve, L.L.; Meister, K.; DeVries, A.L.; Duman, J.G.; Guo, S.; Bakker, H.J.; Voets, I.K. Blocking rapid ice crystal growth through nonbasal plane adsorption of antifreeze proteins. Proc. Natl. Acad. Sci. USA 2016, 113, 3740–3745.

48. Engel, E.A.; Monserrat, B.; Needs, R.J. Vibrational effects on surface energies and band gaps in hexagonal and cubic ice. J. Chem. Phys. 2016, 145, 044703.

49. Shoaib, M.A.; Choi, C.H. Adsorptions of HOCl on ice surface: Effects of long-range electrostatics, surface heterogeneity, and hydrogen disorders of ice crystal. J. Phys. Chem. C 2012, 116, 3694–3701.

50. Park, S.C.; Moon, E.S.; Kang, H. Some fundamental properties and reactions of ice surfaces at low temperatures. Phys. Chem. Chem. Phys. 2010, 12, 12000–12011.

51. Materer, N.; Starke, U.; Barbieri, A.; Van Hove, M.A.; Somorjai, G.A.; Kroes, G.J.; Minot, C. Molecular surface structure of ice (0001): Dynamical low-energy electron diffraction, total-energy calculations and molecular dynamics simulations. Surf. Sci. 1997, 381, 190–210.

52. Eisenberg, D.; Coulson, C.A. Energy of formation of D-defects in ice. Nature 1963, 199, 368–369.

53. de Koning, M.; Antonelli, A. On the trapping of Bjerrum defects in ice-Ih: The case of the molecular vacancy. J. Phys. Chem. B 2007, 111, 12537–12542.

54. Agmon, N. The Grotthuss mechanism. Chem. Phys. Lett. 1995, 244, 456–462.

55. Francis, G.P.; Payne, M.C. Finite basis set corrections to total energy pseudopotential calculations. J. Phys. Condens. Matter 1990, 2, 4395.

56. Liu, T.; Xue, L.; Guo, X.; Zheng, C.G. DFT study of mercury adsorption on α-Fe2O3 surface: Role of oxygen. Fuel 2014, 115, 179–185.

57. Perdew, J.P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865–3868.

58. Tkatchenko, A.; Scheffler, M. Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data. Phys. Rev. Lett. 2009, 102, 073005.

59. Gillan, M.J.; Alfè, D.; Michaelides, A. Perspective: How good is DFT for water? J. Chem. Phys. 2016, 144, 130901.

60. Parq, J.H.; Lee, S.K.; Lee, S.M.; Yu, J. Ab initio study of elastic properties of high-pressure polymorphs of CO2 phases II and V. J. Phys. Chem. C 2016, 120, 23152–23164.

61. Lee, B.H.; Lee, S.K. Effect of lattice topology on the adsorption of benzyl alcohol on kaolinite surfaces: Quantum chemical calculations of geometry optimization, binding energy, and NMR chemical shielding. Am. Mineral. 2009, 94, 1392–1404.

62. Barzilai, J.; Borwein, J.M. Two-point step size gradient methods. IMA J. Numer. Anal. 1988, 8, 141–148.

63. Schwarz, K.; Blaha, P.; Madsen, G.K.H. Electronic structure calculations of solids using the WIEN2k package for material sciences. Comput. Phys. Commun. 2002, 147, 71–76.

64. Hofer, T.S.; Randolf, B.R.; Rode, B.M. The hydration of the mercury(I)-dimer—A quantum mechanical charge field molecular dynamics study. Chem. Phys. 2008, 349, 210–218.

65. Rode, B.M.; Hofer, T.S.; Randolf, B.R.; Schwenk, C.F.; Xenides, D.; Vchirawongkwin, V. Ab initio quantum mechanical charge field (QMCF) molecular dynamics: A QM/MM—MD procedure for accurate simulations of ions and complexes. Theor. Chem. Acc. 2006, 115, 77–85.

66. Vanderbilt, D. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Phys. Rev. B 1990, 41, 7892–7895.

67. Crecelius, E.; Trefry, J.; McKinley, J.; Lasorsa, B.; Trocine, R. Study of Barite Solubility and the Release of Trace Components to the Marine Environment; US Department of the Interior (Minerals Management Service, Gulf of Mexico Region): New Orleans, LA, USA, 2007; Volume OC5 Study MMS 2007-061, pp. 1–176.

68. Denney, D. Fate of mercury in drilling-fluid barite in the marine environment. J. Pet. Tech. 2003, 55, 66–67.

69. Uddin, M.K. A review on the adsorption of heavy metals by clay minerals, with special focus on the past decade. Chem. Eng. J. 2017, 308, 438–462.

70. Churchman, G.J.; Gates, W.P.; Theng, B.K.G.; Yuan, G. Clays and clay minerals for pollution control. In Developments in Clay Science; Bergaya, F., Theng, B.K.G., Lagaly, G., Eds.; Elsevier: Amsterdam, The Netherlands, 2006; Volume 1, pp. 625–675.

71. Shetaya, W.H.; Huang, J.H.; Osterwalder, S.; Mestrot, A.; Bigalke, M.; Alewell, C. Sorption kinetics of isotopically labelled divalent mercury (196Hg2+) in soil. Chemosphere 2019, 221, 193–202.

72. Lee, S.S.; Nagy, K.L.; Park, C.Y.; Fenter, P. Enhanced uptake and modified distribution of mercury(II) by fulvic acid on the muscovite (001) surface. Environ. Sci. Technol. 2009, 43, 5295–5300.

73. Lee, S.S.; Nagy, K.L.; Park, C.; Fenter, P. Heavy metal sorption at the muscovite (001)–fulvic acid interface. Environ. Sci. Technol. 2011, 45, 9574–9581.

74. Tran, L.; Wu, P.; Zhu, Y.; Liu, S.; Zhu, N. Comparative study of Hg(II) adsorption by thiol- and hydroxyl-containing bifunctional montmorillonite and vermiculite. Appl. Surf. Sci. 2015, 356, 91–101.

75. Lumsdon, D.G.; Evans, L.J.; Bolton, K.A. The influence of pH and chloride on the retention of cadmium, lead, mercury, and zinc by soils. J. Soil Contam. 1995, 4, 137–150.

76. Biester, H.; Zimmer, H. Solubility and changes of mercury binding forms in contaminated soils after immobilization treatment. Environ. Sci. Technol. 1998, 32, 2755–2762.

77. Moon, E.S.; Kim, Y.S.; Shin, S.H.; Kang, H. Asymmetric transport efficiencies of positive and negative ion defects in amorphous ice. Phys. Rev. Lett. 2012, 108, 226103.

78. Lee, C.W.; Lee, P.R.; Kang, H. Protons at ice surfaces. Angew. Chem. Int. Ed. 2006, 45, 5529–5533.

79. Lee, D.H.; Kang, H.; Kang, H. Tunneling diffusion of excess protons in amorphous solid water at 10 and 80 K. J. Phys. Chem. C 2019, 123, 3657–3663.

80. Moon, E.S.; Kang, H.; Oba, Y.; Watanabe, N.; Kouchi, A. Direct evidence for ammonium ion formation in ice through ultraviolet-induced acid-base reaction of NH3 with H3O+. Astrophys. J. 2010, 713, 906–911.

81. Lee, C.W.; Kang, H. UV Photolysis of glycine on ice films: Implication for photosynthesis and photodestruction of amino acids in interstellar medium. Bull. Korean Chem. Soc. 2015, 36, 784–788.

82. Kang, H. Chemistry of ice surfaces: Elementary reaction steps on ice studied by reactive ion scattering. Acc. Chem. Res. 2005, 38, 893–900.

83. Shoaib, M.A.; Choi, C.H. Adsorptions of formic and acetic acids on ice surface: Surface binding configurations and a possibility of interfacial proton transfer. J. Phys. Chem. C 2013, 117, 4181–4188.

84. Shoaib, M.A.; Choi, C.H. Na+, F–, Br–, and Cl– adsorptions and penetrations on an ice surface. ACS Earth Space Chem. 2018, 2, 56–63.

85. Moon, E.S.; Yoon, J.; Kang, H. Energy barrier of proton transfer at ice surfaces. J. Chem. Phys. 2010, 133, 044709.

86. Tonigold, K.; Groß, A. Adsorption of small aromatic molecules on the (111) surfaces of noble metals: A density functional theory study with semiempirical corrections for dispersion effects. J. Chem. Phys. 2010, 132, 224701.

87. Fukazawa, H.; Ikeda, S.; Mae, S. Incoherent inelastic neutron scattering measurements on ice-XI; the proton-ordered phase of ice-Ih doped with KOH. Chem. Phys. Lett. 1998, 282, 215–218.

88. Watkins, M.; VandeVondele, J.; Slater, B. Point defects at the ice (0001) surface. Proc. Natl. Acad. Sci. USA 2010, 107, 12429.

89. Watkins, M.; Pan, D.; Wang, E.G.; Michaelides, A.; VandeVondele, J.; Slater, B. Large variation of vacancy formation energies in the surface of crystalline ice. Nat. Mater. 2011, 10, 794–798.

90. He, C.; Zhang, W.; Li, Y. The glass transition behaviors of low-density amorphous ice films with different thicknesses. J. Chem. Phys. 2010, 133, 204504.

91. Chen, J.Y.; Yoo, C.S. High density amorphous ice at room temperature. Proc. Natl. Acad. Sci. USA 2011, 108, 7685–7688.

92. Ojamae, L. Crystalline ice: Amorphous on the surface. Nat. Mater. 2011, 10, 725–726.

93. Refson, K.; Tulip, P.R.; Clark, S.J. Variational density-functional perturbation theory for dielectrics and lattice dynamics. Phys. Rev. B 2006, 73, 155114.