Uranium (VI) speciation as a function of pH at fixed CO2 partial pressure#

Written jointly by Svetlana Kyas (ETH Zurich) and Dan Miron (PSI) on April 4th, 2022

This tutorial shows an example of the calculation of uranium speciation with changing pH values at fixed CO2 partial pressure. This example uses the thermodynamic database psinagra-12-07, which is based on the Nagra/PSI Chemical Thermodynamic Data Base 01/01. The database supports the ongoing safety assessments for the planned low and intermediate-level (L/ILW) and high-level (HLW) radioactive waste repositories in Switzerland.

Note

If your main interest is in calculating thermodynamic properties rather than modeling chemical equilibria and kinetics, you should check out ThermoFun, an excellent project dedicated to this task.

First, we set up the chemical system with the aqueous phase only. We are only interested in the U(VI) species distribution in the aqueous solution.

from reaktoro import *

# Define the Thermofun database
db = ThermoFunDatabase ("psinagra-12-07")

# Define the aqueous phase
solution = AqueousPhase(speciate("C Cl H O P S U"))
solution.setActivityModel(chain(
ActivityModelHKF(),
ActivityModelDrummond("CO2")
))

# Define chemical system by providing a database and an aqueous phase
system = ChemicalSystem(db, solution)


The equilibrium specifications given below indicate that the equilibrium calculations are performed for fixed T, P, pH, and partial pressure of CO2. The corresponding equilibrium conditions and the equilibrium solver are initialized with the instance specs. The conditions will be provided to the chemical solver later when we call it to calculate equilibrium.

Note

A detailed explanation of EquilibriumSpecs and EquilibriumConditions as well as a variety of problems with different boundary conditions are presented in the tutorial Chemical equilibrium with constraints.

# Specify conditions needed to be satisfied at chemical equilibrium
specs = EquilibriumSpecs(system)
specs.temperature()
specs.pressure()
specs.pH()
specs.fugacity("CO2(g)")

# Define the value of conditions needed to be satisfied at the chemical equilibrium state
conditions = EquilibriumConditions(specs)
conditions.temperature(25.0, "celsius")
conditions.pressure(1.0, "bar")
conditions.fugacity("CO2(g)", 0.01, "bar")

# Define the equilibrium solver
solver = EquilibriumSolver(specs)


The initial chemical state is defined according to the following recipe:

• 1000 g of water and

• 1e-5 mol of UO2(OH)2.

Note

When the ChemicalState class is used to define the chemical composition, the names of the species (present in the thermodynamic database) are used. Alternatively, the class Material uses recipes of substances to define the chemical composition (see, for example, thermofun tutorial Addition of limestone to the cement recipe).

# Define initial equilibrium state
state = ChemicalState(system)
state.set("H2O@"     , 1e3,  "g")
state.set("UO2(OH)2@", 1e-5, "mol")


The amount of uranium can be calculated based on the chemical properties of this chemical state. We define the range of pH values and the list of uranium-containing species of interest in our evaluation. We also create the list coeffs_u, which contains the number of uranium atoms in each species needed to correctly scale the amount of U in the species.

# Calculate the amount of uranium element
props = ChemicalProps(state)
bU = props.elementAmount("U")[0]

# Defined an auxiliary array of pH values
import numpy as np
pHs = np.linspace(3, 10, num=71)

# Create list of species names, list of Species objects, and auxiliary amounts array
species_list_str = "UO2+2 UO2OH+ UO2(OH)2@ UO2CO3@ (UO2)3CO3(OH)3+ (UO2)2CO3(OH)3- (UO2)2(OH)+3 " \
"UO2(CO3)3-4 UO2(CO3)2-2 (UO2)2(OH)2+2 (UO2)3(OH)5+ (UO2)4(OH)7+".split()

percentages = np.zeros(len(species_list_str))
amounts     = np.zeros(len(species_list_str))

# Collect the coefficients of U in the selected species
coeffs_u = []
for species_str in species_list_str:
# Fetch species as well as the lists of elements and corresponding coefficients
species = system.species().getWithName(species_str)
symbols = species.elements().symbols()
coeffs = species.elements().coefficients()

# Loop over all the elements and select coefficients for uranium
for [symbol, coeff] in zip(symbols, coeffs):
if symbol == "U":
coeffs_u.append(coeff)

# Define dataframe to collect amount of the selected species
import pandas as pd
columns = ["pH"] \
+ ["amount_" + name for name in species_list_str] \
+ ["perc_" + name for name in species_list_str]
df = pd.DataFrame(columns=columns)


Finally, we run simulations in the for loop on the pH values that we provide to the equilibrium conditions. There we previously set the temperature, pressure and the CO2 gas partial pressure.

for pH in pHs:

# Set the value of pH for the current equilibrium calculations
conditions.pH(pH)

# Equilibrate the initial state with given conditions and component amounts
res = solver.solve(state, conditions)

# If the equilibrium calculations didn't succeed, continue to the next condition
if not res.optima.succeeded: continue

# Otherwise, calculate U(VI) speciation in mols and %
for j in range(0, len(species_list_str)):
amounts[j] = float(state.speciesAmount(species_list_str[j]))
percentages[j] = float(coeffs_u[j] * amounts[j]) / bU * 100

# Update dataframe with obtained values
df.loc[len(df)] = np.concatenate([[pH], amounts, percentages])


In the following, we plot the obtained values of U(VI) speciation in % (using bokeh plotting library).

from bokeh.plotting import figure, show
from bokeh.io import output_notebook
output_notebook()

colors = ['teal', 'darkred', 'indigo', 'coral', 'rosybrown', 'steelblue', 'seagreen', 'red', 'darkred', 'darkkhaki', 'cadetblue', 'indianred']

p1 = figure(
title="",
x_axis_label=r'pH [-]',
y_axis_label='U(VI) SPECIATION [%]',
sizing_mode="scale_width",
plot_height=300)

p2 = figure(
title="",
x_axis_label=r'pH [-]',
y_axis_label='log(U(VI)) SPECIATION [MOL]',
y_axis_type = 'log',
sizing_mode="scale_width",
plot_height=300)

for species, color in zip(species_list_str, colors):
p1.line("pH", "perc_"+species, legend_label=species, line_width=3, line_cap="round", line_color=color, source=df)
p2.line("pH", "amount_"+species, legend_label=species, line_width=3, line_cap="round", line_color=color, source=df)

p1.legend.location = 'center_right'
p2.legend.location = 'center_right'

show(p1)
show(p2)