Microporous and Mesoporous Materials Mercury exchange in zeolites Na-A and Na-Y studied by classical molecular dynamics simulations and ion exchange experiments

Classical molecular dynamics simulations have been employed to study the exchange of Na + for Hg 2+ in zeolite Na-A, with a Si/Al ratio of 1, and zeolite Na-Y, with Si/Al ratios of 2 and 5, in dry and hydrated conditions within the temperature range 330 – 360 K, to understand factors underpinning the performance of zeolites for water decontamination. A classical forcefield based on DFT energies has been developed for the interaction between the Hg 2+ ions and the zeolite O atoms. In terms of water diffusion, zeolite Na-A shows the lowest calculated diffusivity, followed by zeolite Na-Y (Si/Al=2) and Na-Y (Si/Al=5), as a consequence of differing pore dimensions and extra-framework ion loadings. In the absence of speciation anions, the Hg 2+ ions are consistently adsorbed at the supercage windows in both the LTA and FAU framework types. The reduced pore size of zeolite A leads to an average hydration number per Hg 2+ ion of <1.0, whilst the wider pore of zeolite Y exerts less steric hindrance, and thus the Hg 2+ hydration number reaches values between 1.0 and 2.0 in zeolite Y. These observations might indicate that Hg 2+ ions are more strongly immobilized in zeolite A than in zeolite Y. Preliminary measurements of mercury removal using these zeolites, as synthesized from bauxite and kaolin, seem to support these findings.


Introduction
The contamination of rivers by mercury is a serious environmental and health issue, which predominantly affects less developed countries and is derived from the mining of precious metals such as gold [1][2][3]. Even if effective actions are adopted to eliminate the industrial and artisanal methods that produce the undesirable waste, the mercury already in the water bodies must be physically extracted. In this regard, ion-exchange matrices could be a viable method for the removal of the poisonous cations [4,5].
Zeolites are microporous alumino-silicate materials widely known for their outstanding ion exchange characteristics and environmental benignity [6]. Their efficacies as ion-exchange matrices are derived from the replacement of Si(IV) at tetrahedral positions by Al(III), which leaves a net negative charge that is balanced by extra-framework counter-ions commonly belonging to the alkali and alkaline earth groups [6]. Therefore, the lower the Si/Al ratio, the more counter-ions susceptibility of these ions: the formation of partially hydrolysed ions follows the order Hg 2+ < Cd 2+ < Zn 2+ , where the ions occupy positions near the super-cage apertures, obstructing the exchange [10].
More recent reports of zeolite-mediated Hg 2+ removal include a fundamental study of topology effects [5], where a comparison of zeolites Y (Si/Al = 2.5), beta (Si/Al = 17) and mordenite (Si/Al = 7.5) observed superior removal of Hg 2+ by zeolite Y in batch experiments (5 mg/L), although the adsorption capacity was highest for beta. Notably, the Si/Al ratios and certainly particle size distributions were significantly different between the samples (3,6, and 42 μm for Y, beta and mordenite, respectively) which may well have contributed to the conclusion that Hg 2+ adsorption was controlled by intraparticle diffusion in Y and beta, but film diffusion in mordenite [5]. Other studies have focussed on surface functionalization, where it has been observed that modifying Ca-clinoptilolite with cetylpyridinium bromide improved cation exchange by a factor of ∼4 and significantly improved the saturation adsorptive capacity [11], while surface modification with hexadecyltrimethylammonium bromide was found to improve the removal of both Hg 2+ and Pb 2+ by clinoptilolite due to a combination of both ion exchange and complexation processes [12]. Aside from modification with organics, the immobilization of MnO 2 nanoparticles on zeolite NaP was shown to improve its adsorption capacity by at least a factor of 5 [13], where the use of the less energy demanding sonochemical method for their synthesis was of particular interest. The potential of composite materials has also been considered, such as muscovite/zeolite-phillipsite composites [14], which showed the highest decontamination for Hg 2+ at a pH of 4, and a Dubinin-Radushkevich adsorption energy for Hg 2+ lower than As 5+ and U 4+ by ∼25%. A magnetite/mordenite composite [15] was also studied, and while the magnetite component did not add to the material's adsorbent capacity, it allowed for its simple magnetic removal from the aqueous medium.
Considering the importance of the removal of mercury from water, and the clear advantages and capabilities that zeolites present for this application, we have employed classical molecular dynamics simulations to study the exchange of Na + for Hg 2+ in zeolite Na-A (of the LTA framework topology, Si/Al = 1) and zeolite Na-Y (of the FAU framework topology) with two different Si/Al ratios of 2 and 5. The simulations, although performed ignoring the speciation effects of competing anions, such as Cl − , provide insights for the rationalization of macroscopic ion exchange experiments based on adsorption. We have analysed in detail the sites where the Na + and Hg 2+ ions adsorb, describing their interaction with water according to their location in the zeolite framework. The simulations and their comparison to experiment are explained based on the dimension of the pore system, the concentration of extra-framework ions and the mobility of water through the zeolite framework. We also present a preliminary study on the same zeolite topologies, synthesized from abundant and cheap feedstock materials, which shows a variable performance in the removal of Hg 2+ depending on the zeolite topology.

