Arsenic Adsorption onto Minerals: Connecting Experimental Observations with Density Functional Theory Calculations

Abstract

A review of the literature about calculating the adsorption properties of arsenic onto mineral models using density functional theory (DFT) is presented. Furthermore, this work presents DFT results that show the effect of model charge, hydration, oxidation state, and DFT method on the structures and adsorption energies for As^III and As^V onto Fe³⁺-(oxyhydr)oxide cluster models. Calculated interatomic distances from periodic planewave and cluster-model DFT are compared with experimental data for As^III and As^V adsorbed to Fe³⁺-(oxyhydr)oxide models. In addition, reaction rates for the adsorption of As^V on α-FeOOH (goethite) (010) and Fe³⁺ (oxyhydr)oxide cluster models were calculated using planewave and cluster-model DFT methods.

1. Introduction

1.1. Arsenic Chemistry, Geochemistry, Prevalence, and Toxicity

The study of arsenic (As) adsorption on mineral surfaces is necessary to understand both the distribution and mobility of As species in nature as well as to develop remediation strategies for As waste sites. Arsenic is found in a variety of geochemical environments at aqueous concentrations varying from <0.5 to >5000 μg/L, and is found in a variety of geochemical environments [1,2]. Natural and anthropogenically-mediated biogeochemical interactions among arsenic species, biota, and minerals can affect the distribution, mobility, and toxicity of As in the environment [2,3,4,5,6,7,8]. Although recent work has posited that arsenic could be a potential biochemical and astrobiological proxy for phosphorus during biological evolution [9], this hypothesis is controversial [10].

Arsenic can occur in both inorganic (iAs) and organic (oAs) forms; the chemical form of As, or species, as well as the concentration of As, affects the solubility, mobility, reactivity, bioavailability, and toxicity of As [11,12]. iAs occurs predominately as either arsenious acid (H_nAs^IIIO₃ⁿ⁻³; sometimes called arsenous acid) in reducing environments, or arsenic acid (H_nAs^VO₄ⁿ⁻³) in oxidizing environments, where n = 0, 1, 2, or 3 [12,13,14]. Both the oxidation and protonation states of iAs depend on the physiochemical conditions of the sample environment [2]. For example, the three pK_a values for arsenic acid are 2.2., 6.9, and 11.5 [15]; therefore, the pH of the environment will affect the protonation state of H₃AsO₄, which will, in turn, affect the mobility, reactivity, and bioavailability of iAs^V. In addition to iAs, methylated As^III and As^V compounds such as monoarsinic acid (MMA^V) and dimethylarsinic acid (DMA^V) occur both naturally and due to anthropogenic sources [12,13]. Arsenic toxicity depends on the species present; a general trend of decreasing toxicity is: R₃As > H₃AsO₃ > H₃AsO₄ > R₄As⁺ > As⁰, where R is an alkyl group or a proton [12,16].

Arsenic originating from natural water-rock interactions of surface and groundwater [2,6] and from anthropogenic sources such as acid mine drainages [17,18] can lead to drinking water contamination. Biogeochemical processes result in organic As compounds accumulating in oil, shale, and coal [19]. DMA and MMA have been used as herbicides, pesticides, and defoliants and pose a contamination problem to surface and groundwater [20]. Roxarsone (C₆AsNH₆O₆) is used as a supplement in chicken feed and ends up in waste from poultry operations [21]. In addition to natural and anthropogenic groundwater contamination by As, anthropogenic As is contaminating the oceans, and subsequently seafood [16], which could have human health implications.

Human diseases caused by As contamination include various cancers, liver disease, cardiovascular disease, and an increase in mortality from pulmonary tuberculosis [22,23,24,25,26,27,28]. Groundwater that contains As concentrations >10 μg/L limit set by the World Health Organization (WHO) places more than ten million people at risk from arsenicosis. Locales ranging from Southern Bangladesh, India, Argentina, Chile, and Vietnam have groundwater sources with As concentrations >10 μg/L [2,3,6,25,26,28,29,30,31,32,33].

The remediation of arsenic in aquifers can be challenging due to the size of contaminated groundwater systems and varying biogeochemical conditions. For instance, the shallow groundwater of the Chaco Pampean Plain of Argentina spans 10⁶ km² and contains 10 to 5300 μg As/L [34]. A variety of As remediation methods were reviewed recently [35,36]; these methods have shown varying success rates. For example, an in situ study of an alkaline aquifer showed relatively low adsorption of iAs [37]. Conversely, a permeable reactive barrier study showed As removal to <5 μg/L due to induced sulfate reduction and the presence of zero-valent Fe [38], and experiments with household sand filters in Vietnam were able to reduce As concentrations in drinking water to <10 μg/L with a 40% success rate [39]. Laboratory and field experiments with household zero-valent Fe filters showed similarly effective results in Bangladesh [40]. Field tests with Fe-coated zeolites showed 99% removal of As [41]. Mn-Fe oxide-coated diatomites were able to reduce iAs from >40 μg/L to <10 μg/L in a pilot field-scale study [42]. Although many of these studies show promise for field-scale remediation of As, the biogeochemistry, aqueous geochemistry, and mineralogy of the surface and groundwater can affect the efficiency of the methods; therefore, it is necessary to understand the As adsorption process more thoroughly, so that the remediation methods can be applied more effectively.

1.2. Arsenic Treatment Methods

Due to the prevalence of As contamination worldwide and the threat of arsenicosis [43], research and the development of methods to understand As chemistry and to attenuate As in water are imperative [44]. Methods that have been developed to attenuate As such as electrocoagulation and electrodialysis [45,46], treatment with microorganisms to affect the biogeochemical cycling of As [7], and the adsorption of As onto a variety sorbents [47,48,49]. Among the sorbents used for As attenuation are organic polymers [50], and minerals such as dolomite [51], zeolites [52], and Al minerals such as alumina or gibbsite [53,54]. Alumina has been found to effectively remove iAs^V [55], but it is necessary to use activated alumina to efficiently remove iAs^III from solution [56]. A multitude of studies on As adsorbents have been conducted using Fe-based mineral sorbents [57,58], such as magnetite [54,59,60,61,62,63], ferrihydrite [53,54,64,65,66], goethite [48,54,66,67], and zero-valent iron [62].

The focus of the current work is on the chemistry of iAs species adsorbed to Fe-oxide and Fe-hydroxide mineral surfaces. Prior experiments have found that inorganic [64] and organic [64,68,69,70] ligands may adsorb to Fe minerals and compete with As for adsorption sites, but Zn cations [59] may augment As adsorption to Fe sorbents. Therefore, if ligands that inhibit As adsorption can be removed or precipitated and ligands that enhance As adsorption can be added to As-containing site, it could be possible to develop improved As-remediation techniques. Arsenic adsorption may also be enhanced by the addition of magnetite to agricultural waste such as wheat straw [61], but adsorption of As can be reduced by microorganisms [7,66]. If it is possible to control the growth and metabolism of microorganisms present in As-containing water, then it could be possible to enhance As adsorption.

1.3. Studying As Adsorption with Experimental and Modeling Methods

Because of the complex and often uncharacterized biogeochemistry of arsenic-rich environments, it is useful and necessary to use a variety of experimental data and modeling results to characterize these complex matrices. For example, Figure 1 shows models of the monodentate, bidentate, and outer-sphere adsorption of HAsO₄²⁻ to Fe-(oxyhyr)oxide clusters (Figure 2A–C) and periodic models (Figure 2D–F) models. Experimental and modeling techniques have been employed to determine the geometries, energetics, and spectroscopic properties exhibited by As species adsorbed to mineral surfaces; improved knowledge about the adsorption chemistry of As species can aid in the development of methods for attenuating As in the environment. Although the focus of our work is on the use of quantum mechanics (QM) techniques for studying the properties of As adsorption, it is necessary to frame these results in relationship to experimental and other modeling methods, because QM results can be useful for interpreting experimental data and for parameterizing other modeling methods such as classical force fields [71] and surface complexation models [72,73].

Figure 1. (A) Ferric iron (oxyhydr)oxide cluster model of Fe₂(OH)₆(H₂O)₄·4H₂O; (B) monodentate mononuclear (MM) cluster model of Fe₂(OH)₅(H₂O)₄H₂AsO₄·4H₂O; and (C) binuclear bidentate (BB) cluster model of Fe₂(OH)₄(H₂O)₄HAsO₄·4H₂O.

Figure 1. (A) Ferric iron (oxyhydr)oxide cluster model of Fe₂(OH)₆(H₂O)₄·4H₂O; (B) monodentate mononuclear (MM) cluster model of Fe₂(OH)₅(H₂O)₄H₂AsO₄·4H₂O; and (C) binuclear bidentate (BB) cluster model of Fe₂(OH)₄(H₂O)₄HAsO₄·4H₂O.

Figure 2. (A) and (D) monodentate mononuclear adsorption of HAsO₄²⁻; (B) and (E) bidentate binuclear adsorption of HAsO₄²⁻; and (C) and (F) outer-sphere complex of HAsO₄⁻, on clusters of Fe³⁺-(oxyhydr)oxide (A–C) and periodic α-goethite (010) (D–F) models.

Figure 2. (A) and (D) monodentate mononuclear adsorption of HAsO₄²⁻; (B) and (E) bidentate binuclear adsorption of HAsO₄²⁻; and (C) and (F) outer-sphere complex of HAsO₄⁻, on clusters of Fe³⁺-(oxyhydr)oxide (A–C) and periodic α-goethite (010) (D–F) models.

