Acidity constants and redox potentials of uranyl ions in hydrothermal solutions

We report a first principles molecular dynamics (FPMD) study of the structures, acidity constants (pKa) and redox potentials (E) of uranyl (UO2 ) from ambient conditions to 573 K. It is found that UO2 2+ keeps five coordination up to 573 K whereas UO2 + transforms from 5 to 4-coordinate as temperature increases to 573 K. The FPMD-based vertical energy gap method is used to derive pKas and E s. The method is validated by comparing with available experimental data (for E under the ambient conditions and for pKas from ambient conditions to 367 K), with an uncertainty of 1–2 pKa units and 0.2 V for pKa and E . The encouraging results demonstrate that the method may be used to predict the pH–Eh diagrams of f-block elements under the conditions of hydrothermal solutions. The results show that the acidity constants of uranyl decrease with temperature and are lower than 3.0 when the temperature is above 473 K, indicating that hydrolytic forms are dominant for U(VI) in the near neutral pH range. The reduction potential increases with temperature, indicating that the reduced state is more significant at higher temperatures.


Introduction
Knowledge of the aqueous speciation of uranium is crucial for understanding the properties and behavior of uranium in geo-fluids, sediments, and geological disposal sites for nuclear wastes. 1 It has been accepted that changes in acid-base and redox conditions (i.e.pH and Eh) in geological environments play a key role in the transport and deposition of uranium. 2,3Protonation states of U-containing complexes change with the environmental pH, which influences complexation and fixation at mineral-water interfaces. 4The increase in pH and the decrease in Eh facilitate the reduction of U(VI) to U(IV), leading to the precipitation of solid phases.
Numerous studies have been carried out to investigate the aqueous properties of uranyl ions, 1,[5][6][7] as they are believed to be the most common uranium species in the supergene environment.Under ambient conditions, uranyl ions have 5 1st-shell H 2 O ligands, i.e.UO 2 (H 2 O) 5 2+ .The aqueous uranyl ion is a weak acid.A pK a value of 5.2 is recommended by a review 5 and 5.58 has been obtained by the latest experiment. 8The reduction potential of the UO 2 2+ /UO 2 + couple is 0.062 V. UO 2 + is not stable in water because it rapidly disproportionates to UO 2 2+ and U 4+ .
Recently, increasing evidence has shown that UO 2 + can be stabilized by mineral surfaces. 9,10ranium species often exist in the geological environments with a large variation of temperature.In deep geological repositories of nuclear wastes, the temperature can reach 200 1C due to radioactive effects. 11It is common that in crustal fluids, the temperature is several hundred degrees.However, the aqueous properties of uranyl ions have been poorly documented under the conditions of elevated temperatures and pressures. 1The hydration structures of UO 2 2+ and UO 2 + have not been reported as yet at higher T-P.The direct measurement of the acidity constant of uranyl ions is carried out only up to 367.4 K 8 and the reduction potentials have not been measured at elevated temperatures. 1,5][17][18][19][20][21][22] For example, Steele et al. calculated E 0 of UO 2 2+ /UO 2 + being 0.14 V vs. SHE (standard hydrogen electrode). 23Hay et al. predicted the 1st pK a of UO 2 2+ to be 9.6. 24However, it is hard to extend these solvent models to high P-T conditions because they are parameterized under the ambient conditions.Furthermore, since continuum models often ignore the atomic level details of the solvated structures, it is difficult, if not impossible, to predict possible structural transformations of metal complexes at elevated temperature.Density functional theory based first principles molecular dynamics (FPMD) treats the solutes and the solvents on the same level of electronic structure theory 25,26 and is a suitable tool for studying metal complexes under the T-P conditions relevant to geo-fluids, e.g.ref. 27-37.As a hard acid, uranyl ions favor complexing with hard base ligands, e.g.H 2 O, F À , and CO 3 2À .It has been accepted that uranyl ion-carbonate complexes predominate in near-neutral solutions at equilibrium with the atmosphere. 76][47][48][49] Tests on molecular acids and metal cations indicate that for pK a calculations, an accuracy of 2 pK a units can be achieved by GGA functionals (i.e.generalized gradient approximation) 28,[50][51][52] and no noticeable differences have been found for GGA and hybrid functionals. 45][55][56][57] The errors correlate with the proximity of the redox potentials and the valence band position; the more positive the redox potentials, i.e. closer to the water valence band, the larger the errors.The accuracy can be improved by using hybrid functionals mixing a fraction of the Hartree-Fock exchange, such as the HSE06 functional. 58,59For the couples with low reduction potentials (e.g., CO 2 À /CO 2 ), GGA functionals have a similar accuracy to hybrid functionals.Using this method, we calculated the acidity constants of molybdic acid, arsenites and thioarsenites from room temperature to 573 K, and the results showed that the computed pK a s agree with experiment within 2 pK a units. 28,60,61These studies give us confidence that the effect of temperature on acidity constants can be reproduced at a reasonable accuracy with the present setup.
In this study, we employ FPMD techniques to investigate the structures, acidity constants and reduction potentials of UO 2 2+ at elevated temperatures up to 573 K.The hydration structures of UO 2 2+ and UO 2 + have been characterized in detail.The pK a s and redox potentials are calculated using the FPMD based vertical energy gap technique.The favorable agreement with available experimental measurements is encouraging for extending the current method for studying the speciation of f-block elements.