Computational
The codes General Utility Lattice Program (GULP) [18,19] and DL_POLY 4 [20] are used for geometry optimizations and molecular dynamics simulations using classical potentials, respectively. The structure visualizations shown in this work are produced with the software Visualization for Electronic and Structural Analysis (VESTA 3) [21].

Details of MD simulations
The classical description of the zeolites is based on the Born model of ionic solids [22], which employs Coulombic interactions and Buckingham potentials to define the forces acting between pairs of ions. The Coulombic contributions are calculated using the Ewald summation method [23,24], with the short range repulsions and dispersion forces described through Buckingham and Lennard-Jones potentials [25,26] [27]. The interatomic pairwise interaction for (Al 3+ , O 2− ) is included following Catlow et al. after an optimal fitting to the lattice properties of -Al 2 O 3 [28]. In addition to Coulombic forces, the interaction of the Na + ions with the zeolite framework is incorporated via a Buckingham potential to represent the pair (Na + , O 2− ) [29].
The flexible TIP3P model is employed to control the intra-and intermolecular interactions in water, where the flexibility is represented by the adoption of harmonic potentials for the O-H bond and the H-O-H angle, and a Lennard-Jones parametrization to treat the interactions between the O atoms of different molecules [30].
The interaction between the water molecules and the zeolite framework is described by the Lennard-Jones potentials reported by Jaramillo et al. for the pairs (O 2− , O * ) and (O 2− , H * ), where the asterisk denotes water atoms [31]. The interactions (Si 4+ , O * ) and (Al 3+ , O * ) are represented by Buckingham potentials adapted from the values proposed by Sanders et al. for zeolites, following a similar re-scaling procedure for the parameter as reported by Schröder and collaborators [32]. The pair (Na + , O * ) is described with the same Buckingham potential used to represent the interaction (Na + , O 2− ) [29]. The full set of interatomic potentials is listed in section S1, Table S1 of the Supplementary Material; the performance of this set is compared to previous computational work reported by Demontis and collaborators (Supplementary Material, section S2) [33].
The MD simulations, at temperatures of 300, 330 and 360 K, consist of an initial equilibration in the NVE ensemble for 2 ns, followed by a NVT ensemble run for 3 ns, using a Berendsen thermostat with a time constant of 1 ps for thermal energy exchange [34]. Following the equilibration runs, production runs of 12 ns in the NVE ensemble are carried out. A time step of 0.5 fs is used in the simulations, saving the atomic coordinates every 2000 steps. The difference between consecutive stored steps ( ) is thus 1 ps.

