Understanding the Stability of Salt-Inclusion Phases for Nuclear Waste-forms through Volume-based Thermodynamics

Formation enthalpies and Gibbs energies of actinide and rare-earth containing SIMs with silicate and germanate frameworks are reported. Volume-based thermodynamics (VBT) techniques complemented by density functional theory (DFT) were adapted and applied to these complex structures. VBT and DFT results were in closest agreement for the smaller framework silicate structure, whereas DFT in general predicts less negative enthalpies across all SIMs, regardless of framework type. Both methods predict the rare-earth silicates to be the most stable of the comparable structures calculated, with VBT results being in good agreement with the limited experimental values available from drop solution calorimetry.

Nuclear waste sequestration, including legacy materials from weapons programs as well as spent fuel from research reactors and potential commercial fuel recycling remains an important contemporary issue. While many reprocessing techniques exist, and repository solutions have been proposed, there is still a large research focus on how to more effectively and efficiently immobilize certain problematic radionuclides, especially those which are easily volatilized or for which waste glass loading is limited. A novel approach to simultaneously capturing multiple nuclear waste products includes the use of hierarchical architectures of porous materials. The working definition of a hierarchical material is that of a structural motif contained within a larger structure or framework. A class of materials that exhibit this structural characteristic include salt inclusion materials (SIMs).
Salt-inclusion materials exhibit a hierarchical structure that consists of a covalent mixed-oxide framework which contains a void filled with simple ionic salts. While traditional SIMs are characterized by transition metal oxides interconnected with oxyanion units of groups 14 and 15 elements such as Si, Ge, As, P 1-5 more recently, uranyl [6][7][8][9] and lanthanide 10 salt-inclusion phases are being explored for nuclear waste applications due to their porous or "stuffed" nature. The framework allows for structural variability forming uranyl-based silicate, germanate, vanadate, phosphate or borate networks with various 3-D void sizes, which are filled with ionic salts that preferentially contain radionuclides. The general description of uranyl SIMs is the structural formula [A m B n X] [(UO 2 ) p (M q O r ) t ], where [(UO 2 ) p (M q O r ) t ] is the framework consisting of uranyl cations, UO 2 2+ , and M q O r units (M = network forming ion such as Si or Ge), B n X is the salt-inclusion, and A are non-salt-inclusion cations. To widen the class of materials, ion exchange in SIMs can be performed to include targeted isotopic compositions.
Preparation of the framework materials take size and charge variations into account during synthesis; however, little is known about their thermodynamic stability, including formation enthalpies or Gibbs energies. For known phases, calorimetric methods can provide a direct measure of the formation energy of the materials, however to date there is no published literature on the thermodynamic properties of SIMs. Predictive thermodynamics is an attractive technique as it can provide insight into the thermodynamic stability of novel new structures such as SIMs, as well as guide the synthesis of newly formulated materials. Volume-based thermodynamics (VBT) is a tool developed by Glasser et al. [11][12][13] which serves to estimate thermodynamic parameters of a class of newly synthesized or even hypothetical materials when experimental thermochemical data are lacking and other theoretical modeling and simulation techniques are uncertain and can be computationally prohibitive. In this work we aim to provide a library of Gibbs energy values for a set of systems that encompass a multitude of different structural frameworks and potential salt inclusions to effectively inform the sequestration of radionuclides for waste management. To our knowledge this is the first attempt to apply VBT to complex hierarchical structures such as salt-inclusion materials.