Models
The unit cell for all the simulations in this work is a cubic box of 12.43 Å under full 3D periodic boundary conditions.Under ambient conditions, the cell contains 63 H 2 O molecules, which approximately corresponds to the density of liquid water.The number of water molecules at each elevated temperature (Table 1) reproduces the density of the liquid phase at the corresponding saturated vapor pressure. 62The initial configurations were prepared by placing UO 2 (H 2 O) 5 2+/+ complexes at the center of the cell and then inserting appropriate numbers of water molecules around the complexes, and then optimizing and equilibrating with MD simulations.

FPMD details
The FPMD simulations were performed using the freely available CP2K/QUICKSTEP package. 26In QUICKSTEP, the electronic structures are calculated with density functional theory implemented based on a hybrid Gaussian plane wave (GPW) approach. 63The PBE functional was used in this work. 64The core electrons are represented by analytic Goedecker-Teter-Hutter (GTH) pseudopotentials, 65,66 and double-z basis sets augmented with polarization functions were employed for H, O and U.The pseudopotential and the basis set for U have been successfully applied on a uranium dioxide bulk system. 67,68The plane wave density cutoff was set to be 400 Ry.Born-Oppenheimer molecular dynamics (BOMD) simulations were carried out with a time step of 0.5 fs.The temperature was controlled using the Nose ´-Hoover chain thermostat.For each simulation, the production run was performed for at least 5.0 ps following a prior equilibration for at least 2.0 ps.

pK a calculations
Using the vertical energy gap method, the proton of the acid (denoted AH) is gradually transformed into a dummy atom (i.e. a classical particle with no charge) and the free energy of this transformation (denoted D dp A AH ) is evaluated by using the thermodynamic integration technique (Section 2.5).The same procedure is applied to transform a proton of a hydronium located in a simulation cell of the same size into a dummy.The corresponding deprotonation integral is denoted The formula used to calculate pK a reads: is the unit molar concentration and L H + is the thermal wavelength of a proton.The third term k B T ln[c 0 L H + 3 ] accounts for the translational entropy generated by the acid dissociation, which is approximated by the free energy of a free proton at the standard concentration.This term is equal to À0.19 eV under ambient conditions. 69

Redox potential
The redox potential of a redox couple (A À /A ) with respect to SHE is calculated with: The term of D ox A A À stands for the free energy of the reversible removal of an electron from the reduced state, A À (aq) -A (aq) + e À (vac), which is calculated with the thermodynamic integration (Section 2.5).D f G Ã H þ ðgÞ stands for the formation free energy of the gas-phase proton and the values at different temperatures can be found from thermodynamics tables 70 and listed in Table 3. D zp E H(OH 2 ) + denotes the zero point energy correction.