Buckingham potential for (Hg 2+ , O)
We have used DFT calculations, and the (Na + ,O 2− ) parameters already available, to parameterize the Buckingham potential needed to describe the (Hg 2+ ,O) interactions. The DFT calculations are performed with the planewave code Vienna Ab-initio Simulation Package (VASP) [35][36][37][38], while the interatomic potential code General Utility Lattice Programme (GULP) [18,19] is used for the constant volume geometry optimizations.
To perform the DFT calculations, we use the formulation of the generalized gradient approximation (GGA) by Perdew, Burke, and Ernzerhof (PBE) to account for the exchange and correlation contributions to the system electronic energy [39], treating the valence electrons with a basis set of planewaves with a kinetic energy cutoff of 400 eV whilst their nodal features and the inner electronic levels of the atom are described by the projector-augmented-wave (PAW) method [40,41]. The DFT functional is coupled with Grimme's atomic pairwise method using the Becke-Johnson damping function to include the long-range dispersion forces [42][43][44]. The first Brillouin zone is sampled with a k-points mesh of 4 × 4 × 4, which is generated with the Monkhorst-Pack scheme [45]. We employ the Gaussian smearing method with a band width of 0.05 eV to evaluate the occupation of the electronic states and the integration over the k-points during the self-consistent procedure [46,47]. The thresholds for the convergence of the electronic energy and ionic forces are set at 10 −5 eV and 0.03 eV/Å, respectively. The calculations are carried out under periodic boundary conditions, allowing to implicitly include the three-dimensional extension of the solid using only the small set of atoms contained in the unit cell of the crystalline structure. The DFT-optimized cubic unit cell of the all-silica LTA framework, with a lattice parameter of 11.828 Å, is used for both the DFT and IP calculations, with a double Al substitution, and loaded with two Na + ions or one Hg 2+ ion. Following convention, we assume that = 0, and we search for the Buckingham potential parameters and for the interaction between the Hg 2+ ion and the O atoms of the zeolite framework by constructing increasingly finer grids in the ( , ) space. Then, we calculate the Na + /Hg 2+ exchange energies over these ( , ) grids, which are compared to a DFT-calculated reference value. The difference between the classical and DFT-calculated energies is minimized for = 0.350 Å and = 984.5 eV. In the present work, these values of and are used to describe the (Hg 2+ , O) interaction with both the zeolite and water O atoms. The values of and over the finest grid, and further details concerning the models are provided in section S3 of the Supplementary Material.

Processing of the MD trajectories
The method of multiple time origins ( 0 ) is used to improve the statistics of the data. For the mean-square displacement (MSD), the full trajectory of 12 ns is averaged into 10 ns by shifting 0 every 5 ps. A plot with satisfactory linearity is obtained, which is used to derive the self-diffusion coefficients from the Einstein relationship: The contact correlation function ( ) between water and the extraframework ions is calculated with the equation: where is the total number of water molecules. The function ( ) equals 1 if the molecule is in the coordination shell of the ion at time , otherwise it equals 0. The value of ( ) is then averaged over all ions of the same type. The function ( ) thus calculated describes the probability of continuously observing a water molecule within the coordination shell of an ion from time 0 to time .
The residence time is calculated according to the equation: If ( ) shows an exponential decay, then ( ) = − ∕ , and equals the relaxation time [48].
We have used the maximum of the first peak in the RDF profiles of the pairs (Na + ,O * ) and (Hg 2+ ,O * ) to define the cut-off for the coordination shell of each ion. We find that this cut-off depends more on the zeolite framework rather than the nature of the counterion, with LTA having a value of 2.7 Å and FAU of 2.9 Å for both Na + and Hg 2+ ions.
To obtain ( ), the multiple minimum method is applied, averaging the last 3 ns of simulations into 40 ps, shifting t 0 every 120 ps.