Methods
Volume based thermodynamics (VBT). The VBT method incorporates empirical relations to generate estimated quantities of the standard entropy (S°2 98.15 ), enthalpy of formation (Δ f H°2 98. 15 ) and Gibbs energy of formation (Δ f G°2 98. 15 . The method uses crystallographic information from X-ray diffraction or density measurements if the formula mass is known, to obtain the volume per formula unit (V m ). In this work the formula unit volume is calculated by dividing the volume of the unit cell V cell (from a crystallographic information file; CIF) by the number of formula units Z in the unit cell so that V m = V cell /Z. This quantity is then used in conjunction with derived thermodynamic cycles to calculate the formation energetics, as presented in the schematic of Fig. 1.
The standard entropy is calculated with Eq. 1, where the fitted constants k (J/K/mol/nm 3 ) and c (J/mol/K) are applied with the formula unit volume, with the constants varying as to whether the system is organic (liquid or solid) or ionic (hydrous or anhydrous). In this case we take the constants as fitted for anhydrous ionic salts 12 .°= A lattice potential energy is required which is calculated from Eq. 2 and is indicative of the ability of an ionic solid to form from components in the gaseous state, where the ionic strength factor I = ∑ I nz ( 1/2 ) i i i 2 is calculated from the constituents of the salt and the salt-inclusion framework and their respective charges, with n i being the number of ion types, z i their respective charge; and A the standard electrostatic Madelung constant (121.39 nm kJ/mol) 11,12 The lattice energy is then converted into a useable enthalpic value by a multiplicative RT term that includes information on the ion types (s i ) and a constant (c i ) related to whether the ion types are monoatomic, polyatomic (linear or non-linear) as shown in Eq. 3, with R being the ideal gas constant and T the temperature in Kelvin.
The Born-Haber-Fajans cycle, which applies Hess' law is then used to calculate the standard enthalpy of formation in which the constituents of the salt-inclusion material are broken down into their gaseous ionic counterparts, where the salt inclusion components are broken down into their elemental state, and the framework consists of constituents in various oxide forms. Information regarding the gaseous components from the solid phase are obtained from auxiliary information in Table 1 and include enthalpies of sublimation or dissociation, combined with ionization potentials (IP) or electron affinities (EA) for cationic and anionic species respectively, which are found in the literature. The summation of these energies in the gas state along with the lattice enthalpy (Eq. 4) results in a value for the standard enthalpy of formation. The latter value then allows for the calculation of the Gibbs energy of formation by applying auxiliary information for the standard entropy to Eq. 1.  A mixing entropy accounts for the combining of the different components of the salt, where contributions of partially occupied and mixed salts are naturally greater than those with a single cation type. The relation is seen in Eq. 5, where n is the total number of moles and x i is the mole fraction of each constituent.
The VBT approach was applied to three different classes of salt-inclusion frameworks: Uranyl silicates (9 compounds) uranyl germanates (13 compounds) and rare-earth silicates (2 compounds). The compositions were obtained from the literature or synthesized by the methods described in 9,14 , and are listed in Table 2 along with V m values derived from available crystallographic information.
Density functional theory (DFT). The DFT calculations were performed using the code VASP 15,16 , with the Perdew-Burke-Ernzerhof (PBE) generalized-gradient approximation 17 , employing the projector augmented plane wave (PAW) method 18,19 . For calculating the enthalpies of formation of Si 2 O 5 2− and Ge 2 O 5 2− we considered the systems to be composed of a 2D sheet formed by two SiO 4 and GeO 4 tetrahedra, with three corner sharing O atoms and a −2e charge. Considering that the U atoms are surrounded by O atoms, we chose a value of U eff = 4.0 eV, which is a U eff value that is close to that obtained from experimental studies for UO 2 20,21 and has been proven to well-reproduce the structural parameters and band gaps of for UO 3 polymorphs 22-24 . The calculations were performed using 12 × 12 × 1 k-point mesh, 520 eV cutoff energy for the planewave basis set, and 10 −8 eV and 0.001 Å/eV energy and forces convergence criteria, respectively, allowing the systems to fully relax (volume, cell shape and ionic positions). For the SIMs the calculations utilized a 500 eV planewave energy cutoff, 10 −6 energy convergence criteria, k-point mesh with 3000 KPPRA (k-point density per reciprocal atom), and fully relaxed systems.
Thermochemical cycles. Each of the SIMs frameworks are broken down into individual constituents based on the available auxiliary information, where silicate and germanate oxide constituents are initially limited to SiO 2 /GeO 2 and SiO/GeO components with a single negative charge. To obtain a better representation of the silicate SiO 4 , and germanate GeO 4 tetrahedra, which often arrange in Si 4 O 10 and Ge 4 O 10 columns, the components Si 2 O 5 2− (g) and Ge 2 O 5 2− are needed and thus density functional theory (DFT) calculations were performed to calculate the formation enthalpy of these constituents for which no information is available. The anion frameworks are charge-balanced by varying the oxidation state of uranium in the uranyl cations so that the overall salt-framework is neutral. An example of a balanced Born-Haber-Fajans cycle used to calculate the Δ f H°2 98.15 is depicted in Fig. 2. The remaining constituents that make up the various silicate, germanate and rare-earth framework cycles are reported in Table 3, where the single-ion values that make up the salt-inclusions are directly taken from the auxiliary data table.

Results
The lattice potentials calculated using Eq. 2 are plotted as a function of the formula unit volume for the available SIMs in Fig. 3. The uranyl silicate materials include more versatile framework structures, where different charged frameworks and salts are considered. Both the lanthanoid (Ln) silicates and uranyl germanates (except for one structure) have the exact same framework composition. The increased variance of the salt inclusions, including their charge and composition, allows for a range of differently charged uranyl-silicate frameworks, which dictates the lattice stability, which is largely dependent on the ionic strength factor. Conversely, the germanate SIMs have identical frameworks for twelve of the thirteen structures. For both silicates and germanates with self-same frameworks, the lattice potential decreases with increasing V m , as it is inversely proportional to its cube root of the value (see Eq. 2) and the ionic strength factor is less influential due to the similarity of the salt-inclusions. The Δ f H°2 98. 15 are calculated using the auxiliary information in Table 1 and are compared with experiment and  values calculated by DFT in Table 4. Only salt inclusions which did not have partial occupancies were computed by DFT as the significantly larger unit cell required for considering partial occupancies made the calculations prohibitively computationally intensive. The Δ f H°2 98.15 value was also calculated with VBT using volumes derived from DFT relaxed structures, the energies are compared in Fig. 4. The VBT Δ f H°2 98.15 values plus the standard entropy calculated from Eq. 1 provide the Gibbs energy of formation, both of which are listed in Table 4 and the latter depicted in Fig. 5. The energies include the mixing entropy of the salt-components as noted above and as was demonstrated in Juillerat et al. 25 for alkali metals.

Discussion
The results in Table 4  This indicates that for the silicates, the charge on the framework (which contributes to a higher ionic strength factor) and the overall size of the system (such as the total number of atoms per formula unit), influences the thermodynamic stability. More negatively charged frameworks that allow for larger salt inclusions have a more negative enthalpy of formation. With equivalently charged frameworks, the silicon-rich system is found to be more stable than its uranium-rich counterpart, according to both DFT and VBT. The VBT values for [(UO 2 ) 3 (Si 2 O 7 ) 2 ] 6− and [(UO 2 ) 2 (Si 6 O 17 )] 6− framework types with identical Cs 2 Cs 5 F salt inclusion, imply that the silicon-rich composition is more thermodynamically stable (has a more negative formation enthalpy). While the increased negative value in formation enthalpy (+4.3%) might be attributed to the increase in V m (+3.8%) for the silicon-rich framework, it seems more likely that the choice of constituents for the utilized thermodynamic cycle are more influential. In the case of the silicon rich [(UO 2 ) 2 (Si 6 O 17 )] 6− framework the cycle includes the use of + UO 2 2 , which has a greater impact on the formation energetics, since both the first and second ionization potentials are included. The silicon rich framework allows for a better representation of the structure by including both  Table 3. Thermochemical cycles for SIMs framework components. Si and U in their proper Si 4+ and U 6+ oxidation states respectively. This work attempts to use U VI ions in the thermodynamic cycles whenever possible as it is a more realistic description of the system, since all frameworks but one contain this oxidation state of the uranyl cation. Nevertheless, given the limitations of the auxiliary information, + UO 2 2 is not always represented as such in the VBT cycles. As indicated in Table 3, in order to properly charge-balance the system, only a singly charged uranyl cation (UO 2 + ) is often used.  Table 4. Enthalpies of formation (kJ/mol), Gibbs energies of formation (kJ/mol) and standard entropies (J/mol/K) of SIMs from VBT compared with DFT and Experiment. *Monoclinic, ¥ orthorhombic, § hexagonal (distinctions are made for germanates of equal charged frameworks). Equivalent frameworks in both composition and charge differ only in the salt-inclusion which dictates V m , where the Cs 2 Cs 5 F salt-inclusion results in a much larger V m (15.4%) compared to the other four NaK 6 F, KK 6 Cl, Na 0.9 Rb 6.1 F, K 3 Cs 4 F salt compositions. The average framework V m was calculated as 508.0 ± 23.3 Å 3 , where the thermochemical radii of the alkali metals and halides are used to compute the V m of the salt-inclusions. The volume of the salt is then subtracted from the overall formula unit volume of the five identical framework materials, which are then averaged. The larger formula unit volume of the pure cesium containing (Cs 2 Cs 5 F) SIM leads to a formation enthalpy that is less negative than its four counterparts; a similar trend was found in 25 , where larger alkali inclusions (and therefore V m values) resulted in less negative formation enthalpies. For the remaining SIMs of the [(UO 2 ) 3 (Si 2 O 7 ) 2 ] 6− family, both DFT and VBT predict that the chlorine containing KK 6 Cl salt is the least stable structure and the NaK 6 F salt-inclusion is the second most stable structure. DFT predicts the NaRb 6 F to have the most negative formation enthalpy, whereas VBT predicts the mixed K 3 Cs 4 F salt to be the most stable. A similar result was obtained for the mixed KK 1.8 Cs 4.2 F salt in the monoclinic germanate framework presented below.
The uranyl germanate framework, [(UO 2 ) 3 (Ge 2 O 7 ) 2 ] 6− , is analogous to the silicate framework, [(UO 2 ) 3 (Si 2 O 7 ) 2 ] 6− , and twelve different salt inclusions have been incorporated into this framework producing structures in either the orthorhombic or monoclinic setting. A lone hexagonal structure with a different framework, [(UO 2 ) 3 O 3 (Ge 2 O 7 )] 6− has also been synthesized (the experimental results of all uranyl germanate SIMs are detailed in 14 ). The enthalpies of formation of the DFT and VBT values are listed in Table 4 and overall are less negative than those for the silicates with a similar framework composition. DFT values predict the average formation enthalpies of the [(UO 2 ) 3 (Si 2 O 7 ) 2 ] 6− silicates (−9365 kJ/mol) to be more negative by 16.1% than the [(UO 2 ) 3 (Ge 2 O 7 ) 2 ] 6− germanates (−7967 kJ/mol), whereas VBT predicts a difference of 5.6% between the silicates (−14781 kJ/mol) and germanates (−13972 kJ/mol). Yet the effects of the choice of constituents for the thermochemical cycles, i.e., using GeO − /GeO 2 − and SiO − /SiO 2 − reveals that large discrepancies can arise. This highlights the importance and limitations of the auxiliary information when calculating the thermodynamic cycles, especially the need to charge balance the framework components.
VBT predicts the orthorhombic structures in general to be slightly more stable than monoclinic structures. This could in part be due to the symmetry of the structures (i.e., orthorhombic crystal systems have higher symmetry than monoclinic) or the difference in the salt-inclusions. All of the systems with monoclinic symmetry consist of dihalide salts (except for the Cs 2 Cs 5 F) and are cesium rich, whereas the orthorhombic structures generally incorporate less cesium and exclusively include only single halide salts. For the monoclinic structures calculated by DFT, the trends in relative stability are in agreement with the results from VBT, such that the silver containing structure is the least stable, followed by the pure cesium compound. As with the silicates, VBT predicts the K-Cs salt to be the most stable composition, where the salt-inclusion consists of Cs 6 K 2 Cl 2 in the monoclinic form and KK 1.8 Cs 4.2 F in the orthorhombic form. DFT also predicts the monoclinic Cs 6 K 2 Cl 2 salt-inclusion germanate structure to be the most stable. VBT suggests that the increase in silver content leads to less stable structures, as the formation energetics of silver ions is much larger than that of any of the alkali metals. For the orthorhombic structures calculated using both VBT and DFT (which include the following salt structures: KK 6 Cl, KK 6 Br 0.5 F 0.5 and Na 0.9 Rb 6.1 F), the Na 0.9 Rb 6.1 F structure was found to be most stable by both VBT and DFT, where DFT treated the salt-inclusion as fully occupied Na 1 Rb 6 F. The remaining two structures are comparable, differing only in the variation of the halide (KK 6 Cl vs KK 6 Br 0.5 F 0.5 ) with the mixed Br-F halide calculated to be more stable by DFT, which is the reverse for the VBT results, although both methods each predict very similar energies. Note that the DFT calculations for partial/mixed occupancies can be problematic as they demand significantly larger unit cells which are prohibitively computationally expensive, since both structure types include salts that have partial occupancies, only half of both the monoclinic and orthorhombic SIMs could be treated with DFT. The hexagonal structure with lower germanate content is predicted to have the least negative formation enthalpy of the germanate compounds, indicating that the uranium rich composition is significantly less stable than the other synthesized framework compositions. This is analogous to the uranyl silicate results, where Si-rich (or U-poor) frameworks are more stable than the uranium rich compositions for frameworks of identical charge.
With respect to the rare earth SIMs, experimental information regarding the formation enthalpies of one of the Ln-silicate structures, [K 2 K 7 F 2 ][Gd 3 Si 12 O 32 ], is reported 26 . The VBT Δ f H°2 98.15 value from the elements for the SIM is in good agreement with that obtained by drop solution calorimetry (Table 4). Both VBT and DFT predict that the Eu-containing silicate is more stable than its Gd-analogue, however DFT under-predicts the values compared to experiment (for Gd-SIM) and VBT. For all of the SIMs considered here, VBT generally predicts more negative enthalpies of formation compared to DFT, however general trends are in agreement for the silicates, germanates and Ln-silicates.
Note that the formation enthalpies were calculated using V m , from DFT relaxed structures, if experimental data on the crystal structure is lacking. A comparison in the the results from using V m , values from experimental and DFT computed structures for the uranyl silicate and germanate systems is found in Fig. 4. Overall, the volumes calculated with DFT lead to minor differences in the VBT computed energies, with a variation of no more than 2 percent. Most of the values computed with DFT-determined volumes are more negative (more stable) than those computed with experimental values for both silicates and germanate structures.
The formation enthalpies for DFT are calculated in vacuum at 0 K, however to include temperature dependence and entropic contributions are out of the scope of this work as they are too computationally demanding and not every salt-inclusion can be treated since partial occupancies pose a problem when generating the structures. VBT does however, produce entropic values that allow calculating the Gibbs energy of formation of each of the respective compositions (Table 4 and Fig. 5). The trends for the Gibbs energies remain consistent with those calculated for the formation enthalpies in that the silicates are found to be the more thermodynamically stable structures, except for one composition, which has a much smaller V m and salt-inclusion compared to the rest of the structures considered. Similarly, more negatively charged framework structures have increased stability, where the impact of the overall charge of salt-inclusion influences this stability, i.e., more ions within the salt-inclusion increase the ionic strength factor, which contributes to the lattice potential used for these calculations.

Conclusions
In this work we compute relative stabilities of complex hierarchical structures for waste sequestration using computationally inexpensive techniques that rely on sound thermodynamic correlations. The enthalpies and Gibbs energies of formation of 24 SIMs were calculated using VBT methods and compared to the enthalpies of formation from DFT and experimental results when available. VBT and DFT results were in closest agreement for the smaller framework silicate structure, whereas DFT in general predicts less negative enthalpies across all SIMs, regardless of framework type. The uranyl-germanate structures were found to be slightly less thermodynamically stable than their silicate analogues. Both methods predict the Ln-silicates to be the most stable of the comparable structures calculated, with VBT results being in good agreement with an available experimental value from drop solution calorimetry. Additionally, DFT was used to calculate some of the framework components used in the thermochemical cycles for the volume-based methods. This allowed for a more physical representation of the structural units seen in experiment. As auxiliary information on SiO/SiO 2 and GeO/GeO 2 building blocks are limited to singly charged species, DFT aids in obtaining information on higher oxidation states, which are necessary to charge balance these complex systems. While certain thermochemical cycles yield VBT values in better agreement with DFT results, discrepancies still exist between the absolute values of both methods. Similarly, implementation of U eff in DFT, as is standard for f-electron systems, leads to lower (more negative) formation energies, however this does not resolve the disparity as the values calculated with U eff = 4.0 eV are only about ~100 kJ/mol lower than those computed using U eff = 0 eV. Improvements in the thermochemical cycles of VBT and manipulation of the U eff values might produce better agreement.