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

Attention

Always make sure you are using the latest version of Reaktoro. Otherwise, some new features documented on this website will not work on your machine and you may receive unintuitive errors. Follow these update instructions to get the latest version of Reaktoro!

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.set(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 UO₂(OH)₂.

# 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=101)

# 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])

***WARNING***
The chemical equilibrium/kinetics calculation did not converge.
This can occur due to various factors, such as:
  1. Infeasible Problem Formulation: The problem you have defined lacks a feasible solution within the given constraints and assumptions.
  2. Out-of-Range Thermodynamic Conditions: The conditions you have specified for the system may fall outside the valid range supported by the thermodynamic models used.
  3. Numerical Instabilities: Convergence issues may arise from numerical problems during the execution of the chemical/kinetic equilibrium algorithm. Consider reporting the issue with a minimal reproducible example if you believe the algorithm is responsible for this issue.
Disable this warning message with Warnings.disable(906) in Python and Warnings::disable(906) in C++.

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

from reaktplot import *

fig1 = Figure()
fig1.xaxisTitle('pH')
fig1.yaxisTitle('U(VI) Speciation [%]')

fig2 = Figure()
fig2.xaxisTitle('pH')
fig2.yaxisTitle('log(U(VI)) Speciation [mol]')


for species in species_list_str:
    fig1.drawLine(df["pH"], df["perc_"+species], name=species)
    fig2.drawLine(df["pH"], df["amount_"+species], name=species)

fig1.height(800)
fig2.height(800)

fig1.show()
fig2.show()

The plot above shows uranium speciation in % with changing pH values at fixed CO2 partial pressure. We see that for different acidity/alkalinity ranges different species are dominant.