Structural models for MD simulations
The framework type of zeolite A is LTA (space group 3 ) whereas zeolite Y belongs to the FAU topology (space group 3 ). Both framework types share the same Composite Building Unit (CBU) consisting of 24 tetrahedral sites (T-sites) and denominated as the 'sodalite cage' (sod), shown in Fig. 1. In LTA, the framework can be built by connecting the sodalite units via double 4-membered rings (d4r), with a sodalite cage at each corner of an octahedron. In FAU, the sodalite units connect via hexagonal prisms or double 6-membered rings (d6r), placing the sodalite cages at the corners of a tetrahedron. The classical MD simulations are performed using the aluminiumsubstituted structures of these framework types, with a 2 × 2 × 2 expansion of the unit cell of zeolites A and Y, thus obtaining supercells with a dimension of approximately 50 Å.
We employ a Si/Al ratio of 1 for zeolite A, and Si/Al ratios of 2 and 5 for zeolite Y [49]. For the sake of simplicity we refer to zeolite Y with Si/Al ratios of 2 and 5 as ''Y(2)'' and ''Y(5)'', respectively. The location of Al atoms in the zeolite A framework is straightforward, with an alternating sequence of Al and Si atoms -Al-Si-Al-Si-(leaving out the connecting O atoms for clarity), in accordance with Lowenstein's rule. However, the number of possible configurations for zeolite Y with Si/Al ratios of 2 and 5 is much higher. Therefore, we have used the arrangements shown in Fig. 2, which maximize the distance between Al atoms in accordance with Dempsey's rule, with zeolite Y(2) having -Al-Si-Al-Si-and -Al-Si 2 -Al-Si 2 -sequences between nearest and next-nearest Al atoms, while zeolite Y(5) displaying -Al-Si 2 -Al-Si 2and -Al-Si 3 -Al-Si 3 -sequences.
The Na + ions are added to the aluminium-substituted zeolites to balance the framework charge. In the case of zeolite A, we saturate first all S1 positions with Na + and continue adding ions at S2 positions until a total number of 96 Na + per unit cell is reached (see Fig. 1 for site locations). In zeolite Y(2), out of 64 Na + per unit cell, 32 ions are placed at SII positions and the other 32 at SJ positions. Two initial configurations are used for zeolite Y(5); in the first one, the 32 Na + ions per unit cell are located at SII positions, while in the second structure, the ions are placed at SJ positions; thus, the labelling for zeolite Y with a Si/Al ratio of 5 is further split into ''Y(5,SII)'' and ''Y(5,SJ)''. The choice of the location for Na + is based on the most probable positions to be occupied in zeolites Na-A, Na-Y(2) and Na-Y(5) according to previous reports [16,17,50], and to allow a consistent comparison between the two Si/Al ratios used for zeolite Y regarding the effects of placing Na + at different initial positions.
Since in the present work we are particularly interested in the analysis of the Na + /Hg 2+ exchange under both dehydrated and hydrated conditions, but without the need to determine the exact amount of water in the system, we have selected a loading of 120 water molecules per unit cell to study the effect of hydration on the ion exchange, which allows reliable statistics and a clear contrast when comparing to the dehydrated framework. In addition, this water loading enables the validation of our set of interatomic potentials against the computational work by Demontis et al. who also used 120 water molecules per unit cell as an intermediate loading, within a range of 20 to 250 molecules, to study the diffusion of water in zeolite Na-Y [33].
For the study of the exchange of Na + for Hg 2+ , 16 Hg 2+ ions per unit cell are randomly loaded into the systems, while retaining the original concentration of Na + ions acting as counter-ions; this causes a net positive charge of +256 e − , which is neutralized by a background negative charge uniformly distributed through the entire simulation cell and implicitly created when the Ewald method is applied. This background charge ensures a proper convergence of the Ewald sum. Table 1 Diffusion coefficients of water in zeolites Na-A, Na-Y(2), Na-Y(5,SII) and Na-Y(5,SJ) in the absence and presence of Hg 2+ at the temperatures 300, 330 and 360 K (10 −10 m 2 s −1 units), and the corresponding activation energy for translational diffusion (kJ/mol units).
In the absence of Hg 2+ However, it can also lead to an artificial net pressure and energy artefacts that particularly affect systems where the volume of the cell and the charge of the system are allowed to vary [51,52]. DL_POLY amends these anomalies by introducing a volume correction term proposed by Fuchs, who originally aimed to explain why noble metals are facecentred, while alkali metals are body-centred [53]. We ensure that the concentration of Hg 2+ does not affect the equilibration and stability of the system, allowing to observe an explicit exchange of Na + for Hg 2+ . For the sake of simplicity, these calculations neglect the presence of speciation anions, such as Cl − , even though they may have an impact on the adsorption of Hg 2+ in the zeolite. However, the main goal of the present work -our first in this study -comprises the analysis of the Na + /Hg 2+ exchange under ideal conditions, considering only the influence of water. Future research should also include Cl − anions, which would allow a better direct comparison to experiment. We should note here that performing the simulations under dehydrated and hydrated conditions provides initial insight into how competitive interactions, such as the one occurring between Hg 2+ and water, affect the binding of these ions to the zeolite surface, thereby providing an understanding of the Na + /Hg 2+ exchange at the atomic level.