1.4. Studying As Adsorption with Experiments

Vibrational spectroscopy (e.g., Fourier transform infrared spectroscopy (FTIR) or Raman), X-ray absorption near edge structure (XANES), extended X-ray adsorption fine structure (EXAFS) spectroscopies, as well as adsorption isotherm and kinetics experiments are useful for determining the chemistry of As adsorption. FTIR and Raman studies are useful for determining the bonding configurations between As species and mineral surfaces. When As adsorbs as an inner-sphere complex, characteristic vibrational frequency shifts are observable [74,75,76,77,78,79]. XANES spectroscopy can provide information about the oxidation state of As that is adsorbed to mineral surfaces [80,81,82,83,84], and can determine if the oxidation state of As or the surface changes during the adsorption process [78,85]. EXAFS spectroscopy provides information about the coordination chemistry of adsorbed As [46,80,81,82,86,87,88,89,90,91,92]. Coordination state information is useful for determining whether the As adsorbs as a monodentate, bidentate, or outer-sphere complex (Figure 2). Moreover, studies that use two or more instrumental methods such as XANES/FTIR [78] or XANES/EXAFS [82] can increase the reliability of data interpretation. Furthermore, kinetics and isotherm studies provide information about the rates and energetics of As adsorption onto mineral surfaces that further help constrain adsorption mechanisms [55,93,94,95,96,97].

Although experimental techniques have contributed greatly to the understanding of As adsorption chemistry, computational chemistry can be used to help interpret experimental data on the mechanisms of As adsorption to mineral surfaces. In addition, computational chemistry can fill in missing information on the details of adsorption mechanisms and kinetics.

1.5. Studying As Adsorption with Mathematical Models

Surface complexation modeling (SCM) techniques such as the charge distribution multi-site complexation (CD-MUSIC) [76,95,98,99,100,101,102,103,104,105,106], the extended triple-layer model [107], isotherm modeling [108,109], and the ligand and charge distribution model (LCD) [110,111] have been used to model As species in solution and adsorbed to mineral surfaces. In general, SCM models use experimental data or quantum mechanics results such as equilibrium constants, bond lengths, and surface charges to calculate the adsorption isotherms of As on Fe mineral surfaces. When integrated, the CD-MUSIC and LCD SCMs are able to model humic substances interacting with Fe surfaces and the effect they have on As adsorption [111,112]. SCM can be used to interpret interactions among As, mineral surfaces, and organic and inorganic ligands, as well as charge and protonation effects. However, the precision of these models depends on the quality of their parameterization that is obtained from experimental data, which can be difficult to interpret, or with QM results [72,73].

1.6. Studying As Adsorption with Quantum Mechanics Modeling Methods

Numerous studies applying QM, or specifically density functional theory (DFT) methods [113,114], use models similar to those shown in Figure 2A–C. For example, DFT methods have been used to study the thermodynamics [91,115,116,117], vibrational frequencies [75,118,119,120,121], kinetics [120,121,122], ligand effects [87], oxidation-reduction reactions [123], and the coordination of adsorbed As [90,91,117,119,124,125,126,127,128,129,130,131,132].

Many of the previous computational chemistry studies relied on cluster models such as those in Figure 2A–C, but some groups have used periodic models to capture the chemistry of the mineral surface sorbates more precisely [124,125,127]. In this paper, we report structures and kinetics comparisons for the results from both cluster models (Figure 2A–C) and periodic models (Figure 2D–F). These comparisons are necessary to determine if and how the results from the molecular cluster and periodic model calculations differ.

In addition to the use of cluster versus periodic models, other factors can affect the results obtained from DFT calculations. These factors include surface charge, hydration, model convergence criteria, and potentially the software used for the calculation. For example, prior studies used highly charged models [91] to study the adsorption of As to Fe clusters; however, the work that we are presenting used models for the surface clusters with neutral charges, because localized high charges are unlikely. Prior work using the implicit solvation using the conductor reaction field (COSMO) [133,134] suggests that a single water molecule can produce activation energy results that are more precise than the addition of multiple water molecules can [122]. Conversely, other studies have used anhydrous surfaces to model As adsorption to Fe clusters [91]. The work we present herein used both hydrated cluster models and hydrated periodic models in an attempt to model the natural environment of As adsorption more accurately and realistically.

As the processing power and speed of computers increases, the size and complexity of As adsorption models increases, as illustrated by two Fe cluster models from papers dating from 2001 and 2006 [90,129]. The Fe cluster model in the latter paper is larger and likely more realistic than the model in the former paper. Results from prior DFT calculations provide contradictory results about the coordination state of As adsorption. For example, there are DFT results that predict that monodentate [128], bidentate [90,91,119,129], or a mixture of As coordination states are occurring on Fe-mineral surface models [130]. Because both the experimental data and the DFT results are providing ambiguous information about As coordination state, further calculations and experiments are necessary to clarify this topic. Alternatively, As may adsorb in a variety of configurations depending upon which surfaces are present on a given mineral sample [135].

The convergence criteria of the energy minimization calculations for the models could also have an effect on the precision of the calculated results and the ability of the results to reproduce experimental data and provide insight about As adsorption chemistry. For example we use a minimization convergence criteria of 0.03 kJ/mol, whereas, for example, Sherman and Randall [91] used higher tolerance energy minimization criteria 5 kJ/mol; however, tighter convergence criteria again reflects the availability and evolution of computational resources.

In the studies presented here, the thermodynamics, geometries, and kinetics of inner-sphere iAs^III and iAs^V adsorbed as monodentate mononuclear (MM) or bidentate binuclear (BB) complexes to solvated Fe clusters were evaluated. The molecular cluster results show how hydration and the initial Fe cluster charge affect iAs^III and iAs^V adsorption for BB models, and the results compare the of adsorption energies of BB iAs^III and iAs^V onto hydrated neutral Fe clusters. The DFT calculations on the cluster models also compare the calculated As—Fe distances for the BB models with pertinent experimental observations. In addition, inner- and outer-sphere As-Fe complexes were used to determine the activation energies (ΔE_a) of the adsorption/desorption process; these calculations were performed using both molecular cluster models and periodic models of the α-FeOOH (goethite) (010) surface.

2. Methods

2.1. Applied Quantum Mechanics Background

The application of QM with quantum chemistry software allows one to calculate chemical properties such as thermodynamics, kinetics, molecular geometries, spectroscopic parameters, transition states structures that might not be experimentally observable, and potentially hazardous chemistries [136,137,138].

Quantum chemistry calculations begin with an initial input of Cartesian coordinates for the molecule of interest, and then these coordinates are subsequently allowed to change during energy minimization calculations. The bond lengths, bond angles, and dihedral angles of the model are systematically perturbed, followed by the calculation of the relative energy of the model. After each energy calculation, subsequent systematic perturbations of the model geometry take place until the force (F), where F = dE/dr and dE and dr are the change in the energy and change in the model coordinates, respectively, converges to at or near “zero” (i.e., “zero” is defined as the convergence criteria by the modeler). When F is zero, the model then resides at a stationary point on a potential energy surface (PES). The user of computational chemistry software can specify the criterion for energy minimization convergence [136,139].

Subsequent calculation of the vibrational frequencies for the model determines the second derivative of energy with respect to atomic coordinates (d²E/dr²). If the calculated vibrational frequencies are all real (positive), then the model is at a PES minimum; if one vibrational frequency is imaginary (negative), then the model is at a transition state or PES maximum; if the model exhibits > 1 imaginary frequency, then the model is unstable and a new input geometry should be used. Obtaining a PES minimum does not guarantee that the model is at a global minimum, only that the model is at a local minimum. Using multiple input models obtained from a conformational analysis can aid in the attainment of a globally minimized final geometry [136,139].

The calculation of the vibrational frequencies (d²E/dr²) also allows the calculation of thermodynamic properties such as enthalpy, Gibbs free energy, and entropy. If one imaginary frequency is present, the model is at a transition state and the results from this model, the initial structure of the model, and the PES minimum model can be used to calculate rate constants for reactions. Note that throughout the manuscript the output from QM calculations is referred to as results and not data, we refer to the output from experiments as data. The vibrational frequency calculation provides infrared and Raman frequencies for the model, and further calculations can provide NMR chemical shifts, UV-Visible wavelengths, isotope effects, and temperature effects [136].

There are numerous methods available for calculating energies and other chemical properties with QM, among the most widely used are Hartree-Fock (HF) method [140], Møller-Plesset perturbation (MP) theory [141], and density functional theory (DFT) [113,114]. All of these methods arose from the development of quantum mechanics and the Schrödinger equation (ĤΨ(r) = EΨ(r)), where Ĥ is a Hamiltonian operator, Ψ(r) is the wave function, and E is the energy of the model. If Ψ(r) is known for a model, it is possible to solve for E and thus obtain the molecular properties for any model of interest. The Schrödinger equation has only to solve for the electronic energy of the model, because the nuclear energy and positions are relatively low and stationary compared to those of the electrons [142]. Thus far, it has been impossible to solve the Schrödinger equation for models that are more complex than H₂ because Ψ(r) makes the solution of the equation untenable for larger molecules with a greater number of electrons; therefore, it has been necessary to develop approximations for the Schrödinger equation such as HF and MP theory.