Free energy perturbation method
The free energy changes for the electron transfer (ET) and proton transfer (PT) reactions are calculated with a combination of FPMD and thermodynamic integration.According to this implementation, the reactant is gradually transformed into the product using an auxiliary Hamiltonian: Here H R and H P stand for the reactant and product states, respectively.Where Z is the coupling parameter which is increased from 0 (reactant) to 1 (product), that is, for the ET reaction, it is from the reduction state (M + ) to the oxidation state (M 2+ ) and for the PT reaction, it is from the protonated state (HA) to the deprotonated state (A À ).Intermediate values for 0 o Z o 1 correspond to the hybrid systems of reactant and product states, which have no physical counterpart.V r means the restrained harmonic potential, which reads, For ET reactions, V r is usually not needed.For PT reactions, this potential is used to restrain the dummy atom in a location resembling that of the acid proton of the reactant state through the bonding and angle bending whose equilibrium values are d 0 and a 0 , respectively.The equilibrium values used are obtained from the prior free simulations (i.e.without restraints) and the details of V r are shown in Table 2.
The free energy change of the transformation is calculated from the integral of the vertical energy gap with respect to the coupling parameter: The vertical energy gap is defined as the potential energy difference between initial and final configurations, which is calculated from MD trajectories.
For ET reactions, the energy gap is the vertical ionization energy and the oxidation free energy is calculated from, For PT reactions, the energy gap is the vertical deprotonation energy and the deprotonation free energy of the acid HA is calculated with: Similarly, for the deprotonation of hydronium, In practice, the Simpson rules can be used to calculate the integral, e.g. the 3-point formula with Z = 0, 0.5, and 1, respectively In a summary, the pK a and E 0 are calculated with the following two equations, respectively In those equations, the integrals are computed by using the thermodynamic integration.The used value of D zp E H(OH 2 ) + is 0.35 eV, which is obtained from a previous study. 47The values of the other terms including are summarized in Table 3.The vertical energy levels can also be aligned with respect to the SHE by using eqn (11), so that they can be compared with the band edges of liquid water. 48In that case, the EA (electron affinity) and IP (ionization potential) levels of the solutes are written as, As shown in our previous publications, 45,47,52,54 finite size effects are not the dominant source of error due to the effective screening of liquid water.
Screening in polar liquids is taken into account by the Born cavity model.With the extension under periodic boundary conditions by Hummer et al., [71][72][73][74][75] the correction to the free energy has the following form: Here L is the side length of the cubic cell, q is the charge of the ion, x EW is the Madelung constant, e is the dielectric constant of water and R is the ionic radius.
The first term effectively vanishes for aqueous ions due to the high dielectric constant (e 4 80).The second term determines in practice the finite size error, i.e. the error for the oxidation or hydrolysis of M q+ is proportional to (q + 1) 2 À q 2 (or equivalently q 2 À (q À 1) 2 ).The correction to the redox free energy of the reaction with the net charge change from 1 to 0 such as OH /OH À is estimated to be less than 0.1 eV with a box of 9.86 Å, which is within the statistical uncertainty of the FPMD simulations.Notice that in the formulae for pK a and E 0 calculations (eqn ( 10) and ( 11)), there is an effective error cancellation with the deprotonation reaction of H 3 O + (its charge changes from +1 to 0).Such an approach also applies for the reactions for M 2+ /M 1+ or M 2À / M À .This has been evidenced by our calculations, e.g. for pK a for HS À (i.e.HS À /S 2À ), 45 the error is found to be less than 0.1 eV (17.1 vs. 17.0).Here for UO 2

2+
, with a side length of 12.43 Å, L 3 is twice larger than that with a side length of 9.86 Å.The finite size error is further reduced, and therefore we did not include any correction to the results.

Hydration structures UO 2
2+ has 5 water molecules in the first hydration shell from room temperature to 573 K (Fig. 1A and 2A).For UO  2B and the RDF-CN in Fig. 4B) and at higher temperatures, the average CN of H-bonds decreases to 1.

Acidity constants
Table 4 lists the vertical energy gaps, thermodynamic integrals and pK a s for UO 2 2+ .It can be seen that all of the vertical energy gaps converge within 0.1 eV and the statistical errors in pK a s are within 1 pK a unit.Under ambient conditions, the calculated pK a is 5.8, in good agreement with the latest experiment, 5.58. 8he pK a at 373 K is calculated to be 4.0, which is very close to the values measured at slightly lower temperatures, 4.19 at 367.4 K 85 and 4.24 at 358.15 K. 8 The pK a values show a decreasing trend with temperature.A similar trend has also been found for other complexes, such as Fe(H 2 O) 3+ 86 and H 3 AsO 3 .61 At 473 K and 573 K, the pK a s are 2.8 and 1.4, respectively.