Experimental
Zeolites Na-A, Na-Y(2.0) and Na-Y(3.1) were synthesized following a previously published method using bauxite and kaolin as cheaply available feedstock materials [54,55]. Further details concerning the synthesis and structural characterization of the zeolites can be found in the Supplementary Material, section S4.
The removal of Hg 2+ ions was studied by preparing a stock mercury solution (1000 mg/l) dissolving 2.21 g of HgCl 2 in deionized water. The model mercury-polluted water used for the experiment was prepared by spiking deionized water to the required concentration using the stock mercury solution. The adsorption experiments were performed using 200 mL mercury solutions with initial concentrations of 3.00 and 100.00 mg/L for 1440 min, at a pH between 7.0 and 7.3. The Hg 2+ removal was initiated by adding 5 g of the zeolite and stirring continuously at 90 revolutions per minute. A blank mercury adsorption experiment, without adsorbent, was carried out for comparison. The adsorption experiments are carried out in triplicate. Each filtered sample is analysed by inductively coupled plasma-atomic emission spectrometer (ICP-AES). Fig. 3 shows the MSD of water for the four studied zeolite systems in the absence and presence of Hg 2+ ions, with the resulting diffusion coefficients listed in Table 1.  In the absence of Hg 2+ , the slowest water diffusion is observed for zeolite Na-A, followed by zeolite Na-Y(2), with the fastest diffusion occurring in zeolite Na-Y(5), as expected since zeolite Na-A has both the highest concentration of Na + ions and the narrowest windows between larger cages, 8MR (∼ 4 Å) compared to 12MR (∼ 7.5 Å) in framework type FAU. The diffusion coefficients for zeolite Na-A show a minimum of 0.34 × 10 −10 m 2 s −1 at 300 K, increasing to 1.88 × 10 −10 m 2 s −1 when the temperature is raised to 360 K. The difference in diffusion between zeolites Na-Y(2) and Na-Y(5) is a consequence of the variation in concentration of Na + ions between the two frameworks, i.e. 64 ions per unit cell for Na-Y(2) compared to 32 ions per unit cell for Na-Y(5). A higher concentration of Na + will inevitably decrease the mobility of water molecules throughout the micropore system as a result of the interactions between water and the extra-framework ions. At the highest studied temperature of 360 K, the diffusion coefficient of water in Na-Y(2) is a factor of 2.4 lower than in Na-Y(5), with values of 2.64, 6.21 and 6.28 × 10 −10 m 2 s −1 for Na-Y(2), Na-Y(5,SII) and Na-Y(5,SJ), respectively. We can observe that the initial distribution of Na + ions in zeolites Na-Y(5,SII) and Na-Y(5,SJ) does not further affect the diffusivities.

Water diffusion
When Hg 2+ is present in the system, we observe that the diffusion coefficients decrease by at least a factor of 3 in zeolite Na-A, with a more subtle reduction in zeolites Na-Y(2) and Na-Y(5). This observation shows that water diffusion is more sensitive to an increase in the number of extra-framework ions in the already constricted micropores of zeolite Na-A compared to zeolite Na-Y.
In the absence of Hg 2+ , the activation energy of diffusion varies within the narrow range of 19 -26 kJ/mol for the four zeolite systems, highlighting that the interaction with the extra-framework ions is the dominant factor controlling the diffusion of water. However, pore size and counter-ion loading still play a role in the magnitude of the energy barrier, as illustrated by the ranking: In the presence of Hg 2+ , a marginal, but consistent increase in activation energy is observed, although by less than 2 kJ/mol. Fig. 4 provides information regarding the occupancy of each framework site by Na + and Hg 2+ ions, averaged over the entire production run.