HF theory minimizes the energy of each electron iteratively with respect to the average energy of the other electrons in a model [140]. The shortcoming HF is that it does not account for electron correlation (repulsion) between each electron in the model. Using HF leads to an overestimation of model stability. MP theory accounts for electron correlation by systematically perturbing the molecular Hamiltonian [141]; however, the cost of using MP theory is prohibitive for models of geochemical interest. Furthermore, unlike HF theory, MP theory is not variational [143], so, the calculated energy could be lower than the ground state energy.

Unlike HF and MP theories that calculate electron interaction to obtain molecular energies, DFT calculates the electron density of the molecule to determine the energy [113,114,144]. The shortcoming of DFT is that the theory lacks a method for calculating the exact energy of the electron correlation and exchange term (E_xc). The inability to calculate the E_xc results in an underestimation of the total energy. Neglecting E_xc, as HF theory does, or underestimating E_xc, as DFT does, results in an underestimation of the total energy of a given model. For DFT, the lack of a precise correlation energy results because DFT does not account for the columbic interaction (repulsion) between electrons with anti-parallel (opposite) spin. Imprecise exchange energy results because DFT does not account precisely for the fact that electrons with parallel (same) spins cannot reside in the same orbital. Therefore, neglecting or underestimating the E_xc violates the Pauli Exclusion Principle. However, a variety of DFT methods have been developed to approximate E_xc and these methods can produce precise results [138,145]. Significantly, the computational cost of DFT calculations is substantially less than the cost of MP theory calculations, and the precision of DFT calculations is greater than that of HF calculations. Therefore, although DFT methods account imprecisely for E_xc, they are preferable to HF and MP methods.

Many DFT methods are available, and a particular method could be useful for calculating a particular chemical property (e.g., energy) but not for calculating other properties (e.g., NMR chemical shifts); these differences in precision are due to the parameters used for E_xc in the DFT methods [138,145]. Additional work is necessary to evaluate the efficacy of DFT methods for calculating properties such as adsorption energies, rate constants, and structures. Although DFT methods such as B3LYP [146,147] can provide accurate results for a variety of chemical properties, we suggest that computational geochemists begin to explore the use of other DFT methods that could provide improved results.

In addition to using an electron correlation method such as MP or DFT, it is necessary to specify a basis set when using molecular orbital (MO) calculations [136,148,149,150,151]; however, when using planewave calculations [152,153,154,155,156,157,158], basis sets are not used.

For DFT calculations on the clusters, the basis sets are equations that define atomic orbitals and are used in linear combinations to create molecular orbitals. When using basis sets, each atomic orbital for a given atom contains one electron [149], so for the C atom, which has six electrons and an electronic configuration of 1s²2s²2p², it is necessary to have a minimum of five basis functions (i.e., 1s, 2s, 2p_x, 2p_y, and 2p_z). Note that because the p-orbitals are energetically degenerate, the 2p_x, 2p_y, and 2p_z are present as basis functions. The variational principle of quantum mechanics states that E_g ≤ ⟨Ψ(r)∣Ĥ∣Ψ(r)⟩, where E_g is the ground state (lowest) energy of the molecule. The variational principle shows that increasing the accuracy of Ψ(r) will increase the calculated accuracy E relative to E_g; increased basis set size can increase the accuracy of Ψ(r).

Basis set size can be increased by using double zeta (DZ) or triple zeta (TZ) basis sets which double or triple the number of basis sets used, relative to the minimal basis set [148]. For C, the DZ and TZ basis sets have ten and fifteen basis functions, respectively. Another method for increasing the accuracy of a basis set is to use split basis sets, where more basis sets are used for the bond-forming valence electrons, while the non-bond-forming core electrons are treated with minimal basis sets. For C, this would involve one 1s orbital and two each of the 2s, 2p_x, 2p_y, and 2p_z basis sets for a total of nine basis sets.

Further increases in basis set size and precision can be obtained by the addition of polarization and diffuse basis functions [148]. Polarization functions increase the angular momentum basis sets on a particular atomic orbital; for example, in methane each H atom has a 1s orbital, but in the C-H bonds, the H atoms receives some p-orbital character from the C atom; therefore, including p-orbital polarized functions on the H atoms increases the accuracy of the orbital description. In addition, diffuse basis sets are used increase the electronic radius where electrons can reside. Diffuse basis sets are useful for calculations with anions and for weak interactions such as van der Waals interactions.

2.2. Molecular Orbital Theory Calculations with Fe Clusters

For the cluster model DFT calculations, all models were constructed in Materials Studio (Accelrys Inc., San Diego, CA, USA) and the energy minimization, Gibbs free energy, and transition-state calculations were performed in the gas phase using Gaussian 09 software [136]. All energy minimization calculations on the cluster models were performed without symmetry or atomic constraints. Energy minimizations, frequency, and kinetics (transition state) calculations were performed using the hybrid density functional B3LYP with the 6-31G(d) basis set [148,149,150,151]. B3LYP accounts for E_xc and the 6-31G(d) basis set is a DZ basis set with p-polarization functions on the non-H atoms. Energy convergence was set to 0.03 kJ/mol during the energy minimization calculations. The frequency calculations using B3LYP/6-31G(d) ensured that each model was at either a potential energy minimum (no imaginary frequencies) or at a transition state (one imaginary frequency) [136]; however, the frequency calculation does not ensure the model is at a global energy minimum.

For Gaussian calculations, it is necessary to specify an electron correlation method (e.g., B3LYP), a basis set, the type of calculation (e.g., optimization), the Cartesian coordinates of the atoms in the model, the charge of the model, and the spin multiplicity of the model [136]. The electron configuration of Fe³⁺ is [Ar]3d⁵. For the energy minimization calculations, we used high-spin Fe³⁺, where each 3d electron occupies one of the five d-orbitals, and where the electrons are all either spin up or spin down; this means that each of the two Fe atoms in the cluster has five unpaired electrons. The multiplicity is = 2S + 1, where S is the spin of an unpaired electron and can be +½ or –½. Therefore, for our high-spin clusters there are 10 unpaired electrons, each having a spin of +½, so the multiplicity = 2(10/2) + 1 or 11. For the rate constant calculations, we experimented by using high-spin multiplicity for both Fe atoms (i.e., 11), and a combination of up spin for one Fe atoms (i.e., +5) and down spin for the other (i.e., −5).

