Quantum-mechanical Methods for Quantifying Incorporation of Contaminants in Proximal Minerals

Incorporation reactions play an important role in dictating immobilization and release pathways for chemical species in low-temperature geologic environments. Quantum-mechanical investigations of incorporation seek to characterize the stability and geometry of incorporated structures, as well as the thermodynamics and kinetics of the reactions themselves. For a thermodynamic treatment of incorporation reactions, a source of the incorporated ion and a sink for the released ion is necessary. These sources/sinks in a real geochemical system can be solids, but more commonly, they are charged aqueous species. In this contribution, we review the current methods for ab initio calculations of incorporation reactions, many of which do not consider incorporation from aqueous species. We detail a recently-developed approach for the calculation of incorporation reactions and expand on the part that is modeling the interaction of periodic solids with aqueous source and sink phases and present new research using this approach. To model these interactions, a systematic series of calculations must be done to transform periodic solid source and sink phases to aqueous-phase clusters. Examples of this process are provided for three case studies: (1) neptunyl incorporation into studtite and boltwoodite: for the layered boltwoodite, the incorporation energies are smaller (more favorable) for reactions using environmentally relevant source and sink phases (i.


Incorporation Reactions Important in Low Temperature Mineralogy
The incorporation of atoms into a host phase is an important phenomenon in mineralogy and materials science.Understanding the fundamental thermodynamics and kinetics of the incorporation mechanism, including charge-balanced substitutions, final incorporated geometries, and the thermodynamic stability of substituted phases, is of importance for many geochemical and mineralogical studies including biomineralization, acid mine drainage, and long-term predictions of geologic alteration of nuclear waste.
Incorporation of contaminants into minerals is a primary mechanism for immobilization of radionuclides at a contaminated site (e.g., Hanford Site in Richland, Washington) or potential immobilization after failure of a geologic repository (e.g., Waste Isolation Pilot Plant in Carlsbad, New Mexico).The study of incorporation of radionuclides will enable more refined risk assessment models for long-term evaluation of geologic repositories.Risk assessment of a geologic repository requires an understanding of the breakdown of natural and engineered barriers leading to the transport of radionuclides.Significant mineral sink phases occur as a result of the oxidation of the fuel itself, oxidation of engineered spent fuel canisters, occur in the surrounding geology, or are added to a geologic repository as a backfill.Table 1 lists minerals available as a sink for radionuclides along with the source of the mineral.The examples used in the current review are related to immobilization of radioactive waste, particularly actinides.Rutile TiO 2 Natural occurrence, ceramic waste forms Tc, Ru [12] Notes: a This item comprises the sulfate minerals anhydrite (Ca), barite (Ba), anglesite (Pb), celestine (Sr); b Carbonates with aragonite structure: aragonite (Ca), witherite (Ba), cerussite (Pb), strontianite (Sr).

Quantum-Mechanical Incorporation Calculations
Quantum-mechanical calculations can be used to identify atomistic scale phenomena that impact incorporation, including coupled substitutions and charge transfer.The use of quantum-mechanical calculations to study low-temperature geochemical and mineralogical phenomena has been ever increasing since the implementation of density functionals, such as the local density approximation (LDA) and the generalized gradient approximation (GGA) [13], proving the Hohenberg-Sham theorem that the total energy of a system can be determined according to the electron density.Taking advantage of Bloch's theorem and the periodicity of mineral structures, density functional theory (DFT) calculations using planewave basis sets and pseudopotential approximations have become common place in computational mineralogy [14].However, dealing with incorporation at a quantum-mechanical level is challenging due to the combination of aqueous and solid systems.That is, the incorporated species typically comes from an aqueous solution, rather than a solid phase, while the final doped mineral is a solid.Thus, the subject of incorporation (i.e., the charged, aqueous source species) would be best treated using a quantum-chemical cluster model, whereas the solid host phase is best treated using periodic boundary conditions.A similar challenge is posed by the reaction products which are the solid periodic host with the incorporated molecule or ion and the, typically charged, aqueous sink phase of the species that was released from the host mineral into the surrounding water.The challenge is, then, how to accurately model an incorporation reaction such that a uniform computational theory is used throughout the incorporation reaction.Details regarding computational challenges for calculating accurate incorporation energies are in Sections 2.1 and 4.
Incorporation energies are the reaction energy describing the incorporation reaction, where the final incorporation may range from a slight impurity to a complete solid solution.For the latter, the incorporation does not result in any structural changes sans modification of the unit cell parameters according to Vegard's law (e.g., [15]).In this case, the Gibbs free energy of mixing can be considered the incorporation energy.More traditionally, incorporation energies refer to the incorporation of an impurity into a bulk mineral phase.For example, Se is associated with pyrite and arsenopyrite [16]; thus, the energy associated with the "defect" Se in the Fe-sulfide is the incorporation energy.
From quantum-mechanical calculations, the kinetics of the incorporation (i.e., the precise pathway of the incorporated ion) is not considered.For example, this type of calculation does not consider whether a cation sorbs on the mineral surface and diffuses into the mineral or sorbs on the surface and co-precipitates with the mineral, with the latter being a more likely scenario.Rather, the computational methodology presented here, allows one to evaluate the overall thermodynamic stability of the final incorporated phase.In contrast, the experimental incorporation limit may be dependent on the dominating process of the ones described above, because, e.g., co-precipitation maybe far from equilibrium with difficulty reaching equilibrium subsequently by diffusing out excess incorporation cations.
The incorporation energy is influenced by the mechanism, e.g., the incorporation site.On the experimental side, one might assume that the incorporated ion or molecule would find its most favorable position; however, computationally, different incorporation mechanisms/sites must be assessed separately due to the static nature of the quantum-mechanical incorporation calculations.For example, the incorporation of Ti into zircon, an important reaction for accurate thermochronological dating, is dependent on the location of the incorporated Ti.That is, for pressures below 3.5 GPa, Ti substitution at the Si site is more thermodynamically favored over substitutions at the Zr site and vice versa at higher pressures [17].On the atomic-scale, incorporation of Ti into zircon is significantly altered by pressure conditions and corrections are required for the Ti-in-zircon geothermometer.
Substitutions of cations with different oxidation states than the host require a charge-balancing mechanism.A common mineralogical example is the plagioclase feldspars (M x+ (Al x ,Si 4−x )O 8 , M = Na + , Ca 2+ ), where the substitution of Al 3+ for Si 4+ associated with the Ca 2+ for Na + in the anorthite-albite solid solution [18,19].The excess charge at the A site (Ca 2+ for Na + ) is compensated in the tetrahedral framework (Al 3+ for Si 4+ ).One can imagine, however, for incorporations that do not form a complete solid solution (e.g., Y 3+ incorporation into zirconia, ZrO 2 [20]), the charge-balance mechanism is more complex.In addition, coupled-substitutions at other lattice positions are not always possible, depending on the mineral structure.
In the absence of secondary cation sites (e.g., SiO 4 tetrahedra) for coupled substitutions, interstitial impurities or oxygen vacancies are added to balance the total charge.The addition of an H + atom is commonly used as a theoretical charge-balancing mechanism owing to its small size and availability.For a quantum-mechanically calculated incorporation energy, the position of the additional H + atom in the substituted structure is an important parameter to consider.For example, for NpO 2 + substitution of , three positions for the additional H + atom were considered, including bonded to the axial oxygen, peroxide oxygen, and interlayer water.The additional H + atom was found to favor the site on the axial oxygen by 0.8 to 1.3 eV with respect to the other possible sites.This H + substitution mechanism is also employed in the example below detailing NpO 2 + incorporation into sulfate and carbonate minerals (see Section 3.1).
In the introduction of this study, we highlight previous approaches from the literature for calculating incorporation energies using quantum-mechanical methods (Section 2.1).In addition, a detailed description of our new method for calculating incorporation energies is presented (Section 2.2) along with examples of its implementation (Section 3).Subsections denoted by an * describe new research presented in this study.Finally, computational challenges that are addressed in this new method and those that remain are presented (Sections 4 and 5).Our intent is to help advance methodologies implemented in low-temperature computational mineralogy such that we can more accurately and predictively model environmental mineralogical phenomena, namely incorporation.

Previous Approaches
Traditionally, optimized geometries and charge densities, rather than energies, have been evaluated for incorporated systems (e.g., [9,21,22]).Quantum-mechanically optimized geometries are compared with experimental measurements, such as extended X-ray absorption fine structure (EXAFS) measurements, to identify coordination environments and distances to nearest neighbors.In this way, impurity incorporation, as opposed to adsorption, and the impurity geometry within a host phase can be confirmed (e.g., [9]).In addition, changes in charge or electron density about an impurity can be identified along with associated changes in the electronic band structure.However, quantum-mechanical incorporation energy calculations are less common due to some inherent computational challenges.
Demichelis et al. [23] studied the incorporation of water into calcium carbonate to form hydrous carbonate minerals.Using a DFT-based approach, the authors characterized the thermodynamics of water incorporation reactions starting from low water content carbonates (calcite, aragonite) and adding water to create hydrated phases (monohydrocalcite, ikaite).A generalized carbonate hydration reaction is: While the structures of both hydrated and anhydrous carbonates have been accurately modeled by DFT methods, the thermodynamics of the hydration reactions have not been reproducible with standard DFT techniques.A range of density and hybrid functionals (e.g., Perdew-Becke-Ernzerhof (PBE), Adamo-hybrid PBE (PBE0), Becke three-parameter Lee-Yang-Parr (B3LYP) [13,[24][25][26]), as well as dispersion corrections for van der Waals interactions (e.g., B3LYP-D2, PBE-D2 [27,28]), were used to determine the most accurate computational approach.In addition, this study evaluated the treatment of water in two phases: gas and solid (ice).The computed values of enthalpy, entropy, and Gibbs free energy for the reactions were determined based on quantum-mechanical geometry and vibrational frequency calculations.By combining their calculated internal energies with experimental values for heat of fusion and heat capacity, the enthalpies of carbonate hydration reactions were significantly improved.There was considerable variability when different functionals were implemented, but dispersion corrections provided the closest match to experimental data.
Tsuchiya and Wang [29] used ab initio determination of thermodynamic properties to understand ferric iron incorporation into MgSiO 3 perovskite.In this study, the reaction energy of the incorporation of Fe 3+ into Mg-perovskite was not explicitly determined.However, the, calculation of phonon frequencies, using a quasi-harmonic approximation, allowed for derivation of properties such as Helmholtz free energy, Gibbs free energy, bulk modulus, and heat capacity of the phase with the most favorable incorporated geometry.For these calculations, the authors employed an LSDA + U approach (the "S" in LSDA stands for the treatment of unpaired spins, whereas the "U" denotes the application of a Hubbard U parameter as an approximation for the localization of valence orbitals, [30]) combined with a direct LDA method [31] on an 80-atom supercell of Fe-incorporated MgSiO 3 .
Another study by Kuo et al. [12] investigated solid phase reactions of technetium and ruthenium incorporation into rutile (TiO 2 ).This study examined the possible dopant configurations within the rutile host and the extent to which incorporation may occur as a function of temperature.Using two DFT-based codes, local Density functional calculations on Molecules (DMol 3 ) [32,33] and Vienna Ab initio Simulation Package (VASP) [34][35][36][37], the authors assessed a variety of clustering schemes for Tc and Ru within a large rutile supercell.The energetics of binary dopant clustering in adjacent and near neighbor configurations were compared to single defect incorporation.Overall, Tc single defects were shown to have lower energy than Ru ones, indicating that Tc would more readily enter the rutile structure.The extent of dopant incorporation and preferred cluster configuration was modeled over a range of temperatures and the computed results correlate well with experimental observations [38][39][40].
Technetium was also the subject of a study by Skomurski et al. [8] who examined the potential incorporation of Tc 4+ and TcO 4 − into hematite.Incorporation of these two Tc species was modeled by coupled substitution of Tc 4+ and Fe 2+ for two ferric iron atoms and by TcO 4 − occupying a lattice vacancy.Energies were determined by using unrestricted Hartree-Fock calculations in Crystal06 [41].
Incorporation energies in this study consider only calculated electronic energies and do not include any thermochemical energy contributions.Ionic species in the reaction are calculated as charged clusters in the gas phase with no hydration component.Incorporated Tc 4+ and Fe 2+ geometries were determined over a range of Tc 4+ -Fe 2+ distances.The computed energies are lower at shorter distances, which is expected given a lower Coulomb repulsion between Tc 4+ and Fe 2+ than between two ferric irons.Overall, Tc 4+ incorporation was shown to be energetically feasible, with ∆E incorp on the order of −0.43 eV per Tc 4+ ion for up to four Tc within the 4 × 1 × 1 hematite supercell.In contrast, TcO 4 − incorporation into hematite is unlikely as it was shown to be energetically unfavorable and requiring a defect in the hematite lattice.These types of incorporation reaction are good candidates for using a combination of hydrated cluster and solid phases because the sources of Tc and sink of Fe is the surrounding aqueous environment.This computational treatment may provide more accurate information about the likelihood of this particular Tc immobilization pathway.Contaminant incorporation has also been explored in relation to clay minerals.In the case of clays, incorporation typically occurs via cation exchange, where interlayer cations are replaced when in contact with aqueous solutions.This type of reaction has been modeled in several studies including one by Rosso et al. [11].They simulated the replacement of interlayer K + with Cs + in muscovite, KAl 2 (AlSi 3 )O 10 (OH) 2 .Muscovite has a 2:1 tetrahedral-octahedral (T-O-T) layering ratio.In the ideal structure, the tetrahedral sites host Si 4+ , while the octahedral sites host Al 3+ .However, these sites can become disordered and lead to a net charge on the layer surface, commonly called layer charge (LC).This negative surface charge is then balanced by interlayer cations, like K + in the case of muscovite.
To examine the ability of muscovite to accept Cs into its structure, the authors set up the following reaction to explore using computational methods: X(K) + Cs 0 (g) → X(Cs) + K 0 (g) , where X represents the two-formula unit (Z = 2) primitive unit-cell of muscovite with one of the two potential K interlayer sites being considered for exchange.Geometry optimizations and total energy calculations for these species were performed using the planewave-based DFT code Cambridge Serial Total Energy Package (CASTEP) operating with a GGA/PBE setup.Ultrasoft pseudopotentials were used to lessen the computational expense.Calculation of the gas-phase K and Cs species was conducted by placing the lone atoms in a periodic box.The box size was varied to ensure that energy obtained was unaffected by the periodicity.The energetics of the reactions were considered using gas-phase K and Cs, but the authors also applied experimentally-derived hydration energies and ionization potential from the literature.This change to K + (aq) and Cs + (aq) lowered the energy of exchange by ~3 kJ/mol for all reactions.Exchanging Cs for K was explored for four different scenarios designed to separately evaluate the influence of LC, cation radii/interlayer geometry, and Al/Si substitution on the exchange energetics.Overall, they found that increasing LC promotes Cs exchange for K.The authors also discovered that K and Cs have opposing effects on the interlayer distance when the LC is high (K causes the distance to shrink, while Cs allows the layers to spread apart).The most favorable reaction, with an exchange energy of −5.3 kJ/mol, was achieved using aqueous cations and muscovite layers with no Al/Si substitution creating a net LC of −1.
These examples indicate that a comprehensive method is required to treat incorporation from an aqueous phase that allows for either a complete computational treatment of sub-equations or for usage of, e.g., experimental hydration energies, where available.

Method/Research for Calculating Incorporation Reactions Using Aqueous Source and Sink Phases
Quantum-mechanical calculations of reaction energies require calculations of each reactant and product using the same computational theories and parameters.Thus, as seen in previous studies, all reactants and products are calculated using either periodic boundary conditions or cluster models, but not a combination of the two.However, the reactant and product describing the source and sink phases for the substitution couple are typically charged aqueous species.Thus, calculations that are more focused on the actual geochemical process of incorporation at the mineral-water interface should combine periodic and cluster calculations such that the solid host phase and aqueous source and sink phases are each modeled with the best accuracy.
We have developed a method for evaluating incorporation of aqueous species using a combination of periodic and cluster calculations.One governing rule for this method is that all calculations within a single reaction are performed with the same computational method.Scheme 1 outlines the computational steps taken in our first study using charged aqueous source and sink phases for incorporation reactions.Reaction (1) describes the traditional reaction, where the sources for Np 5+ and H + are Np 2 O 5 and water (i.e., an H 2 O molecule calculated in a 10 × 10 × 10 Å 3 unit cell) and the sink for U 6+ is UO 3 .Reaction (2) describes the reaction using neutral source and sink molecules calculated using periodic boundary conditions in CASTEP [42].While the overall reaction energy is negative, indicating favorable incorporation, the source and sink molecules are environmentally uncommon, even in carbonate-free water.Therefore, an extra step is taken to consider ionization and hydration of the neptunyl and uranyl (Reaction (3)).The ionization and hydration step is performed using Conductor-like Screening Model (COSMO) calculations as implemented in DMol 3 [43,32,33].The energy contributions from the aqueous species in Reaction (3) are a combination of the total energy for the gaseous molecule and the hydration energy based on the COSMO calculation.Thus, the final combined reaction shows the cancelation of the periodic and cluster calculation for the gas phase molecule, as well as the hydration of that molecule.An example where ionization and hydration are treated separately and the influence of both of these can be evaluated is shown in Scheme 2 and, in more detail, in Appendix Table A1.Scheme 1. Equations describing the method for calculating incorporation energies using charged aqueous species, where Reactions (1) and ( 2) are performed on neutral species using CASTEP and Reaction ( 3) is performed on neutral and charged aqueous molecules in DMol 3 .In the reactions below, K 2 (UO 2 ) 2 (SiO 4 ) 2 (H 2 O) 3(s) stands for boltwoodite, whereas K 2 (UO 2 )(NpOOH)(SiO 4 ) 2 (H 2 O) 3(s) stands for the neptunyl-incorporated boltwoodite.The combined reaction is Reaction (2) minus Reaction (3) [1].
In our more recent investigations ( [7], this study), this methodology has been further developed such that five primary reactions are calculated individually and combined for the total reaction energy (Scheme 2 modified from [7]).The series of reactions in Scheme 2 starts with using the mineral host (PbSO 4 ) and UO 3 as a solid host for UO 2 2+ and PbO as a solid sink for Pb 2+ (Reaction ( 1)).This setup of the equation with only solid phases on both sides of the equation describes the relative stability of these solid phases and can be used for further solid solution calculations.The example given in Reaction (1) of Scheme 2 would be equivalent to going from the Pb end member to a quarter fraction of the uranyl end member.Theoretically, the incorporation energy using such a periodic approach should be independent of the actual program used (e.g., CASTEP, DMol 3 , Crystal) as long as suitable computational parameters (such as, basis set, density functional, k-point density) are used.Our calculations indicate, however, that small differences are observed depending on the program used, on the order of a few tenth of an eV.In order to make the transition from solid to aqueous sources and sinks, solid phase and cluster calculations must be added together to get to an overall reaction of the periodic solid phase with aqueous contaminants.The first step (Reaction (1)) of this process describes incorporation using the traditional periodic solid oxide source and sink phase approach.All calculations for this reaction are performed using periodic boundary conditions.In addition, any H 2 O, added as a source for hydrogen, is calculated as gaseous H 2 O with periodic boundary conditions.Literature values were used for the condensation energy of water.
Scheme 2. Reactions describing the method for calculating incorporation energies using charged aqueous species.Reactions ( 1) and ( 2) are performed on neutral species using periodic boundary conditions in CASTEP and DMol 3 .Reaction ( 3) is performed on neutral molecules with and without periodic boundary conditions in DMol 3 .Reactions (4) and 5) are performed using cluster calculations in DMol 3 .Hydration energy is implicitly included using a bulk dielectric fluid as implemented with COSMO in DMol 3 , or better, by a combination of explicit water molecules in the 1st and 2nd hydration sphere and a dielectric fluid around.The combined reaction is the sum of Reactions ( 1)-( 5) (modified from [7]).A more detailed set of reactions for neptunyl incorporation is provided in Appendix Table A1.
Reaction (2) describes the conversion of the solid oxide source and sink phases to neutral gaseous species in a periodic arrangement.That is, the neutral molecules are calculated in a large unit cell (≥10 × 10 × 10 Å 3 to minimize molecule-molecule interactions especially when they carry a dipole moment) using periodic boundary conditions.Reaction (3) describes the transition from a periodic to a neutral molecular species, where all calculations are performed using a program that can use consistent computational parameters (e.g., basis sets, density functionals, spin treatment, spin-orbit coupling if applied) for cluster and periodic approaches.The work presented here is done using DMol 3 , but other programs are available that could be used in a similar fashion (e.g., Crystal09 [41,44]).Ideally, and this can be used as a test for computational consistency, the energy of transition from periodic to non-periodic treatment (ΔE 3 in Scheme 2) should be much less than the overall reaction energy (ΔE in Scheme 2).The magnitude of ΔE 3 is due to any dipole-dipole interaction between molecules in periodic unit cells.For most reactions in this study on actinyl incorporation, these neutral molecules are metal-hydroxide complexes (e.g., NpO 2 (OH) 2 ).Reaction (4) describes the dissociation and ionization of the gaseous molecular species using cluster calculations.This reaction was previously considered as the ionization reaction, but in fact, an energetic component of dissociating the ligands from the cations is inherent to the calculation.Thus, it has been redefined as the dissociation and ionization energy.Reaction (5) describes the hydration of the ionized species.In this example, hydration was considered implicitly using the COSMO bulk dielectric fluid model as implemented in DMol 3 .For even higher accuracy [45], the first and second hydration sphere are considered explicitly.Reactions (4) and ( 5) could be combined (and are in Scheme 1); however, their separation allows one to observe the impact of the hydration on the overall reaction.The sum of reaction energies (1) through (5) is the energy for the incorporation reaction based on the charged aqueous source and sink molecules.

Np-Incorporation into Uranyl Minerals
Quantum-mechanical calculations have been used to gain atomistic understanding of neptunyl incorporation into uranyl minerals, namely studtite + and NpO 2 2+ substitution for UO 2 2+ in studtite were evaluated using the traditional solid oxide source and sink phases (i.e., Np 2 O 5 and UO 3 ), resulting in incorporation energies equal to 1.12 and 0.42 eV, respectively.Note that a high spin state was used for Np 5+ (Np 5+ has two unpaired 5f electrons) in all calculations and the spin ordering was evaluated in previous studies for Np 2 O 5 [2].For NpO 2 + incorporation, an additional H + atom was added to the axial neptunyl oxygen for charge balance, and H 2 O, calculated with periodic boundary conditions, was used as the source.A comparison of source and sink phases (e.g., oxide, silicate) was performed in the study of Np-incorporation into boltwoodite, which led to the development of the methodology presented in this paper for calculating incorporation reactions using charged aqueous molecules as source and sink molecules.Scheme 1 shows a subset of those results, where the incorporation energy for reactions using the charged aqueous molecules was 0.13 eV compared to 0.97 eV for the reaction based on solid oxide source and sink phases.Boltwoodite is a uranyl sheet silicate composed of layers of uranyl hexagonal bipyramids connected by silica tetrahedra (Figure 1).Monovalent cations (K + and/or Na + ) and water are contained between the uranyl silicate sheets (in the interlayer).Due to the complexity of the boltwoodite structure, multiple charge-balancing incorporation mechanisms were evaluated, including (1) addition of a H + atom; (2) coupled substitution of an interlayer cation (divalent for monovalent); and (3) coupled substitution within the uranyl silicate sheet or intra-layer (PO 4 + for SiO 4 ).The incorporation limit was estimated for the interlayer and intra-layer coupled substitutions by calculating the enthalpy of mixing for the boltwoodite/Np-modified boltwoodite solid solution series.That is, simple solid solution calculations were able to be performed because the theoretical Np-boltwoodite could be constructed in which all of the U sites were replaced by Np.Thus, the Gibbs free energy of mixing for the solid solution was approximated based on a Margules fit for the enthalpy of mixing and the −TΔS part of ΔG was approximated by using the configurational entropy (Equation (2).The configurational entropy consists of two terms, the basic configurational entropy for a binary solid (i.e., )) and a term describing the probability of a configurational relationship between the two substituted cations.That is, for the single boltwoodite unit cell, there are two possible locations for Np and two locations for P substitution, in the case of intra-layer substitution.This entropy term can be fit based on the difference in the calculated enthalpy for the different configurations.(1 ) ( ln( ) (1 ) ln(1 ) (1 )((1 ) ln ln )) H 1 − H 2 is the difference in total energies of the host mineral with the substitution in sites 1 and 2, respectively.Figure 2 shows the approximated Gibbs free energy of mixing curves for the intra-layer and interlayer coupled substitution mechanisms.Based on our assumptions, the Gibbs free energy of mixing is always negative near the end-member compositions due to the infinite entropy contribution at that point (see inset of Figure 2).
The incorporation limit is calculated as the minimum of the Gibbs free energy of mixing.Thus, the resulting expression for the limit of incorporation is shown in Equation (3).Based on this approximation for the thermodynamics of mixing in the U-Np-boltwoodite solid solution, the approximate limit of Np incorporation into boltwoodite is 3.12 × 10 −4 molar fraction Np at 300 °C (172 ppm Np/(boltwoodite)) for the interlayer coupled substitution and 1.39 × 10 −3 molar fraction Np at 300 °C (768 ppm Np/(boltwoodite)) for the intra-layer coupled substitution.The resulting incorporation limit was in relative agreement with experimental results, which showed hundreds of ppm of Np incorporated into uranyl sheet silicates with charged interlayer cations [46].(1 2 ) ln (1 4 3 ) ln (2 3 ) 2 inc (1 ) Figure 2. Gibbs free energy of mixing (ΔG mix ) vs. composition at 300 °C for the Np-U boltwoodite solid solution series based on the intra-layer (solid) and interlayer (dashed) coupled substitution mechanisms.The inset highlights the negative ΔG mix near the end-member composition (modified from [1]).

U-Incorporation into Sulfate and Carbonate Minerals *
One study implementing this new method was conducted by Walker and Becker [7] and investigates uranyl (UO 2

2+
) and neptunyl (NpO 2 + ) incorporation into carbonate and sulfate mineral phases.Five sulfates (anglesite, anhydrite, barite, celestine, gypsum) along with five carbonates (aragonite, calcite, cerussite, strontianite, witherite) were considered as host phases for the radionuclide-bearing ions.The reactions are set up such that uranyl and neptunyl replace cations in the crystal lattice.For neptunyl incorporation, the charge imbalance created by substitution of a monovalent species for a divalent cation is compensated by a nearby H + .
Calculations to determine the energetics of incorporation reactions using solid sources of the actinide ions and solid sinks for the replaced cations were carried out using two quantum-mechanical codes: CASTEP and DMol 3 .Both codes were used to calculate geometries and energies of periodic species (solid and molecular), while DMol 3 alone was used for cluster calculations.Computational parameters were held consistent throughout the calculations and included the use of the generalized gradient approximation (GGA) with PBE functional.For CASTEP calculations, a planewave cutoff energy of 500 eV was selected (for DMol 3 calculation, a double numerical basis set with polarization d-functions (DND) was applied [32]).In addition, ultrasoft pseudopotentials were used to approximate core electrons and lessen computational expense [47].For neptunyl-bearing phases, a spin-polarized approach was used to account for the two unpaired 5f-like spins on Np 5+ .For the aqueous cluster calculations in DMol 3 , a dielectric continuum model, COSMO, was used in combination with increasing numbers of explicit water molecules added to simulate hydration of charged aqueous species.It was found that when used as the sole hydration mechanism, COSMO underestimates hydration energies for the ions relevant to this study (e.g., Ca 2+ , Ba 2+, Sr 2+ , Pb 2+ , etc.) by approximately 2-4 eV.This difference from experimental results is significant and can have a major influence on the overall ∆E of incorporation.The accuracy of computed hydration energies is improved when some degree of explicit hydration is applied in addition to COSMO.This approach gave hydration values within ~0.3 eV of experimental findings for UO 22+ and NpO 2+ Additionally, literature values were used for the condensation energy of water and the hydration energy of smaller ions, like H + .
Overall, calculated incorporation reaction energies show that gypsum and aragonite are the two most favorable hosts for uranyl and neptunyl.The incorporation energies for these reactions are as follows: ∆E gyp = 0.19 eV and ∆E arag = 0.27 eV for uranyl, ∆E gyp = 0.17 eV and ∆E arag = −0.37 eV for neptunyl.The least favorable structures for uranyl and neptunyl incorporation from these calculations are calcite, barite, and witherite.Optimized geometries for the incorporated structures were studied in detail and some trends were observed.In sulfate mineral hosts, sulfate groups rotate to accommodate uranyl and neptunyl species.Sulfate-group O atoms are moved out of alignment with the axial O on the actinyl ion.The coordination environment of the actinide species can also be observed and is composed by either four or five sulfate oxygen atoms.For carbonates, similar relationships are observed with one notable example being the uranyl-incorporated aragonite structure.In this structure, carbonate groups arrange themselves such that they achieve six-fold oxygen coordination in the equatorial plane.The stability of this structure gives aragonite the lowest (most favorable) calculated uranyl incorporation energy of 0.27 eV.
Incorporation reactions were also studied for selected systems with two types of lattice defects: vacancies and impurities.Vacancies in the carbonate and sulfate hosts were created by removing an adjacent cation-anion pair.Uranyl or neptunyl was then substituted for a cation near the vacancy.The presence of these vacancies increases the energy of the pure carbonate and sulfate phase and thus, makes incorporation energies more favorable, even negative for two of the tested uranyl incorporation reactions (celestine, cerussite).In addition to lattice vacancies, impurity ion defects were tested for a series of neptunyl incorporation reactions.For this type of reaction, a cation-anion pair in the sulfate or carbonate host was replaced by NH 4 Cl.Ammonium chloride was chosen as it is a common reagent for precipitation experiments [48] because ammonium sulfates/carbonates and chlorides of divalent cations are more soluble than the respective host minerals chosen in this study.The incorporated ammonium chloride can be substituted with an actinyl and H + ion in solution.These reactions showed an increase in incorporation energy for some minerals (anglesite), while others became slightly more favorable (anhydrite, celestine).
Orbital projections and partial density of state (PDOS) spectra can be used to analyze the electronic configuration of incorporated mineral structures.For example, Figure 3 illustrates the electronic interaction between neptunium and the carbonate O atoms in neptunyl-incorporated aragonite.The neptunyl coordination environment is shown in Figure 3a.The neptunium is coordinated by 5 O atoms in the equatorial plane, and the hydrogen bond that forms between the neptunyl O and the charge-balancing H + is 1.94 Å.In Figure 3b, the highest occupied molecular orbital (HOMO) shows the eight-lobed shape of the Np 5f orbitals, not overlapping the O 2p orbital dumbbell.The HOMO is approximately 0.12 eV below the Fermi level (E F ) and is a non-bonding orbital.The orbital shown in Figure 3c, however, is a σ-bonding orbital, which is approximately 1.08 eV below the Fermi level (E F = 2.24 eV) and shows the partial covalent character of the bond formed between the neptunium and the carbonate O atoms.Both the HOMO and the bonding orbital are located on the PDOS projection in Figure 3d (their contribution to the DOS is indicated by arrows).The contribution from the Np 5f orbitals fills in the top of the valence band just below the Fermi level.As compared with aragonite without neptunyl incorporation, these states with partial Np 5f character fill a portion of the aragonite bandgap of about 4 eV and reduce it to about 1.5 eV.

U-Incorporation into Fe-Oxide Minerals *
Aqueous uranium, commonly found as the UO 2 2+ uranyl ion, can form hydroxide and water complexes which have been shown to sorb strongly to Fe 3+ -oxide surfaces, including magnetite [49,50].Magnetite (Fe 3 O 4 ) is an interesting mineral to consider as it is a common corrosion product of steel, formed during anoxic water-steel interaction via the following reaction: As a result, it is conceivable that magnetite may be one of the first materials to potentially immobilize uranium that has been previously released from a containment vessel.
This study initially considered incorporation of both U 6+ and U 4+ via a variety of solid phase reactions (i.e., the source of U and the sink of Fe are both solid phases) using a quantum-mechanical periodic planewave code named CASTEP.A high-spin state was applied for U 4+ , Fe 2+ , and Fe 3+ .The computational parameters used were GGA/PBE as a density functional, a planewave cutoff energy of 600 eV, 0.06 Å k-point spacing, ultrasoft pseudopotentials, and a convergence tolerance of the self-consistent field (SCF) procedure of 10 −5 eV.Since both U 6+ and U 4+ have a higher positive charge than the replaced Fe 2+ or Fe 3+ , charge-balancing of the incorporated oxide phase was achieved by a structural vacancy of either Fe 2+ or Fe 3+ near the incorporation site (Figure 4).By comparing the energetics of vacancies in the tetrahedral and octahedral site, vacancies were found to be more energetically favorable in the octahedral lattice position.The energies of these reactions are comparable to previous studies involving uranium incorporation into hematite, and may represent vacancy formation processes as they are likely to occur in nature [51].Within the tested configurations: U 6+ in octahedral site with octahedral Fe 3+ vacancy, U 6+ in octahedral site with tetrahedral Fe 3+ vacancy, U 4+ in octahedral site with octahedral Fe 2+ vacancy, U 4+ in octahedral Fe 3+ site along with reduction of neighboring octahedral Fe 3+ to Fe 2+ to balance charge; U 6+ incorporation into the octahedral site of the magnetite structure balanced by a octahedral ferric iron vacancy was the most favorable reaction.
while these reactions provide a first-order understanding of the immobilization pathways for uranium in magnetite, a more involved approach combining periodic solid and cluster calculations can better model the true environmental interaction of aqueous uranium and iron oxides.We present calculations using this new method for one of these incorporation reactions in which U 6+ is incorporated into the ferric iron octahedral site of magnetite, balanced by an octahedral Fe 3+ vacancy.For these calculations, we have used DMol 3 , an ab initio DFT-based code, to model incorporation into a 14-atom primitive cell of inverse-spinel magnetite.Iron cations in magnetite were all high spin and set in a ferrimagnetic arrangement.Tetrahedral sites were spin down, while all octahedral Fe atoms were spin up, giving a net moment of 8 for the primitive cell that has a formula of Fe 6 O 8 [52].Calculations were run using the generalized gradient approximation (GGA) in conjunction with the PBE functional.DFT semi-core pseudopotentials approximated the interactions of core electrons with valence shells.For oxygen, iron, and uranium, 6, 16, and 14 electrons, respectively, were calculated explicitly.An energy convergence tolerance of 2.7 × 10 −4 eV was selected for both SCF and geometry optimization cycles.For our incorporation reactions, DMol 3 total energies were used to represent overall ΔE incorp .One of the most critical terms in the overall thermodynamics for incorporation processes from the aqueous to the solid phase is the hydration energy of the aqueous species.While there are different ways to calculate these (combinations of explicit water molecules and homogeneous dielectric fluids surrounding them, static vs. molecular-dynamics approaches, force-field vs. quantum-mechanical methods), it is highly advisable to compare these calculated values with experimental results.However, an additional complication is that some experimental results are reported as ΔG values of hydration/solvation, others as ΔH.In the example described in Scheme 3, for some hydration energies, literature values from experiment have been used.For example, experimental values from the literature were used to describe the hydration energies of UO 2 2+ , H 2 O, and H 3 O + [53,54].The values for the uranyl ion and H 3 O + are ∆G solvation values, while the water value is the enthalpy of vaporization.Scheme 3. Reactions describing U incorporation into magnetite using charged aqueous species.Reactions (1) and ( 2) are performed on neutral species using periodic boundary conditions in DMol 3 .Reaction ( 3) is performed on neutral molecules with and without periodic boundary conditions in DMol 3 .Reactions ( 4) and ( 5) are performed using cluster calculations in DMol 3 .Hydration energy is implicitly included using a bulk dielectric fluid as implemented with COSMO in DMol 3 , or better, by a combination of explicit water molecules in the 1st and 2nd hydration sphere and a dielectric fluid around.The combined reaction is the sum of Reactions ( 1)-( 5).Periodic molecules and non-periodic clusters are labeled (pbc m) and (clus) respectively.A more detailed breakdown of the individual reactions is provided in Appendix Table A2.
The overall reaction of periodic magnetite with aqueous uranyl is achieved by moving systematically from solid source and sink phases, to neutral clusters and then to hydrated, charged clusters.In this reaction, in the presence of hydrated uranyl, water, and hydronium ion, magnetite is able to incorporate uranium into its structure and release hydrated ferric iron (∆E incorp = 2.84 eV).The multiple steps required to get to this overall solid-aqueous reaction are shown in Scheme 3. It is important to note that this energy cannot be considered the ∆G of the overall reaction, which would be necessary in order to convert the reaction energy into an equilibrium constant, because not all aspects of the entropy change are included.The entropy of the water molecules surrounding an ion in solution decreases during hydration and part of this process is considered in the treatment of hydration energies, while the treatment of changes in the vibrational entropies of all species involved is extremely computationally expensive.Work by Walker and Becker [7], which has been summarized in the previous section, has shown that the −TΔS contribution to the ΔG change of incorporation reaction is on the order of 0.1 to 0.3 eV (typically the entropy of the solids goes up during incorporation because the perfect order of the host mineral phase is disturbed by the incorporation process).We can assume for our solid phases containing lattice vacancies that the variations may be slightly larger, but overall, the energy contribution may be small but significant if the conversion to equilibrium constants is made.

Challenges in Validating Calculated Results against Experiments
One particular challenge is the comparison between calculation and experiment.Ideally, one would like to compare the thermodynamics of incorporation, a closely related property which is the thermodynamic limit of incorporation at equilibrium, the structure and its distortion of the lattice around the incorporation site, and, e.g., potential changes to the electronic structure ("doping") as a result of incorporation.However, for most systems, hardly any thermodynamic data are available from experiment for incorporation reactions, in particular at the ppm (and lower) scale.For some systems, data are available for incorporation limits, either from experiment or from natural samples (see, e.g., Table 1 in [55] for actinide incorporation into zircon, or U incorporation into natural garnets [4]).Incorporation limits from experiment, e.g., from co-precipitation experiments, often suffer from the fact that equilibrium has not been reached during the time of incorporation.Another potential difficulty in co-precipitation experiments is that it is not always trivial to distinguish structural incorporation of single ions from incorporation of a different, e.g., nanoparticulate phase.While nature typically gives "samples" more time to reach equilibrium, even natural materials may not be fully equilibrated, especially if no "lubricating" medium, such as water, is available to help with the transport of ions that have been incorporated in excess.On the lower end of incorporation concentration, i.e., below the thermodynamic limit, it is possible that nature did not provide enough of the species to be incorporated to reach the highest possible incorporation concentration.Thus, in many cases, computational results need to be compared with structural data from, e.g., EXAFS experiments.While a good comparison between properties such as coordination number and coordination distances are a necessary requirement to judge if the calculations are on the right track, they are by no means sufficient to evaluate changes in the thermodynamics upon incorporation.
For incorporation of foreign species, interesting and measurable changes of electronic properties may occur.For example, incorporation may lower the width of the bandgap, and it may also change the character of the semiconducting host, i.e., if it is a p-or n-type semiconductor [3][4][5][6].If the computations compare well with the experiment in terms of the electronic parameters that can be verified by experiments, such as conductivity, Hall effect, or inverse ultraviolet photoelectron spectroscopy, then, the calculations can offer significantly more insight into electronic property changes resulting from incorporation, such as changes in the band structure, the density of states, redistribution of the electron and spin density, the shape of individual bonding, non-bonding, and anti-bonding orbitals between the incorporated species and the host phase.Incorporation at higher concentrations may also generate specific vibrational modes of the incorporated species in the host that are different from host modes or of the species in solution.If these are calculated, they can be compared with Raman or IR spectroscopy results.
Since thermodynamic data, e.g., as obtainable from calorimetric measurements may be hard or impossible to determine, in order to compare experimental and theoretical results as comprehensibly as possible and to test the validity of the approach described in this study, experiments are most promising on systems that allow for incorporation at the percent scale, reach equilibrium within a suitable amount of time, and structure, electronic, and vibrational properties should be obtained.The combination of these would allow for the highest degree of validation.
Incorporation of uranium into iron (hydr)oxides has been investigated experimentally (e.g., [10,49,56]) and the results of these studies can help us validate or improve our computational approach.Unfortunately, these studies have yet to tackle the determination of incorporation energetics, a key focus of our quantum-mechanical calculations.At this point, these experimental studies can confirm certain minerals as possible hosts for contaminant species and also provide comparison for the incorporated geometries.For example, Nico et al. [10] studied the incorporation of U 6+ into goethite and magnetite.Their batch chemistry experiments showed evidence for the precipitation of solid uranyl phases (hydroxides and carbonates), but also for incorporation (and sorption) of U into the structures of iron (hydr)oxides.
The incorporation of oxidized uranium into iron (hydr)oxide phases was confirmed by analyzing the mineral substrates after a thorough extraction process to remove adsorbed U.The incorporated geometries were characterized using X-ray Absorption Near-Edge Spectroscopy (XANES) and EXAFS analyses.The fitting of XANES and EXAFS is complex; for example, it is not trivial to distinguish the exact percentages of, e.g., aqueous vs. adsorbed vs. incorporated U in goethite or magnetite [10].The derived bond distances and coordination numbers in the three simultaneously-occurring phases come with a certain amount of error, such that any comparison to computationally derived geometries must be mindful of potential discrepancies in bond distances.The U-incorporated magnetite structures calculated in this paper with U being in an octahedral Fe site have U-Fe distances with four Fe surrounding a U site (for charge compensation, two Fe in the coordination sphere were removed) are 3.08 Å compared to the "ideal" structure in [10] listed as 2.97 Å which would mean a difference of about 0.1 Å.For U-O distances in the first coordination sphere (Figure 4), the calculated values (2.13 Å for four equatorial O, 2.26 Å for the two axial O) compare well with "ideal" values in [10] of 2.09 Å (in the experimental paper, no distinction was made between equatorial and axial oxygen).).Uranium has replaced a Fe 3+ in an octahedral site and a second Fe 3+ octahedral site is now vacant to balance charge.Bond distances (in Angstroms) for the first U-O sphere are indicated.Colors for the atoms are as follows: U 6+ blue, Fe 2+ purple, Fe 3+ green, O 2− red.