Site occupancy
If the simulation is carried out without water and in the absence of Hg 2+ ions, the Na + ions retain the initial positions set at the beginning of the simulation, with zeolite Y(5,SJ) as the only exception, where a partial migration is observed from SJ sites to the vacant SI' and SII sites. Thus, under dry conditions, the extra-framework ions tend to adsorb in sites where they can maximize the interaction with framework O atoms.
Upon Hg 2+ adsorption, but under the same dry conditions, a stronger redistribution of the Na + ions is observed. In the case of zeolite Na-A, the site preference by the Hg 2+ ions follows the order S3 < S2 < S1 with ratios of 1.00:1.22:1.66. Consequently, between 10 and 20% of Na + ions migrate from S1 and S2 sites and re-adsorb preferentially at the S3 site. In zeolite Na-Y(2), Hg 2+ favours SII and SJ sites, adsorbing approximately 40 and 47% of a total of 128 Hg 2+ ions in the simulation cell, respectively. Additionally, just under 20 Hg 2+ ions are able to diffuse inside the sodalite cages to adsorb at SI' sites; here the Na + /Hg 2+ exchange removes Na + ions from SII and SJ sites, with most Na + ions diffusing into the sodalite cages, ending up at SI'. Zeolite Na-Y(5,SII), which has more SJ sites readily available, shows a different profile for the Hg 2+ ions, with 60% of Hg 2+ adsorbing at SJ, and the remainder at SI' and SII. Zeolite Na-Y(5,SII) also shows the strongest Na + re-arrangement upon addition of Hg 2+ (when compared to simulations performed in the absence of Hg 2+ ). Under dehydrated conditions and without Hg 2+ , Na + ions do not move significantly from their initial positions at SII sites over the entire simulation. In contrast, in the presence of Hg 2+ , approximately half of Na + ions leave the SII sites and re-adsorb at the SI' and SJ sites. In zeolite Na-Y(5,SJ), Hg 2+ ions are almost equally distributed between SII and SJ sites, in a similar fashion to zeolite Na-Y(2), while Na + ions show a distribution that resembles that of zeolite Na-Y(5,SII), with a higher adsorption at SII sites followed by a fairly similar occupancy of SI' and SJ sites. These results indicate that, in the absence of water, Hg 2+ ions favour relatively equally the S1 and S2 sites in zeolite Na-A, and the SII and  Fig. 1; when the ion cannot be considered adsorbed at any framework site, it is classified as ''Free'', which is represented by ''F'' in the horizontal axis of the graph). The counting is performed for each ion under different simulation conditions: (first column) Na + in dehydrated conditions in the absence of Hg 2+ , (second column) Na + in dehydrated conditions in the presence of Hg 2+ , (third column) Hg 2+ in dehydrated conditions, (fourth column) Na + in hydrated conditions (120 water molecules per unit cell) in the absence of Hg 2+ , (fifth column) Na + in hydrated conditions in the presence of Hg 2+ , (sixth column) Hg 2+ in hydrated conditions. The temperatures 300, 330 and 360 K are represented by blue, green and red colours, respectively, with the initial configuration of Na + ions indicated in grey. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) SJ sites in zeolites Na-Y(2) and Na-Y(5). Zeolite Na-Y(5,SII) is the only exception owing to an initial Na + configuration that leaves a large number of available SJ sites ready for Hg 2+ adsorption.
The main consequence of including water in the microporous systems is the unambiguous preference by the Hg 2+ ions to adsorb at the supercage windows, i.e. S2 sites in zeolite Na-A, and SJ sites in zeolite Na-Y.
The profiles for the distribution of Na + ions among the different sites in the presence of water follow similar trends to those observed in the simulations performed under dry conditions. In particular, the data for zeolite Na-Y(2) in the absence of Hg 2+ can be compared to the site occupancies reported by Hunger et al. for zeolite Na-X (Si/Al = 1.2) using neutron powder diffraction [17]. For a loading of 72 water molecules per unit cell, occupancies of 28%, 34% (fixed for 72 water molecules per unit cell, but it remains between 34 and 37% for other water loadings when allowed to change during refinement, see Table 3 in Ref. [17]) and 33% for the SI', SII and SJ sites, respectively, were obtained. Our corresponding Na + ion distribution in zeolite Na-Y(2) at 360 K is 22, 31 and 47% for SI', SII and SJ sites, respectively, which is certainly comparable to experiment considering that Hunger and coworkers used a zeolite with a Si/Al ratio of 1.2, i.e. approximately half the value used in our simulations. We should also note that these authors observed the formation of a new phase with a symmetry different from 3 when loading 120 molecules per unit cell [17]. In general, we do not observe any sustained correlation in the occupancy of the different framework sites with temperature. There are a few exceptions, such as the adsorption of Na + ions at SII and SJ sites in hydrated conditions and in the absence of Hg 2+ ions for zeolite Na-Y(5,SII). Tables S5 and S6 of the Supplementary Material list the full set of occupancies plotted in Fig. 4.