Redox potentials
The vertical energy gaps for reduction potential calculations are listed in Table 5.The energy levels of the solutes and liquid water are summarized in Table 6 and illustrated in Fig. 5.As shown in our previous study, PBE places the VBM (valence band maximum) of liquid water 3.2 V too high while it estimates the CBM (conduction band minimum) much better (by only 0.5 V too high). 48Due to the hybridization with the too high VBM of liquid water predicted by PBE, the IP level of the solute would also be at a too high position, which then leads to underestimation of the reduction potential.This error is more severe for the more oxidative solutes because the pining with the VBM of liquid water is worse. 53For UO 2 2+ /UO 2 + , the IP level predicted by PBE is 0.95 V (Table 6), over 1.3 V above the PBE water VBM, indicating that the negative effect of the band misalignment is relatively small on the solute levels.This is further confirmed by the energy level diagram (Fig. 5), where it clearly shows a symmetric distribution of IP and EA levels around the redox level (i.e.normal linear response from the Marcus picture 88 ).As previously shown, 53 pronounced asymmetry in reorganization energies (energy differences between vertical and redox levels) is a sign of the energy levels suffering from large delocalization errors.
We believe that this explains why the redox potential computed by PBE functional is reasonably in line with experiment (À0.05 V vs. 0.062 V) at room temperature, and the favorable agreement lends us confidence in the PBE results at elevated temperatures.It is also clear that the symmetric alignment of the energy levels remains at elevated temperatures (Table 5).As temperature increases, the redox potential gradually increases to 0.25 V at 573 K.The increasing trend is the same as that found for the Fe 2+ /Fe 3+ couple. 60However, no experimental data are available for uranyl ions at elevated temperatures.

Discussion
It is believed that the crustal fluid is reductive because of the lack of oxygen.So, the increasing trend of redox potentials with the temperature indicates that the reduced state will get more populated at elevated temperatures compared to under ambient conditions.By considering the disproportionation nature of UO 2 + , aqueous U(IV) could be abundant in geo-fluids at higher temperatures.The acidity constants of UO 2 2+ get very low beyond 373 K (i.e.pK a s o 3.0 at 473 K and 573 K).This indicates that under those conditions, the dominant species are hydrolytic complexes instead of UO 2 2+ in the near-neutral pH (the pH of neutral water is around 5.6 at both temperatures 89 ).These findings should be taken into account in future research, such as the surface complexation of uranium.For example, a temperature of 473 K is possible at deep geological storage sites for high level nuclear wastes.However, previous studies only focus on the adsorption of UO 2 2+ whereas little has been done on its hydrolytic forms or the reduced states. 90t has been suggested that at elevated U concentration, polymeric complexes such as (UO 2 ) 2 (OH) 2 2+ and (UO 2 ) 3 (OH) 5 2+ may be formed. 5,91The entropic contribution gets more significant with increasing T, which is unfavorable to polymeric structures.However, there is no data available for temperatures over 373 K, and therefore it is still unclear if polymerization plays an important role in hydrothermal solutions.The free energies of uranyl ion complexes formed with other ligands (e.g.F À and CO 3 2À .)are also lacking at high temperatures.As shown in the previous studies, [38][39][40][41] FPMD based methods can provide reasonable estimates for thermodynamic data, and thus we expect the method presented here can be extended to elevated T-P conditions.Such complexing may change pK a s and U 0 and this can also be addressed by using the FPMD vertical energy gap technique.These unresolved issues will merit future investigations.
Many key thermodynamic parameters are still lacking for important actinides even under ambient conditions, e.g., hydrolysis constants of Ac(IV) are still poorly documented.At elevated temperatures, only limited pK a data are available for actinides and no measurement of reduction potentials has    12) and ( 13), respectively.E 0 is taken from Table 5. b The VBM is obtained by substituting the integral in eqn ( 13) with the vertical ionization potential of the water.The CBM is estimated based on orbital energies, i.e. by replacing the integral in eqn (12) with minus the average energy of the LUMO (the lowest unoccupied molecular orbital).c Ref. 48. , IP: vertical ionization potential of UO 2 been performed. 92Another issue is the pH-dependent stoichiometry of Ac(VI)-carbonate complexes.It has been found that as pH increases, carbonate complexes get more significant for UO 2 2+ , but the detailed mechanism has not been quantified. 7 the other hand, kinetic factors can play an important role and need to be taken into account to fully understand the speciation.The FPMD-based technique provides a feasible way for quantifying these data and mechanisms.
The results above have shown that acidity constants and reduction potentials of uranyl ions are reproduced within 1.0 pK a unit and 0.2 V, respectively.This error margin is similar to that found for the species of main group elements and transition metal cations in our previous calculations. 45For pK a s, due to the closed shell nature of proton transfer reactions, GGA functionals can give good prediction and hybrid functionals do not show any obvious difference. 45So, the error in pK a calculation is almost all the statistical error.For E 0 , the GGA result is reasonable for uranyl ions because their reduction potential is relatively low.For the redox couples of higher reduction potentials, GGA causes serious underestimation whereas hybrid functionals can give better results.This is because hybrid functionals improve the alignment of the VBM of liquid water, 87 which effectively reduces the errors in computed redox potentials.For example, the Ag + /Ag 2+ couple having a reduction potential as high as 1.98 V is an extreme test case for the present techniques.The GGA result is 1.15 V whereas HSE06 is rather satisfactory, being 1.72 V. 54 In light of this promising progress, we are optimistic on the application of the current technique for predicting pH-Eh diagrams of the elements of geological interest under the conditions of hydrothermal fluids.The present work on the U complexes is the first application to f-block species, and the encouraging results indicate that the method is ready to calculate the acid-base and redox chemistry of f-block species at a reasonable accuracy.Finally, we are aware that the UO 2 2+ /UO 2 + couple has electronic configurations of f 0 /f 1 , a relatively safe situation for DFT.Complex electron correlation will be a challenge for the electronic structure theory when there are multiple electrons in f orbitals, which will have to be addressed in future work.