Theoretical Challenges
While the method presented in this paper is consistent in terms of combining aqueous charged species from solution and neutral host and incorporated phases with periodic boundary conditions, it may be worth trying to analyze the potentially greatest sources of error.In quantum-mechanical calculations, such as those presented for incorporation of species into mineral phases, it is always advantageous to apply the highest possible accuracy in terms of basis sets (high energy cutoff for planewave approaches or high-quality basis set for atomistic basis sets), k-point density, suitable DFT or hybrid functional, and other corrections (e.g., relativistic and spin-orbit coupling effects, localized f orbital treatment that can be treated using the Hubbard U formalism [3][4][5][6] or computationally-expensive hybrid DFT functionals such as the Heyd-Scuseria-Ernzerhof (HSE) functional [57]).All of these computational settings control incorporation energy changes as long as incorporation is treated using solid sink and source phases, the incorporated and released ion have the same charge and form similar complexes (e.g., actinyl ions are linear dioxides while Fe ions that may be released are just single ions, either in solids or in the aqueous phase, such that the latter criterion for replacement is not fulfilled), and the incorporation mechanism for two replacements to be compared is the same.However, typically, not all of these criteria are fulfilled such that other aspects of the overall suite of equations to be considered may control the error of the overall incorporation energy.
If the source and sink phases to be considered are aqueous species, the most significant source of error tends to be the treatment of hydration energy.First, there are, in principle, a number of computational treatments of hydration energy.One way is to "bond" a certain number of water molecules to the ion of interest in a static way according to: Equation ( 6) can be evaluated using a quantum-mechanical (which is advised due to the peculiar behavior of actinide valence orbitals) or a classical force-field approach.The reaction energy in Equation ( 6) tends to become more and more negative with an increasing number of hydrating water molecules used (n) due to the static character of the calculation (at a finite temperature, water molecules further away from the central cation would have an increasingly less ordered orientation).One way to overcome the static character is using a (quantum-mechanical or classical) molecular-dynamics calculation [58][59][60].Since there are hydrogen bonds between the original water molecules (Note: This is liquid water and not water from the gas phase!), the n(H 2 O) part on the left can be corrected by adding in the heat of condensation of water (on the order of −0.4 eV per water molecule) or by simulating the water as a droplet or water in a periodic box.Similar to the approach shown in Equation ( 6), the ion can also be "dipped" in a dielectric fluid using, e.g., a COSMO, polarizable continuum model (PCM) or conductor-like PCM (C-PCM) [61][62][63][64] approach.Ideally, one would combine the first or first two hydration spheres using explicit water molecules with a dielectric-fluid model in a static, or even better, dynamic approach.In either case, it may be necessary, and highly advisable, to benchmarks the atomistic calculations against experimental values, where available, or other computational values from the literature because discrepancies of a few eV (which is a lot especially once energies are converted to equilibrium constants) are not uncommon.When comparing work from different sources, special attention has to be applied to the nature of such hydration energy values because they may be ΔG, ΔH, or values of some other aspect of the hydration/solvation energy.This often depends on the nature of the calculation (e.g., static vs. dynamic) or the way experimental hydration energy values were obtained.The question of the thermodynamic nature of the hydration energy, in particular the entropy changes involved, leads to a more general consideration about the nature of the overall thermodynamics of incorporation with all related aspects of changes in the entropy.Ideally, one would like to obtain ΔG values in order to judge, how much incorporation is thermodynamically possible.In addition, in order to obtain ΔG values for the hydration energy, other sources of entropy changes need to be evaluated.Even if solid source and sink phases are used, there may be a change in the vibrational entropies.Since pure mineral phases may have less vibrational modes than phases with incorporated species, the −TΔS term tends to decrease by a few tens of eV, making incorporation somewhat less energetically uphill and, thus, somewhat more likely.Another (sometimes major) source of entropy changes may be a change from the number of reactants to the number of products; if the number of species goes up, so does the entropy, and vice versa.
In this paper, we have presented an approach to combine charged aqueous species as source and sink phases for incorporation reactions into solid mineral phases.This approach is widely applicable to a plethora of petrological, mineralogical, materials, environmental, or nuclear science applications.While the usual thoroughness of choosing computational parameters is indicated, special care has to be applied to the treatment of hydration energies, their distinction in terms of ΔG and ΔH values, and the treatment of other sources of entropy changes (vibrational, change in the number of species on both sides of the reaction).With these issues in mind, it should also be possible to create a database that allows for a widespread usage of this modular way of combining equations.This is because a lot of the equations (e.g., hydration energies) can be "re-used" with future users being able to focus on the energetics of the host, the host with the incorporated ion, and the (easier to calculate) solid source and sink phases.With such a publically available dataset, this approach may be a powerful tool to evaluate the thermodynamics of incorporation reactions, and aid in the benchmarking of experiments that tend to be subject to kinetic limitations.

