CO2 solubility in NaCl solutions from 25 to 90 °C : Reaktoro versus PHREEQC
Reaktoro is "an open source computational framework developed in C++ and Python to simulate chemically reactive processes" (reaktoro.org). Its website provides an application example entitled : "Carbon-dioxide solubility in brines with different salinity and temperature", for which the Jupyter notebook is given below left.
In this example, 10 moles of CO2(g) at 100 bar (the total/overall pressure of the system being 100 bar) are equilibrated with 1 kg of water, with [NaCl] either of 1, 2 or 4 mol/kgw, and at temperature varying from 25 to 90 °C. The execution of the notebook generates a plot showing the CO2 concentration (actually the C concentration) in the three aqueous solutions of NaCl, as function of temperature (see here).
With Reaktoro, a Python nested for loop over NaCl molality and temperature is used (3rd cell, see here). The same example, calculated with PHREEQC, needs a script (below right) with 3 subsequent simulations including the following keywords : SOLUTION, REACTION_PRESSURE, EQUILIBRIUM_PHASES, REACTION_TEMPERATURE, and a USER_GRAPH to generate a plot incorporating the Reaktoro's results (exportable with print(T, molNaCl, float(aqprops.elementMolality("C"))) at the end of the 3rd cell and downloadable here).
In order to the Reaktoro and PHREEQC calculations to be comparable, the following conditions are fullfilled :
- Same database (phreeqc.dat) (the Reaktoro embedded phreeqc.dat database is available in the Reaktoro's GitHub repository (link here))
- The
aqueousphase.set(ActivityModelPhreeqc(db))line in the Reaktoro script "ensures consistency with PHREEQC by using the same activity model PHREEQC would use with the phreeqc.dat database" (information given here and here) - Same activity (fugacity) model for the gas phase : the Peng-Robinson Equation of State (EoS) that is implemented in PHREEQC (Appelo et al., 2014) ; ensured with the
gaseousphase.set(ActivityModelPengRobinsonPhreeqcOriginal())line in the Reaktoro script (information given here and here)
The results given by Reaktoro and PHREEQC are quite similar (see graph below). The slight differences may be due to the different methods/approaches to calculate equilibrium : "Reaktoro is based on the Gibbs Energy Minimization (GEM) methods as well as a revised Law of Mass Action (rLMA) approach" (geg.ethz.ch) while PHREEQC is based on the LMA method only (aquion.de).
Reaktoro script
from reaktoro import *
# Initialize a thermodynamic database
db = PhreeqcDatabase("phreeqc.dat")
# Create an aqueous phase automatically selecting all species with provided elements
aqueousphase = AqueousPhase(speciate("H O C Na Cl"))
aqueousphase.set(ActivityModelPhreeqc(db))
# Create a gaseous phase with CO2(g)
gaseousphase = GaseousPhase("CO2(g)")
gaseousphase.set(ActivityModelPengRobinsonPhreeqcOriginal())
# Create the chemical system
system = ChemicalSystem(db, aqueousphase, gaseousphase)
# Create the equilibrium solver
solver = EquilibriumSolver(system)
import numpy as np
import pandas as pd
temperatures = np.arange(25.0, 90.0, 5.0)
molsNaCl = np.array([1.0, 2.0, 4.0])
P = 100.0
df = pd.DataFrame(columns=["T", "amountNaCl", "amountCaq"])
for molNaCl in molsNaCl:
for T in temperatures:
# Initial amount of the CO2 gas
n0CO2g = 10.0
# Define initial chemical state corresponding to the NaCl-brine of the given concentration
state = ChemicalState(system)
state.setTemperature(T, "celsius")
state.setPressure(P, "bar")
state.set("H2O" , 1.0 , "kg")
state.set("CO2(g)", n0CO2g, "mol")
state.set("Na+" , molNaCl , "mol")
state.set("Cl-" , molNaCl , "mol")
# Calculate equilibrium state
res = solver.solve(state)
# Stop if the equilibration did not converge or failed
if not res.succeeded(): continue
# Fetch resulting aqueous properties of the chemical state
aqprops = AqueousProps(state)
# Update value ["T", "amountNaCl", "amountCaq"] in the dataframe
df.loc[len(df)] = [T, molNaCl, float(aqprops.elementMolality("C"))]
from reaktplot import *
fig = Figure()
fig.title("SOLUBILITY OF CO2 IN NACL BRINES")
fig.xaxisTitle('TEMPERATURE [°C]')
fig.yaxisTitle('AMOUNT OF DISSOLVED CO2 [mol/kgw]')
for molNaCl in molsNaCl:
df_NaCl = df[df['amountNaCl'] == molNaCl]
fig.drawLineWithMarkers(df_NaCl["T"], df_NaCl["amountCaq"], name=f'{molNaCl} mol of NaCl')
fig.show()
PHREEQC script
DATABASE phreeqc.dat_Embedded-in-Reakt.dat # https://github.com/reaktoro/reaktoro/blob/main/embedded/databases/phreeqc/phreeqc.dat
SOLUTION
Na 1000
Cl 1000
REACTION_PRESSURE
98.6923 # 100 bar
EQUILIBRIUM_PHASES
CO2(g) 1.99428 10 # 2nd column = log 98.6923
REACTION_TEMPERATURE
25 30 35 40 45 50 55 60 65 70 75 80 85 90
USER_GRAPH
-headings 1.0_mol_of_NaCl_PHREEQC
-chart_title "SOLUBILITY OF CO2 IN NACL BRINES"
-axis_titles "TEMPERATURE [°C]" "AMOUNT OF DISSOLVED CO2 [mol/kgw]"
-axis_scale x_axis 25 90 auto auto
-axis_scale y_axis 0.2 1.2 auto auto
010 plot_xy tc, tot("C") symbol_size = 0
-plot_tsv_file Reaktoro_data_1mNaCl.tsv
END
# ================================
SOLUTION
Na 2000
Cl 2000
REACTION_PRESSURE
98.6923 # 100 bar
EQUILIBRIUM_PHASES
CO2(g) 1.99428 10 # 2nd column = log 98.6923
REACTION_TEMPERATURE
25 30 35 40 45 50 55 60 65 70 75 80 85 90
USER_GRAPH
-headings 2.0_mol_of_NaCl_PHREEQC
010 plot_xy tc, tot("C") symbol_size = 0
-plot_tsv_file Reaktoro_data_2mNaCl.tsv
END
# ================================
SOLUTION
Na 4000
Cl 4000
REACTION_PRESSURE
98.6923 # 100 bar
EQUILIBRIUM_PHASES
CO2(g) 1.99428 10 # 2nd column = log 98.6923
REACTION_TEMPERATURE
25 30 35 40 45 50 55 60 65 70 75 80 85 90
USER_GRAPH
-headings 4.0_mol_of_NaCl_PHREEQC
010 plot_xy tc, tot("C") symbol_size = 0
-plot_tsv_file Reaktoro_data_4mNaCl.tsv