Summary
In this study, we investigate the structures, acidity constants and redox potentials of uranyl ions from ambient conditions to 573 K by using FPMD techniques.It is found that UO 2 2+ holds five H 2 O ligands up to 573 K, whereas UO 2 + keeps 5-coordinate up to 473 K and loses one H 2 O ligand at 573 K.The FPMD based vertical energy gap method is applied to calculate pK a s and E 0 s.The results show that the pK a s and E 0 can be reproduced with an accuracy of 1 pK a unit and 0.2 V, respectively.This demonstrates that the current methodology is able to predict the aqueous speciation of f-block elements at a good accuracy under the conditions of geological interest.The pK a s of the aqueous uranyl ions decrease with temperature and are lower than 3 at temperatures above 473 K, indicating that hydrolytic forms should be dominative.The reduction potentials show an increasing trend with temperature, implying that the reduced states may play more significant roles as T increases.The computed data can be directly used in future studies for understanding the transport and fixation of uranium.

Fig. 1
Fig. 1 RDFs (radial distribution functions) and CNs (coordination numbers) for U-O derived from the simulations of UO 2 2+ and UO 2 + .

Fig. 2
Fig. 2 Snapshots of UO 2 2+ and UO 2 + derived from the simulations at 573 K. Some water molecules have been removed for clarity.O = red, H = white and U = blue.

Table 1
The T-P conditions and numbers of water molecules for simulations.The numbers include the water molecules in the first shell of uranyl ions

Table 2
The parameters used in the harmonic potentials (eqn (4)) restraining the dummy particles for UO 22+.Equilibrium bond lengths (d 0 ) are in Bohr and equilibrium angles (a 0 ) are in radians.All the coupling constants are in a.u.H d in the table means the dummy particle

Table 3 D
dp A H 3 O +, k B T ln[c 0 L H + 3 ] and ÀD f G Ã H þ ðgÞ at different temperatures.The energy unit is eV

Table 4
8he computed vertical energy gaps (in eV), thermodynamic integrals (in eV) and pK a s for UO 2 2+ at various temperatures.The numbers in Statistical error in the vertical energy gap is calculated as the semi-difference between the values using the first half or the second half of the trajectory only.aThe values were measured at slightly lower temperatures than 373 K, 4.19 at 367.4 K 85 and 4.24 at 358.15 K.8

Table 5
The computed vertical energy gaps (in eV), thermodynamic integrals (in eV) and redox potentials (in V) for the UO 2 Statistical error in the vertical energy gap is calculated as the semi-difference between the values using the first half or the second half of the trajectory only.

Table 6
a EA and IP of UO 2 2+ /UO 2 + are obtained by using eqn (