Figure 1 .
Figure 1.Atomic model showing the uranyl silicate sheets of boltwoodite K(UO 2 )SiO 3 OH(H 2 O) 1.5 , where the interlayer (between the sheets) is filled with K and H 2 O. Colors for the atoms are as follows: U 6+ blue, Si 4+ yellow, K + purple, O 2− red, H + white.

Figure 3 .
Figure 3.The coordination environment and electronic structure of neptunyl in aragonite.(a) Coordination environment of neptunyl with nearest-neighbor distances.One H + per unit cell is added to charge compensate the replacement of Ca 2+ by NpO 2+ The second H + ion comes from the adjacent unit cell (red = O, gray = C, white = H, blue = Np); (b) Wavefunction showing the non-bonding character of the HOMO (no overlap of the eight-lobed Np 5f and O 2p orbitals; (c) σ-bonding orbital (E = 1.15 eV or 1.08 eV below EF) between Np 5f and O 2p orbitals from two O atoms of the same carbonate molecule (the other three σ-bonding orbitals, not shown, have comparable electron binding energies, with one bond each from the three other CO 3 2− ; and (d) Partial density of states (PDOS) projection for the energy range from −15.0 to 10 eV.

Figure 4 .
Figure 4. Atomic model showing the primitive cell for U 6+ incorporated magnetite (Fe 4 UO 8 ).Uranium has replaced a Fe 3+ in an octahedral site and a second Fe 3+ octahedral site is now vacant to balance charge.Bond distances (in Angstroms) for the first U-O sphere are indicated.Colors for the atoms are as follows: U 6+ blue, Fe 2+ purple, Fe 3+ green, O 2− red.

Table 1 .
Examples of incorporation reactions studied for better understanding of the transport of radionuclides in the environment.