Water distribution per adsorption site
The features of the hydration shell around each ion provide relevant information in terms of the adsorption properties of the zeolite. Therefore, we have analysed in more detail the interaction of water with the Na + and Hg 2+ ions. Fig. 5 shows the average number of water molecules coordinating to the Na + and Hg 2+ ions, splitting the number according to the sites where the ions are adsorbed. Water favours coordination to ions located at the supercage windows, which is particularly obvious for Na + in zeolite Na-A. Here, 60 to 70% of Na + ions are adsorbed at S1 sites (centre of the 6MR), whereas only 10 to 20% are located at S2 sites (entry of the supercage). However, less than 6% of the total number of water molecules coordinating Na + do so when the ions are at the centre of the 6MR, with 60 to 70% of the water binding Na + ions at the supercage windows. This behaviour shows that the higher the coordination of the ions by framework O atoms, the lower their hydration tends to be.
We observe a consistent increase in the number of water molecules coordinating the Hg 2+ ions when the loading of Na + ions decreases, which may have implications for the retention of Hg 2+ ions by the zeolite framework. In zeolite Na-A, with the highest extra-framework loading, 60 to 70 molecules directly bind to the Hg 2+ ions, leading to an average hydration number of approximately 0.5. This value increases by a factor of 2 -3 in zeolite Na-Y(2), with a total of approximately 170 molecules involved in the coordination of Hg 2+ , raising the average hydration number to 1.3. In zeolites Na-Y(5,SII) and Na-Y(5,SJ), the total number of water molecules coordinating to the Hg 2+ ions is actually higher than for Na + , with values that vary between 220 and 260 for Hg 2+ , compared to 180 to 200 molecules for Na + , yielding hydration number ranges of 1.7 -2.0 and 0.7 -0.8 for Hg 2+ and Na + , respectively. At least 75% of this Hg-coordinated water is around SJ sites in zeolites Na-Y(2) and Na-Y(5). The full set of data plotted in Fig. 5 can be found in Table S7 of the Supplementary Material.
The comparison between zeolites Na-A and Na-Y(2) emphasizes the role of the pore diameter in the hydration of the extra-framework ions. With a Si/Al ratio of 1, zeolite A can load 1.5 times more extraframework ions than zeolite Y(2). Nevertheless, in the absence of Hg 2+ , the total number of water molecules hydrating the Na + ions in zeolite Na-A is 240, against 370 molecules for zeolite Na-Y(2) (although we set cutoff radii of 2.7 and 2.9 Å for zeolites Na-A and Na-Y(2), respectively, the number is still lower for zeolite Na-A if a radius of 2.9 Å is assumed for both zeolites). These numbers decrease proportionally when Hg 2+ is added to the system, down to 180 and 260 -280, respectively. This trend is visually supported by the analysis of the radial distribution function (RDF) of water O atoms around the Na + ions. Fig. 6 shows the RDF of the pair (Na + , O * ) in the absence of Hg 2+ , for zeolites Na-A and Na-Y(2), at 360 K. The RDF profile of zeolite Na-Y(2) displays a single peak with a maximum just under 3.0 Å that flattens after 3.5 Å. In contrast, zeolite Na-A shows two consecutive peaks at 2.7 and 3.7 Å, both smaller in intensity compared to the single peak of zeolite Na-Y(2). If we consider that a large proportion of the total number of water molecules are coordinating ions adsorbed at the supercage entry (60 to 70% in zeolite Na-A, as stated earlier, and almost 80% in zeolite Na-Y(2)), then the RDF analysis is additional evidence that ions at S2 sites in zeolite Na-A cannot directly coordinate as many water molecules as the ions at SJ sites in zeolite Na-Y(2). We can ascribe this outcome to the appreciably smaller diameter of the supercage window in framework type LTA compared to type FAU (8MR vs. 12MR).
These trends could provide insight into what we could expect to occur if speciation anions are also present in the simulation system. For example, although the Hg 2+ ⋯Cl − interaction is stronger, leading to the formation of chloride complexes, smaller supercage windows and higher concentrations of Na + ions may be more likely to create the conditions conducive to the cleavage of the Hg 2+ ⋯Cl − coordination, resulting in a better immobilization of Hg 2+ .
We have calculated the residence time of the water molecules within the coordination shell of the extra-framework ions located at the supercage windows of zeolites A and Y(2), with the values listed in Table 2. In the absence of Hg 2+ , the stronger confinement in zeolite Na-A forces the water molecules to remain coordinated to the Na + ions for an average time that ranges between 2.3 and 2.6 ps, depending on the temperature. This time is reduced to 1.4 -1.5 ps in zeolite Y (2). When Hg 2+ is present in the pore system, the residence time C. Hernandez-Tamargo et al.   Table 3 lists the efficiency of zeolites Na-A and Na-Y towards the removal of Hg 2+ from the solution. Our zeolite Na-A is able to extract 98.6 and 96.4% of Hg 2+ from solutions with initial concentrations of 3.0 and 100.0 mg/L, respectively. Similar levels of extraction have been reported for nanostructured zeolite Na-A [56]. However, our zeolites Na-Y(2.0) and Na-Y(3.1) perform poorly, removing only 37.7 and 11.5% of Hg 2+ from the 3.0 mg/L solution, respectively, and between 2.5 and 2.6% from the 100.0 mg/L solution. Zeolites Na-Y(2.0) and Na-Y(3.1) show a low tolerance for Hg 2+ uptake, being rapidly saturated when the concentration of Hg 2+ is raised, which is in striking contrast to previous experimental reports [5,10]. For example, Murthy and collaborators optimized several parameters for the removal of mercury (pH, zeolite dosage, contact time), and were able to extract 98% of Hg 2+ from solution by employing zeolite Na-Y (see Table 3) [5]. We can tentatively rationalize our unsatisfactory results by considering the stability of different Hg 2+ ⋯Z − aggregates. Murthy and collaborators used a stock solution of HgSO 4 [5], while we employed HgCl 2 to study the removal of Hg 2+ from water. If we compare the dissociation constant of the aggregates HgSO 4 (aq) and HgCl 2 (aq) in solution, we note that it is much more probable to find free Hg 2+ ions in the presence of SO 2− 4 (aq) compared to Cl − (aq) [57]. Therefore, a lower proportion of free Hg 2+ ions in the stock solution could explain why our zeolites Na-Y(2.0) and Na-Y(3.1) show a lower exchange efficiency when compared to the reports in the literature. Nevertheless, even with a stock solution of HgCl 2 , our zeolite Na-A is able to extract most of the Hg 2+ ions from water. As discussed above, a reasonable explanation for the efficiency of Na-A is the observed effect of the zeolite topology and Na + concentration on the hydration number of the Hg 2+ ions, i.e. the smaller supercage windows and the higher concentration of Na + ions in Na-A are more effective in disrupting the Hg 2+ ⋯Cl − interaction when these aggregates are adsorbed in the zeolite, thereby freeing the Hg 2+ ions and increasing the chances of their immobilization in the exchange matrix. In contrast, zeolite Na-Y, with larger supercage windows and lower Na + concentration, should be able to adsorb the HgCl 2 aggregates without causing a similar disruption in the Hg 2+ ⋯Cl − interaction, thus hindering the retention of Hg 2+ . We may also consider that the density of negative charge is higher in the framework of zeolite A, and, thus, a stronger electrostatic repulsion of the Cl − ions should be expected compared to zeolites Y(2) and Y(5). This repulsion, together with a larger concentration of extra-framework ions and smaller supercage windows, should increase the propensity towards Hg 2+ removal from the solution. However, we could only achieve a definitive analysis of this hypothesis by explicitly including the Cl − ions in the simulation. Yet, the initial trends identified here can already help to rationalize the variable efficiency with which different zeolites remove Hg 2+ from water.