Fe model surface complex clusters were designed to minimize surface charge because charge buildup on the actual mineral surfaces is believed to be relatively small. However, surface models with charges were also calculated to demonstrate the effect of model charge on calculated energetics, and to compare with prior studies that used charged models [91]. In addition, explicit hydrating H₂O molecules were included for the aqueous species, the surface Fe-OH groups, and the adsorbed arsenic acid (iAs^V) molecules. Figure 1 shows examples of a tetrahydrated Fe (oxyhydr)oxide cluster, (Fe₂(OH)₆(H₂O)₄, A, without adsorbed HAsO₄²⁻, B, monodentate mononuclear (MM) adsorbed HAsO₄²⁻, and C, bidentate binuclear (BB) adsorbed HAsO₄²⁻. Hydration can be a key in accurately predicting structures and frequencies of anionic surface complexes because they are both a function of the hydration state of the sample [75,157]. The interatomic distances calculated in each model surface complex were compared to As-O bond lengths and As—Fe distances obtained from EXAFS spectra [46,80,81,82,86,87,88,89,90,91,92].

Gibbs free energies (G) of each species were estimated by calculating G in a polarized continuum the permittivity for water of 78.4 at 298.15 K. The integral equation formalism of the polarized continuum model (IEFPCM) [159] was used to calculate the Gibbs free energy in solution. The polarized continuum places the model in the cavity of a field with a constant permittivity, in this case for water; this field is used as a proxy for the solvent of interest. Single-point energy calculations were performed with the B3LYP functional and the 6-31+G(d,p) basis set [148,149,150,151]; the TZ basis set with augmented functions on the non-H atoms, and d-polarized and p-polarized functions for the non-H atoms and H-atoms, respectively, was used in order to improve the accuracy of the energy that was calculated using B3LYP/6-31(d) during the energy minimizations; this is a standard practice. The ΔG_ads was then determined from stoichiometrically balanced reactions. Configurational entropy terms are neglected in this approach; hence, we emphasize that these are Gibbs free energy estimates. We do not expect the precision of ΔG_ads to be better than ±10 kJ/mol.

To illustrate the potential effect that a chosen DFT method can have on geometry and thermodynamics, we report the As—Fe, As-OH, As=O, and As-OFe bond distances for the Fe₂(OH)₄(OH₂)₄HAsO₄·4H₂O model and ΔG_ads for the reaction:

H₂AsO₄⁻·8H₂O + Fe₂(OH)₆(OH₂)₄·8H₂O → Fe₂(OH)₄(OH₂)₄HAsO₄·4H₂O + OH⁻·4H₂O + 5H₂O

(1)

For this reaction, we used either B3LYP, PBE0 [160,161,162], and M06-L [163] DFT methods to minimize each structure in the reaction with the 6-31G(d) basis set, and the 6-311+G(d,p) basis with the self-consistent reaction field (SCRF) IEFPCM and the solvent water to calculate the single-point energy of each structure. The PBE0 and M06-L methods were chosen to compare with the results from the often-used B3LYP method, because the PBE0 function was the method used for the periodic planewave calculations for this work, and the M06-L method was specifically parameterized for use with transition metals such as Fe [162].

For the transition-state calculations on the clusters, the outer-sphere complexes of AsO₄³⁻ on an Fe³⁺-(oxyhydr)oxide clusters were obtained by constraining one Fe—As distance and allowing all other atoms to relax. The constrained distance was increased incrementally, allowing for energy minimization of the system at each new distance. The reaction path was graphically visualized. The change in energy for the adsorption reactions (ΔE_ads) were inferred by using the total electronic energy plus the zero-point correction obtained from the inner-sphere frequency calculations [163].

2.3. Planewave Calculations Using α-FeOOH (010)

The starting configuration for the periodic bidentate, binuclear HAsO₄²⁻ on the goethite (α-FeOOH) (010) surface was taken from previous simulations of HPO₄²⁻ on the same surface [135]. Phosphate and arsenate structures and chemistries are similar, so this starting configuration is a reasonable approximation. An energy minimization was performed on this starting configuration to allow the atoms to relax as necessary for the As for P substitution. Energy minimizations were carried out with the lattice parameters constrained the experimental values (9.24 × 9.95 Ų) [164] and with a vacuum gap between surface slabs of 10 Å. The model stoichiometry was 24FeOOH, HAsO₄²⁻, 29H₂O, and 2H₃O⁺ (Fe₂₄O₈₃H₈₉As). The small model system size and the high percentage of H⁺ per H₂O molecules severely limits the realism of the model compared to experimental systems, so the results from these calculations should be considered exploratory of model system behavior rather than an accurate portrayal of arsenate adsorption thermodynamics and kinetics.

Projector-augmented planewave calculations [152,165] with the Perdew-Erzenhof-Burke exchange correlation functional [160] were performed with the Vienna Ab-initio Simulation Package (VASP 5.2) [153,154,155,156,157]. The PBE0 exchange correlation functionals were Fe pv (14 valence e⁻), O (6 valence e⁻), H (1 valence e⁻) and As (5 valence e⁻) as labeled in the VASP exchange correlation functional library. Energy cut-offs (ENCUT in VASP input files) of 500 eV and 400 eV were used for energy minimizations and molecular dynamics simulations, respectively. The precision of the self-consistent field calculation of electron density was (PREC = Accurate = 700 eV/ROPT = −2.5 × 10⁻⁴) for energy minimizations and (PREC = Medium = 700 eV/ROPT = −2.5 × 10⁻³) for molecular dynamics simulations. The PREC-flag determines the energy cutoff (ENCUT) when no value is given for ENCUT in the central input file of VASP, INCAR, and the ROPT-tag controls the real-space optimization. The lower accuracies of the molecular dynamics (MD) simulations were chosen for practical reasons. Thousands of configurations and their energies need to be calculated for MD simulations, so the less stringent electron density grid speeds up the energy calculation at each step. The assumption here is that although the MD simulations are less accurate, they are not dramatically in error for predicting atomic structure. Thus, the MD simulations can be used to relax the atomic positions to achieve an approximate configuration and energy, and then energy minimizations can be performed to obtain structures and energies that are more precise. Without the MD simulations, the likelihood of the energy minimizations becoming trapped in a local potential energy minimum is much greater.

Periodic DFT calculations were run with 1 k-point created with the Monkhorst-Pack mesh. The DFT+U correction [166,167] was employed with a U = 4 eV for Fe and 0 eV for all other elements. Spin states were ordered according to the experimental observed magnetic ordering of goethite [168]. These selections have worked reasonably well in a previous study on goethite and goethite-water [73]. No dispersion-corrections were employed in these calculations although the results of DFT calculations of adsorption onto mineral surfaces may be affected by Van der Waals forces and how they are approximated [86]. The MD simulations were run at a temperature of 300 K maintained by the Nosé-Hoover thermostat [169]. Time steps were 0.5 fs (POTIM = 0.5). POTIM is a time-step variable and its value depends on the type of calculation being performed. Note that because some DFT methods may over structure and freeze water at 300 K, some authors have used higher temperatures to overcome this problem [170]. Another method is to use D instead of H [171], so that a 1 fs time step can be used instead of a 0.5 fs time step. Both are accepted practices, but we prefer to use the actual temperature and H atoms. This may cause error, but these errors are intrinsic to the method and as such not different in character from other computational uncertainties. Introducing errors by giving the atoms extra kinetic energy or mass may mask the discrepancies with experiment. Instead, highlighting these discrepancies is important and points to the need to improve computational methods.

Calculation of the periodic surface complex structures also allows for creation of more realistic molecular clusters that are surface specific. Figure 3 shows how an extended cluster (Figure 3B) can be extracted from the periodic model (Figure 3A). By selecting all the O atoms bonded to the Fe atoms connected to the As surface complex, a molecular cluster is created that retains surface-specific structure. The O atoms at the edge of the cluster are then terminated with H⁺ in order to satisfy valence and adjust the overall charge of the cluster as desired. Positions of hydration H₂O molecules can be included in the extended cluster (Figure 3B) to better mimic the aqueous phase. The combination of periodic and molecular cluster DFT results can take advantage of the strengths of each approach. For example, the periodic calculations should provide more accurate surface structures and adsorption energies, but molecular cluster models can be used to predict IR, Raman, and NMR spectra [135].

Figure 3. (A) Periodic model of goethite and (B) an extended cluster extracted from the periodic model.

Figure 3. (A) Periodic model of goethite and (B) an extended cluster extracted from the periodic model.

3. Results and Discussion

The molecular cluster results show how the initial cluster charges affect Gibbs free energy of adsorption (ΔG_ads) of iAs^V and iAs^III, the effect of hydration on ΔG_ads for the iAs^V models, and compare the of adsorption energies of triprotic iAs^III, and monoprotic or diprotic iAs^V onto hydrated neutral Fe clusters. The molecular cluster calculations also compare the calculated As—Fe distances with EXAFS data. In addition, calculations to determine the rate constant for iAs^V adsorption on Fe clusters and periodic models are discussed.

3.1. Effect of Cluster Charge on ΔG_ads

Reactions (1)–(7) in Table 1 show the ΔG_ads of H₃AsO₃ (Reaction (1) and (2)), HAsO₄²⁻ (Reaction (3) and (4)), and H₂AsO₄⁻ (Reaction (5)–(7)) onto either a neutral (Fe₂(OH)₆(OH₂)₄⁰) or +4 charged (Fe₂(OH)₂(OH₂)₈⁴⁺). Reaction (7) is stoichiometrically equivalent to one previously reported; however, for our calculation we used an energy minimization convergence criteria of 0.03 kJ/mol, whereas Sherman and Randall [91] used 5 kJ/mol. Significantly, Sherman and Randall [91] reported that the reaction:

Fe₂(OH)₂(OH₂)₆H₂AsO₄³⁺ (BB) + H₂O → Fe₂(OH)₂(OH₂)₇H₂AsO₄³⁺ (MM)

(2)

was endothermic and required +95 kJ/mol of energy; however, our results predict that this conversion would require +17 kJ/mol of Gibbs free energy. Both calculations predict that the BB structure is energetically favorable, but our results show energy difference between the BB and MM structures is not as large, and that these structures could co-exist in nature. The possibility of the presence of both BB and MM agree with prior research [130], but the lower ΔG_ads of the BB is inconsistent with the claim that MM adsorption is dominant [128]. The methodologies used to calculate the conversion of BB to MM could account for the calculated energy differences; the methods differ by:

• Model convergence criteria;
• Implicit solvation (our work) and gas-phase results from Sherman and Randall [91];
• Electronic energies (ΔE) from Sherman and Randall [91] and ΔG for our work;
• Use of a single, gas-phase minimized H₂O model to balance Equation (2) [91], whereas we used 1/8 the energy of an implicitly solvated model with eight H₂O molecules.

Throughout this work, unless otherwise noted, we used explicitly hydrated models for all of the reactants and products for the energy minimization calculations, and those same models for the implicitly solvated (i.e., IEFPCM) single-point energy calculations.

Table 1. Effect of Fe cluster charge on the ΔG_ads of iAs^III and iAs^V.

Reaction #	Reaction	ΔGads (kJ/mol)
(1)	H3AsO3 + Fe2(OH)6(OH2)40 → Fe2(OH)4(OH2)4HAsO3 + 2H2O	−60
(2)	H3AsO3 + Fe2(OH)2(OH2)84+ + 16H2O → Fe2(OH)2(OH2)6HAsO32+ + 2H3O+·8H2O	−159
(3)	HAsO42− + Fe2(OH)6(OH2)40 + 8H2O → Fe2(OH)4(OH2)4HAsO4 + 2OH−·4H2O	+14
(4)	HAsO42− + Fe2(OH)2(OH2)84+ → Fe2(OH)2(OH2)6HAsO42+ + 2H2O	−263
(5)	H2AsO4− + Fe2(OH)6(OH2)40 + 3H2O → Fe2(OH)4(OH2)4HAsO4 + OH−·4H2O	−309
(6)	H2AsO4− + Fe2(OH)2(OH2)64+ + 7H2O → Fe2(OH)2(OH2)6HAsO42+ + H3O+·8H2O	−336
(7)	H2AsO4− + Fe2(OH)2(OH2)84+ → Fe2(OH)2(OH2)6H2AsO43+ + 2H2O	−338

Table 1. Effect of Fe cluster charge on the ΔG_ads of iAs^III and iAs^V.

Reaction #	Reaction	ΔGads (kJ/mol)
(1)	H3AsO3 + Fe2(OH)6(OH2)40 → Fe2(OH)4(OH2)4HAsO3 + 2H2O	−60
(2)	H3AsO3 + Fe2(OH)2(OH2)84+ + 16H2O → Fe2(OH)2(OH2)6HAsO32+ + 2H3O+·8H2O	−159
(3)	HAsO42− + Fe2(OH)6(OH2)40 + 8H2O → Fe2(OH)4(OH2)4HAsO4 + 2OH−·4H2O	+14
(4)	HAsO42− + Fe2(OH)2(OH2)84+ → Fe2(OH)2(OH2)6HAsO42+ + 2H2O	−263
(5)	H2AsO4− + Fe2(OH)6(OH2)40 + 3H2O → Fe2(OH)4(OH2)4HAsO4 + OH−·4H2O	−309
(6)	H2AsO4− + Fe2(OH)2(OH2)64+ + 7H2O → Fe2(OH)2(OH2)6HAsO42+ + H3O+·8H2O	−336
(7)	H2AsO4− + Fe2(OH)2(OH2)84+ → Fe2(OH)2(OH2)6H2AsO43+ + 2H2O	−338

Comparing Reactions (1) with (2), (3) with (4), and (5) with (6) shows that the adsorption of iAs onto the more highly charged surfaces is more energetically favorable. However, it is unlikely that a +4 localized charged would occur in nature, and the results for ΔG_ads for the 0 charged Fe clusters are more realistic for exergonic adsorption of H₃AsO₃ and the endergonic adsorption of HAsO₄²⁻. The adsorption of H₂AsO₄⁻ is energetically favorable, regardless of the initial Fe cluster charge (Reactions (5)–(7)).

3.2. Effect of Fe Cluster Hydration on ΔG_ads for Anhydrous and Octahydrated H₂AsO₄⁻

Reactions (8)–(13) in Table 2 give examples of the effect of neutral Fe cluster hydration on the ΔG_ads of H₂AsO₄⁻. The reactions differ only by the number of H₂O molecules of hydration that are present on the initial Fe cluster (i.e., 0, 4, 8, or 8 for Reactions (8)–(11), respectively) and the HAsO₄²⁻-Fe cluster (i.e., 0, 4, 4, and 8, for Reactions (8)–(11), respectively). Reactions (12) and (13) used octahydrated H₂AsO₄⁻ as the reactant. Note that we used a tetrahydrated hydroxide model to mass and charge balance Reactions (8)–(11).

Table 2. Effect of Fe cluster hydration (Reactions (8)–(11)) and iAs^V hydration (Reactions (12) and (13)) on ΔG_ads of iAs^V.

Reaction #	Reaction	ΔGads (kJ/mol)
(8)	H2AsO4− + Fe2(OH)6(OH2)4 + 3H2O → Fe2(OH)4(OH2)4HAsO4 + OH−·4H2O	−186
(9)	H2AsO4− + Fe2(OH)6(OH2)4·4H2O + 3H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + OH−·4H2O	−195
(10)	H2AsO4− + Fe2(OH)6(OH2)4·8H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + OH−·4H2O + H2O	−217
(11)	H2AsO4− + Fe2(OH)6(OH2)4·8H2O + 3H2O → Fe2(OH)4(OH2)4HAsO4·8H2O + OH−·4H2O	−223
(12)	H2AsO4·8H2O + Fe2(OH)6(OH2)4·4H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + OH−·4H2O + 5H2O	−64
(13)	H2AsO4·8H2O + Fe2(OH)6(OH2)4·8H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + OH−·4H2O + 5H2O	−86

Table 2. Effect of Fe cluster hydration (Reactions (8)–(11)) and iAs^V hydration (Reactions (12) and (13)) on ΔG_ads of iAs^V.

Reaction #	Reaction	ΔGads (kJ/mol)
(8)	H2AsO4− + Fe2(OH)6(OH2)4 + 3H2O → Fe2(OH)4(OH2)4HAsO4 + OH−·4H2O	−186
(9)	H2AsO4− + Fe2(OH)6(OH2)4·4H2O + 3H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + OH−·4H2O	−195
(10)	H2AsO4− + Fe2(OH)6(OH2)4·8H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + OH−·4H2O + H2O	−217
(11)	H2AsO4− + Fe2(OH)6(OH2)4·8H2O + 3H2O → Fe2(OH)4(OH2)4HAsO4·8H2O + OH−·4H2O	−223
(12)	H2AsO4·8H2O + Fe2(OH)6(OH2)4·4H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + OH−·4H2O + 5H2O	−64
(13)	H2AsO4·8H2O + Fe2(OH)6(OH2)4·8H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + OH−·4H2O + 5H2O	−86

As the H₂O molecules of hydration increases from Fe₂(OH)₆(OH₂)₄ to Fe₂(OH)₆(OH₂)₄·8H₂O (Reactions (8)–(11)), the ΔG_ads of H₂AsO₄⁻ becomes more negative. Reactions (10) and (11) differ only by the numbers of H₂O molecules present on the cluster product, the HAsO₄²⁻-Fe cluster, four for Reaction (10) and eight for Reaction (11). The ΔG_ads for Reactions (10) and (11) differ by 6 kJ/mol, which is less than the ±10 kJ/mol error associated with thermodynamics calculations, so the results are indistinguishable. Using eight hydrating H₂O molecules on the HAsO₄²⁻-Fe cluster rather than four significantly increased the time need to minimize the model. Significantly, Reactions (12) and (13), where iAs^V is present as octahydrated H₂AsO₄⁻ exhibit ΔG_ads that are likely more realistic than those seen for ΔG_ads of anhydrous H₂AsO₄⁻ (Reactions (8)–(11)). Reactions (8)–(11) are shown here to emphasize the cluster hydration results. Therefore, we used the anhydrous H₂AsO₄⁻ reactant in Reactions (8)–(11) because these calculations are focusing on the hydration state of the clusters and are used here as a teaching tool. To attain results that are meaningful with respect to nature and experimental conditions, all of the products and reactants should be hydrated.

Recently [122], a claim was made that using a single explicit H₂O molecule and implicit solvation with the self-consistent reaction field (SCRF) COSMO [133,134] could produce superior results to using multiple H₂O molecules. This argument is based on calculations that used small, monohydrated organic and inorganic molecules in the COSMO SCRF that showed better agreement with experiment when a single, rather than multiple H₂O molecules were used to hydrate the models [172]. However, the work that modeled iAs^V interacting with ferric hydroxide clusters assumed that the result from simple organic and inorganic molecules would be applicable for the cluster calculations [122]; they did not test models with more than one H₂O molecule.

One argument for the addition of multiple H₂O molecules is that the model would better approximate the aqueous environment in which As adsorption occurs. However, if the multiple H₂O molecules are arranged in a way that does not lead to an observable PES minimum (i.e., becoming trapped in a local minimum), then using more than one H₂O molecule could lead to errors. On the other hand, the results for Reactions (12) and (13) suggest that multiple H₂O of hydration for all products and reactants could lead to results that are chemically more realistic than reactions with anhydrous reactants and products. Furthermore, chemical properties such as vibrational frequencies [75] are dependent on hydrogen bonding, so the inclusion of additional explicit H₂O molecules could be necessary to calculate precise spectroscopic and structural results.

Our work used the IEFPCM, and not the COSMO reaction field used by Farrell and Chaudhary [122]. A particular quantum chemistry method such as the SCRF, DFT method, or basis set can be useful for precisely calculating particular chemical properties, such as energies, but may produce other chemical properties that are imprecise. In this instance, the COSMO SCRF was parameterized to work most successfully with limited explicit hydration, but other SCRF such as IEFPCM may require the addition of more H₂O molecules to obtain precise results. Furthermore, particular DFT methods have been developed that are useful for calculating energies, geometries, and kinetics [138,173,174], whereas other DFT methods are useful for calculating spectroscopic properties such as NMR chemical shifts for H and C [175,176]. Because the exact E_xc for DFT is unknown, it is not yet possible to use one DFT method to calculate every chemical property. Therefore, it is necessary when doing DFT calculations to read the literature to find methodologies that are efficient for calculating the chemical properties of interest and to be willing to experiment with variations of those methods, if the calculated results are imprecise when compared with experimental data. This procedure is similar to those an experimentalist takes when deciding how to study a chemistry of interest.

3.3. Effect of As Oxidation State and DFT Method on ΔG_ads

Reactions (14)–(19) in Table 3 show the ΔG_ads of octahydrated H₃AsO₃, HAsO₄²⁻, and H₂AsO₄⁻ onto tetra- and octahydrated Fe₂(OH)₆(OH₂)₄⁰. The ΔG_ads for iAs^III is more favorable than it is for iAs^V; this is fortunate, if correct, because iAs^III is more toxic than iAs^V is. The results for the preferential adsorption of iAs^III over iAs^V at the point-of-zero charge for the Fe cluster are supported by experimental data that shows the same trend [101]. Under acidic or basic conditions, the Fe clusters would have charges, and neutral H₃AsO₃ might not adsorb to Fe surfaces as favorably. Conversely, although these reactions show less favorable adsorption of octahydrated HAsO₄²⁻ and H₂AsO₄²⁻ to the neutral Fe cluster than for iAs^III, under acidic or basic conditions when the cluster is charged, the charged iAs^V ions could adsorb more favorably than uncharged iAs^III. Although the products and reactants are anhydrous, Reactions (1)–(7) support this assertion, where neutral H₃AsO₃ adsorbs more strongly to the neutral cluster (Reaction (1)) than to the +4-charged cluster (Reaction (2)). Conversely, the charged iAs^V reactants in Reactions (3)–(7) adsorb more strongly to the +4-charged clusters than they do to the neutral clusters. These results that show weaker adsorption of H₃AsO₃ and stronger adsorption of charged iAs^V to charged Fe clusters are also supported by experimental data [101]. Furthermore, the results for Reactions (14) and (15) show that HAsO₄²⁻ would adsorb less favorable to the neutral Fe cluster than does H₂AsO₄⁻. Note that a tetrahydrated hydroxide model was used to mass and charge balance Reactions (14)–(17); the hydroxide was not necessary to mass charge balance Reactions (12) and (13).

Table 3. Effect of As oxidation state on BB adsorption to neutral Fe clusters. For Reaction (16), the density functional theory (DFT) methods used to calculate ΔG_ads were: ^a, B3LYP; ^b, PBE0; and ^c, M06-L.

Reaction #	Reaction	ΔGads (kJ/mol)
(12)	H3AsO3·8H2O + Fe2(OH)6(OH2)4·4H2O → Fe2(OH)4(OH2)4HAsO3·4H2O + 10H2O	−124
(13)	H3AsO3·8H2O + Fe2(OH)6(OH2)4·8H2O → Fe2(OH)4(OH2)4HAsO3·4H2O + 10H2O	−146
(14)	HAsO42−·8H2O + Fe2(OH)6(OH2)4·4H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + 2OH−·4H2O	+15
(15)	HAsO42−·8H2O + Fe2(OH)6(OH2)4·8H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + 2OH−·4H2O + 4H2O	−6
(16)	H2AsO4−·8H2O + Fe2(OH)6(OH2)4·4H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + OH−·4H2O + 5H2O	−64 a, −35 b, −3 c
(17)	H2AsO4−·8H2O + Fe2(OH)6(OH2)4·8H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + OH−·4H2O + 5H2O	−86

Table 3. Effect of As oxidation state on BB adsorption to neutral Fe clusters. For Reaction (16), the density functional theory (DFT) methods used to calculate ΔG_ads were: ^a, B3LYP; ^b, PBE0; and ^c, M06-L.

Reaction #	Reaction	ΔGads (kJ/mol)
(12)	H3AsO3·8H2O + Fe2(OH)6(OH2)4·4H2O → Fe2(OH)4(OH2)4HAsO3·4H2O + 10H2O	−124
(13)	H3AsO3·8H2O + Fe2(OH)6(OH2)4·8H2O → Fe2(OH)4(OH2)4HAsO3·4H2O + 10H2O	−146
(14)	HAsO42−·8H2O + Fe2(OH)6(OH2)4·4H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + 2OH−·4H2O	+15
(15)	HAsO42−·8H2O + Fe2(OH)6(OH2)4·8H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + 2OH−·4H2O + 4H2O	−6
(16)	H2AsO4−·8H2O + Fe2(OH)6(OH2)4·4H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + OH−·4H2O + 5H2O	−64 a, −35 b, −3 c
(17)	H2AsO4−·8H2O + Fe2(OH)6(OH2)4·8H2O → Fe2(OH)4(OH2)4HAsO4·4H2O + OH−·4H2O + 5H2O	−86

For Reaction (16), we compared the ΔG_ads calculated with the B3LYP, PBE0, or M06-L DFT methods. The ΔG_ads calculated with B3LYP (−64 kJ/mol), M06-L (−3 kJ/mol), and PBE0 (−35 kJ/mol) results all predict favorable, exergonic reactions. These results show that calculated thermodynamic results are dependent on the chosen DFT method. Because the adsorption of iAs^V is experimentally observed over a wide pH range [101], the results from B3LYP, PBE0, or M06-L could be correct. Thermodynamic results from DFT calculations are typically precise within ±10 kJ/mol; therefore, the B3LYP results would range from −74 to −54 kJ/mol and the PBE0 results would range from −45 to −25 kJ/mol; these results do not overlap, so the precision of the results from B3LYP and PBE0 differ. Within the ±10 kJ/mol error range of DFT methods, the M06-L results would range from −13 to +7 kJ/mol; therefore, because iAs^V adsorption is favorable, the M06-L results could be erroneous.

3.4. As—Fe Distance and As-O Bond Length Data from Experiments Compared with Cluster and Periodic Model Results

Table 4 shows the BB As—Fe distances and As-O bond distances calculated from this work for As^III and As^V, the calculated BB results of Sherman and Randall [91] for As^V, and EXAFS data for As^III and As^V on four mineral surfaces.

Notably, that the Fe₂(OH)₂(OH₂)₆H₂As^VO₄³⁺ (BB) model [91] exhibits four As-O bonds, two of which are As-OH bonds that are 1.73 Å. Sherman and Randall [91] observed a 1.62–1.64 Å As-O bond with EXAFS that is not present in their models where H₂AsO₄⁻ has two As-OH single bonds to the As atom. We argue that the 1.62–1.64 Å As-O bond is an As-O double bond (As=O) that our calculations predict to be 1.63 Å due to HAsO₄²⁻ being adsorbed to the Fe surface (model Fe₂(OH)₄(OH₂)₄HAs^VO₄ (BB)), rather than a As-OH single bond, and that iAs^V is not present on Fe surfaces as H₂AsO₄⁻ with two As-OH single bonds. However, we also note that the Fe₂(OH)₄(OH₂)₄HAs^VO₄ (BB) model is not hydrated with explicit H₂O molecules and that when the model is hydrated with either four or eight explicit H₂O molecules, the As=O bond length increases from 1.63 to 1.67 Å, which is slightly longer than the observed 1.62–1.64 Å As=O length, but is still shorter than the 1.73 Å bond lengths from the model of Sherman and Randall [91]. In Table 4, we report the results for the +3-charged model that was energy minimized using the methods described in this work. The As-OFe and As-OH bond distances calculated here differ by 0.01 Å from those reported by Sherman and Randall [91], whereas the As—Fe distances calculated here are both 0.05 Å shorter than those reported by Sherman and Randall [91]. The difference in As—Fe distances likely occur due to differences in methodology, but the results from both models lie within the range of experimental uncertainty.

Table 4. As—Fe interatomic distance, As-OFe bond distance, and As-O bond distance results from this work compared to previous calculations and extended X-ray adsorption fine structure (EXAFS) data for As^V and As^III adsorption onto ferrihydrite (Fh), hematite (Hm), goethite (Gt), and lepidocrocite (Lp). For Fe₂(OH)₄(OH₂)₄HAs^VO₄·4H₂O (BB), the DFT methods used to calculate ΔG_ads were: ^x, B3LYP; ^y, PBE0; and ^z, M06-L.

AsV Complex	As—Fe (Å)	As—Fe (Å)	As-OFe (Å)	As-OFe (Å)	As-OH (Å)	As-OH (Å)	As=O (Å)
Fe2(OH)4(OH2)4HAsVO4 (BB)	3.13	3.25	1.71	1.73	1.83		1.63
Fe2(OH)4(OH2)4HAsVO4·4H2O (BB)	3.20 x, 3.19 y, 3.08 z	3.28 x, 3.24 y, 3.29 z	1.69 x, 1.68 y, 1.69 z	1.72 x, 1.71 y, 1.72 z	1.76 x, 1.75 y, 1.79 z		1.67 x, 1.66 y, 1.65 z
Fe2(OH)4(OH2)4HAsVO4·8H2O (BB)	3.30	3.30	1.70	1.70	1.76		1.67
Goethite (010) periodic model (BB)	3.56	3.68	1.72	1.72	1.78		1.73
Fe2(OH)2(OH2)6H2AsVO43+ (BB)	3.24	3.24	1.70	1.70	1.72	1.72
Fe2(OH)2(OH2)6H2AsVO43+ (BB) a	3.29	3.29	1.71	1.71	1.73	1.73
AsV on Fh (BB) a	3.27	3.38	1.70	1.70	1.67		1.64
AsV on Gt (BB) a	3.30	3.30	1.70	1.70	1.70		1.63
AsV on Lp (BB) a	3.30	3.32	1.71	1.71	1.66		1.63
AsV on Hm (BB) a	3.24	3.35	1.70	1.70	1.70		1.62
AsV on Fh (BB) b	3.25 (±0.02)
AsV on Gt (BB) b	3.28 (±0.01)
AsV on Lp (BB) c	3.31 (±0.014)		1.69 (±0.004)
AsV on Gt (BB) c	3.30 (±0.008)		1.69 (±0.004)
AsV on Fh (BB) d	3.27
Goethite (010) periodic model (MM)	3.54	5.00 †	1.78	1.75	1.71 ‡		1.68
AsV on Gt e	3.25 §		1.689				1.679
AsIII Complex	As—Fe (Å)	As—Fe (Å)	As-O (Å)	As-O (Å)		As-O (Å)	As=O (Å)
Fe2(OH)4(OH2)4HAsIIIO3·4H2O (BB)	3.26	3.41	1.77	1.74		1.90	na
Fe2(OH)4(OH2)4HAsIIIO3·8H2O (BB)	3.29	3.39	1.78	1.72		1.90
AsIII on Lp (BB) c	3.41 (±0.013)		1.78 (±0.014)
AsIII on Gt (BB) c	3.31 (±0.013)		1.78 (±0.012)
AsIII on Fh (BB) d	3.41–3.44
AsIII on Fh and Hm (BB) f	3.35 (±0.05)
AsIII on Gt and Lp (BB) f	3.3–3.4
AsIII on Gt (BB) g	3.378 (±0.014)

Notes: ^a Sherman and Randall [91]; ^b Waychunas et al. [89]; ^c Farquhar et al. [82]; ^d Gao et al. [87]; ^e Loring et al. [128]; ^f Ona-Nguema et al. [92]; ^g Manning et al. [81]. ^† This As—Fe distance does not agree with that reported by Loring et al. [128] for a MM. ^‡ For the MM periodic model, there is one As-OH bond and two As partial double bonds (1.71 and 1.68 Å), because HAsO₄²⁻ is adsorbed to the surface. ^§ For the Loring et al. [128] model, both As-O bonds are aprotic and should have partial double bonds.

Table 4. As—Fe interatomic distance, As-OFe bond distance, and As-O bond distance results from this work compared to previous calculations and extended X-ray adsorption fine structure (EXAFS) data for As^V and As^III adsorption onto ferrihydrite (Fh), hematite (Hm), goethite (Gt), and lepidocrocite (Lp). For Fe₂(OH)₄(OH₂)₄HAs^VO₄·4H₂O (BB), the DFT methods used to calculate ΔG_ads were: ^x, B3LYP; ^y, PBE0; and ^z, M06-L.

AsV Complex	As—Fe (Å)	As—Fe (Å)	As-OFe (Å)	As-OFe (Å)	As-OH (Å)	As-OH (Å)	As=O (Å)
Fe2(OH)4(OH2)4HAsVO4 (BB)	3.13	3.25	1.71	1.73	1.83		1.63
Fe2(OH)4(OH2)4HAsVO4·4H2O (BB)	3.20 x, 3.19 y, 3.08 z	3.28 x, 3.24 y, 3.29 z	1.69 x, 1.68 y, 1.69 z	1.72 x, 1.71 y, 1.72 z	1.76 x, 1.75 y, 1.79 z		1.67 x, 1.66 y, 1.65 z
Fe2(OH)4(OH2)4HAsVO4·8H2O (BB)	3.30	3.30	1.70	1.70	1.76		1.67
Goethite (010) periodic model (BB)	3.56	3.68	1.72	1.72	1.78		1.73
Fe2(OH)2(OH2)6H2AsVO43+ (BB)	3.24	3.24	1.70	1.70	1.72	1.72
Fe2(OH)2(OH2)6H2AsVO43+ (BB) a	3.29	3.29	1.71	1.71	1.73	1.73
AsV on Fh (BB) a	3.27	3.38	1.70	1.70	1.67		1.64
AsV on Gt (BB) a	3.30	3.30	1.70	1.70	1.70		1.63
AsV on Lp (BB) a	3.30	3.32	1.71	1.71	1.66		1.63
AsV on Hm (BB) a	3.24	3.35	1.70	1.70	1.70		1.62
AsV on Fh (BB) b	3.25 (±0.02)
AsV on Gt (BB) b	3.28 (±0.01)
AsV on Lp (BB) c	3.31 (±0.014)		1.69 (±0.004)
AsV on Gt (BB) c	3.30 (±0.008)		1.69 (±0.004)
AsV on Fh (BB) d	3.27
Goethite (010) periodic model (MM)	3.54	5.00 †	1.78	1.75	1.71 ‡		1.68
AsV on Gt e	3.25 §		1.689				1.679
AsIII Complex	As—Fe (Å)	As—Fe (Å)	As-O (Å)	As-O (Å)		As-O (Å)	As=O (Å)
Fe2(OH)4(OH2)4HAsIIIO3·4H2O (BB)	3.26	3.41	1.77	1.74		1.90	na
Fe2(OH)4(OH2)4HAsIIIO3·8H2O (BB)	3.29	3.39	1.78	1.72		1.90
AsIII on Lp (BB) c	3.41 (±0.013)		1.78 (±0.014)
AsIII on Gt (BB) c	3.31 (±0.013)		1.78 (±0.012)
AsIII on Fh (BB) d	3.41–3.44
AsIII on Fh and Hm (BB) f	3.35 (±0.05)
AsIII on Gt and Lp (BB) f	3.3–3.4
AsIII on Gt (BB) g	3.378 (±0.014)

Notes: ^a Sherman and Randall [91]; ^b Waychunas et al. [89]; ^c Farquhar et al. [82]; ^d Gao et al. [87]; ^e Loring et al. [128]; ^f Ona-Nguema et al. [92]; ^g Manning et al. [81]. ^† This As—Fe distance does not agree with that reported by Loring et al. [128] for a MM. ^‡ For the MM periodic model, there is one As-OH bond and two As partial double bonds (1.71 and 1.68 Å), because HAsO₄²⁻ is adsorbed to the surface. ^§ For the Loring et al. [128] model, both As-O bonds are aprotic and should have partial double bonds.

For the iAs^V models, the two calculated As—Fe distances within each configuration differ by 0.12, 0.08, and 0.00 Å from each other for the Fe₂(OH)₄(OH₂)4HAs^VO₄, Fe₂(OH)₄(OH₂)₄HAs^VO₄·4H₂O, and Fe₂(OH)₄(OH₂)4HAs^VO₄·8H₂O models, respectively. Sherman and Randall [91] reported two As—Fe distances for ferrihydrite (Fh), lepidocrocite (Lp), hematite (Hm), and goethite (Gt), whereas the other data report a single As—Fe distance for the adsorption onto Fh, Gt, or Lp. The As^V—Fe distance results from the Fe₂(OH)₄(OH₂)₄HAs^VO₄·8H₂O model agree within experimental uncertainty with the Gt and Lp data of Sherman and Randall [91], the Gt and Fh data of Waychunas et al. [89], and the Gt and Lp data of Farquhar et al. [82]. The data for the As^V—Fe distances overlap for the minerals used for these studies; therefore, determining the type of mineral to which the iAs^V is bonding could be difficult; however, we can state that the Fe cluster models are predicting As—Fe distances that are indicative of BB adsorption of As^V. Similarly, the calculated experimental and As^V—OFe bond distances agree precisely.

For the BB results for the periodic structures of α-FeOOH (010), the As-OFe bonds (1.72 Å) show precise agreement with experiment and correlate well with the cluster results. However, the As—Fe distances (3.56 and 3.68 Å) both overestimate the experimental data, and the As-OH and As=O bonds both overestimate the experimentally observed bond distances. For the MM periodic model, the calculated As—Fe distance is overestimating the 3.25 Å distance measured distance of Loring et al. [128] by 0.29 Å. The errors associated with the periodic calculations could be due to systematic errors in the planewave calculations or could be because the α-FeOOH (010) surface is not the surface where As adsorption predominately occurs.

The two As^III—Fe bond distances in the Fe₂(OH)₄(OH₂)₄HAs^IIIO₃·4H₂O model differ by 0.15 Å and by 0.10 for the octahydrated version of that model., The longer calculated As—Fe distances (ca. 3.4 Å) agree with most of the As^III data within uncertainty, whereas the shorter calculated As—Fe distance agrees well with the data of Gao et al. [87]. Again, because of the data overlap and due to the uncertainty in the As—Fe distances observed for As^III adsorption onto Hm, Gt, Lp, and Fh, it is difficult to resolve different sorption mechanisms of various Fe-oxide and Fe-hydroxide minerals; however, it is possible to state the BB adsorption is occurring. Furthermore, it is possible to differentiate As^III—Fe and As^V—Fe BB adsorption, because the former distances are approximately 3.4 Å, whereas the latter are approximately 3.3 Å. Significantly, these distances are seen both experimentally and computationally. Moreover, there is good agreement between the calculated As^III-OFe distance and the EXAFS data of Sherman and Randall [91].

For the Fe₂(OH)₄(OH₂)₄HAs^VO₄·4H₂O (BB) model, we compared the As—Fe distance results obtained from the B3LYP, PBE0, and M06-L methods. For the As—Fe distances, the B3LYP (3.20 Å) and PBE0 (3.19 Å) results correlate for one of the As—Fe distance, and the B3LYP (3.28 Å) and M06-L (3.29 Å) results agree for the other As—Fe distance. The As-O and As=O bond lengths agree well among the DFT methods and agree precisely with the experimental data in Table 4. This type of DFT method testing helps eliminate the possible effects of exchange-correlation functional errors on the results.

The Fe₂(OH)₄(OH₂)₄HAs^VO₄·4H₂O (BB) model used for these DFT method comparisons is the adsorption product of Reaction (16), and the M06-L ΔG_ads results for Reaction (16) ranged from −13 to +7 kJ/mol, suggesting thermodynamic adsorption results from M06-L that are potentially unfavorable, relative to the results from B3LYP and PBE0 (Table 3). In addition, the 3.08 Å As—Fe distance from the M06-L minimized model (Table 4) is predicts is significantly shorter than the experimental data and the results from B3LYP and PBE0. Therefore, because the As—Fe distance calculated with M06-L is imprecise, it is likely that the ΔG_ads results from M06-L is also imprecise.

Notably, the PBE0 As—Fe distance results from the cluster and periodic calculations differed distinctly. The cluster results underestimated the As—Fe distance data by approximately 0.1 Å, whereas the periodic calculation results overestimated those data by approximately 0.25 Å. The PBE0 method, like many DFT methods, may contain different parameters, depending on which software package implements it (e.g., VASP, Gaussian 09, etc.), so the results obtain with a particular DFT method using different software packages might not be directly comparable. In addition to the potential differences in the DFT methods, model sizes could also contribute to the discrepancies in the calculated distances obtained from the periodic and cluster models. One would presume that the larger periodic models would provide results that are more precise relative to the data than the smaller cluster models do; however, neither model size produced precise As—Fe distances. Differences between periodic and cluster model results have been discussed previously [73,135].

3.5. Sorption Kinetics for iAsV on Cluster and Periodic Models

Calculations were completed to show a possible reaction pathways for desorption of the monodentate inner-sphere complex of HAsO₄²⁻ from Fe³⁺-(oxyhydr)oxide cluster and from the periodic goethite (010) surface. Although the bidentate binuclear complex is likely to be more stable [88,91], the monodentate species is an intermediate between the bidentate and outer-sphere species. Figure 4 shows desorption of iAsV from a Fe₂(OH)₄(H₂O)₅-HAsO₄ cluster model, where the model begins as a MM structure with a As—Fe distance of 3.27 Å (Figure 4A), moves through a transition state structure (Figure 4B), and reaches the outer-sphere structure (Figure 4C) where the As—Fe distance is 4.36 Å. The Fe—As distances were increased manually and then held constant in each calculation until there ceased to be a bonding interaction. The energies of the monodentate reaction pathway are portrayed as a function of the Fe—As distance in the Figure 4 for both periodic and cluster models.

Figure 4. Desorption of iAs^V from Fe clusters showing (A) the initial MM model, (B) the transition state model, and (C) the outer-sphere, final structure model.

Figure 4. Desorption of iAs^V from Fe clusters showing (A) the initial MM model, (B) the transition state model, and (C) the outer-sphere, final structure model.

Based on the results shown in Figure 5, ΔE_a for the breaking of the first bond in monodentate complex requires approximately +133 and +70 kJ/mol in the periodic and cluster models, respectively. (Note that there is a small increase in energy of the model system near a Fe—As distance of 4.2 Å, but this energy increase is insignificant compared to the first barrier.) The energy barrier for the reverse reaction is higher, +148 kJ/mol, in the periodic model because the outer-sphere complex is lower in energy in this model. The higher energy of the inner-sphere complex and the high-energy barrier of adsorption suggest that adsorption would not occur to the goethite (010) surface under these conditions; this result corroborates with the discussion the long As—Fe distances reported in the previous section for the periodic models. In the molecular cluster models, the adsorption reaction barrier is insignificant, i.e., +1 kJ/mol. We strongly remind the reader, however, that the conditions of this model are not realistic compared to experimental conditions where lower H⁺ activity, lower arsenate concentrations and greater volume of water would exist and affect the results. On the other hand, the cluster models exhibit almost no energy barrier to adsorption from outer-sphere to monodentate with the inner-sphere complex having a lower energy (Figure 5).

Figure 5. Constrained scan of Fe—As distances starting from monodentate configuration to outer-sphere using periodic and molecular cluster DFT calculations results in ΔE_a of adsorption of +148 and +1 kJ/mol, respectively and desorption of +133 and +70 kJ/mol, respectively.

Figure 5. Constrained scan of Fe—As distances starting from monodentate configuration to outer-sphere using periodic and molecular cluster DFT calculations results in ΔE_a of adsorption of +148 and +1 kJ/mol, respectively and desorption of +133 and +70 kJ/mol, respectively.

The discrepancies between the two types of models are illustrative of some of the problems that can be encountered by each type of approach. First, although the periodic models were run for short (i.e., 6000 steps × 0.5 fs = 3 ps) molecular dynamics simulations at 300 K to relax the atoms, the complex nature of the periodic model all but ensures that a global minimum configuration will not be obtained. This is an example of a general problem, i.e., adding more atoms to the simulation may make it more realistic but increases the number of potential energy minima dramatically. Thus, the transition state may overestimate the ΔE_a because the system is not in the lowest possible potential energy configuration, especially with respect to the configuration of H₂O molecules.

There are numerous possibilities for overcoming the metastable minimum problem. Longer simulation runtimes are one option. These longer simulations could be performed with tight-binding DFT (i.e., DFT-B; see REF for a review of DFT-B) or classical force fields. However, one problem with classical simulations is that it is difficult to created accurate parameterizations that allow for bond-making and bond-breaking, especially for configurations far from equilibrium such as transition states (see [71] for a review). Replicate MD simulations are another option for exploring configuration space. These require multiple simulations at different temperatures to be run simultaneously such that higher and lower temperature configurations can be switched with potential energies overlap. Again, however, this method requires a significant expansion of computational effort to run the multiple MD simulations.

Alternatively, the cluster model allows the “surface” atoms to relax without constraint from the remainder of the crystal, which may help explain the lower ΔE_a. The “outer-sphere” configuration is higher in energy in this case because the HAsO₄²⁻ is not completely solvated by the six extra H₂O molecules in the model. In addition to the loss of solvation energy, protonation of the HAsO₄²⁻ does not occur in the cluster whereas in the periodic model H₂AsO₄⁻ is the product.

Even with these discrepancies in the energies and products between the periodic and cluster models, the structures of the reactants and transition states are similar. A combined approach using insights from both periodic and cluster models is useful at this time because each has strengths and weaknesses. These first simulations of the reaction path can be reiterated by performing longer MD simulations at various steps and by searching for lower energy points along the reaction path. In addition, lower energy transition states determined via the cluster approach can be used to guide construction of transition states in the more realistic (but more complex) periodic simulations.

The main problem in this case, however, is likely to be that the (010) surface is not the dominant face responsible for adsorbing arsenate onto goethite. Thus, extensive calculations of any type could prove futile in terms of reproducing the observed ΔE_a of adsorption/desorption. Other surfaces could be examined, because the (010) surface may not be the most preferred surface for adsorption of arsenate onto goethite. We had selected the (010) surface for arsenate adsorption onto goethite based on the analogy with chromate adsorption [177,178]. Recent DFT calculations have been used to suggest that the (210) surface adsorbs phosphate more strongly than the (010) surface [135]. None of the three just-mentioned studies included arsenate, however, so we are using the other oxyanions as analogs for arsenate. Our future computational work will focus on arsenic sorption reaction mechanisms onto the (210) surface, and this work would benefit from experiments similar to those that used chromate [177]. This point emphasizes that realistic model construction is one of the most important considerations in performing computational geochemistry. Too often, missing components are some inaccuracy in the original model creation leads to discrepancies with observation that cannot be resolved with even the most accurate quantum mechanical calculations.

4. Conclusions

This work explored the effect of cluster charge, hydration, As oxidation state, and DFT methods on the Gibbs free energy of adsorption (ΔG_ads) of inorganic arsenic (iAs) species onto Fe³⁺-(oxyhydr)oxide models. In general neutral clusters and hydrated models produced ΔG_ads results that are likely more realistic than models with charged clusters and anhydrous models. As shown with experiments [101], iAs^III adsorption onto neutral Fe³⁺-(oxyhydr)oxide cluster models was more exergonic than iAs^V adsorption onto the same cluster models. For the DFT calculations on the clusters, the results showed that both ΔG_ads and As—Fe distances depend on the DFT method used to calculate those properties; however the As—Fe distance results from these calculations generally agreed precisely with the experimental data cited.

The cluster model As—Fe distance and As-O bond distance results showed relatively precise agreement with the experimental data. Conversely, the periodic planewave calculation results for iAsV adsorption onto α-FeOOH (010) generally overestimated the As—Fe distance and As-O bond length data for iAs^V adsorption onto goethite. Other α-FeOOH surfaces could produce results that are more precise.

Sorption kinetics calculations using DFT with cluster and a periodic model of α-goethite (010) showed discrepancies in the calculated activation energies of iAs^V adsorption. One major difference for the discrepant results could be that the relatively large periodic model did not reach an energy minimum during the DFT MD simulation, whereas the cluster model was smaller than the periodic model and did attain a PES minimum. Although the calculated activation energies for the two methods differed, the initial and transition state structures for both calculations were structurally similar. Longer DFT MD simulations and periodic structures other than the (010) surface of α-FeOOH could produce results that are more precise.

The calculated reaction rates, thermodynamics, and structural results presented in this work provide results that could lead to a better understanding of the adsorption of arsenic to Fe (oxyhydr)oxide minerals. However, further studies are necessary to better determine which DFT methods produce the most precise results, the effect of model size on model precision, and the effects of model hydration and surface charge on As adsorption to Fe (oxyhydr)oxide models. Furthermore, basis set size, which was not addressed herein, could potentially affect the precision of the results for the cluster models; therefore, future studies should include the evaluation of basis set effects. Furthermore, increased collaborative efforts among experimental and computational (geo)chemists could lead to improved knowledge about arsenic adsorption on Fe minerals.