Conclusions
In the present work we have combined classical molecular dynamics simulations and ion-exchange experiments to study the efficiency of zeolites Na-A and Na-Y to remove Hg 2+ from water. To enable the simulations, a new forcefield to describe the interaction between the Hg 2+ ions and the zeolite O atoms was derived based on DFT calculations. The computer simulations allowed a comprehensive analysis at the atomic scale of the Na + /Hg 2+ exchange, based on the adsorption sites of the extra-framework ions, the Si/Al ratio, the hydration conditions and the topology of the zeolites, while our experiments provided preliminary results regarding the efficiency of zeolites Na-A and Na-Y, with Si/Al ratios between 1.0 and 3.1, to extract Hg 2+ from batch solutions of HgCl 2 .
The simulations have shown that the Hg 2+ ions tend to adsorb at the supercage windows of both LTA and FAU framework types in the presence of water, with their hydration number increasing with the size of the pore diameter and the decrease in the number of extraframework Na + ions. We observed an average of 0.5 water molecules coordinating the Hg 2+ ions in zeolite Na-A, with this number increasing to 1.3 and at least 1.7 in zeolites Na-Y(2) and Na-Y(5), respectively. The analysis of the radial distribution function of water coordinating the extra-framework ions shows that the much smaller 8MR windows in zeolite A hinder the hydration of the ions located at S2 sites (supercage windows) compared to the equivalent SJ sites in zeolite Y.
We did not consider speciation ions, such as Cl − , in the simulations, aiming to simplify the models and the analysis of the data. Nevertheless, if the observed effect of the zeolite topology and the number of extra-framework Na + ions on the hydration of Hg 2+ is tentatively applied to the Hg 2+ ⋯Cl − interaction, we can conclude that zeolite A should be more effective than zeolite Y for the removal of Hg 2+ from solutions with strong speciation of the Hg 2+ ions. This line of reasoning was supported by preliminary experiments, where our zeolite Na-A was able to extract almost all the Hg 2+ from batch solutions of HgCl 2 . However, zeolite Na-Y, with larger supercage windows and lower concentration of Na + , performed very poorly under the same experimental conditions. Simulations where Cl − ions are explicitly included are yet to be performed, but the initial results presented in this work are promising, and already suggest important implications for the design of effective zeolite materials for the removal of heavy metals from water.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.