Cation exchange competition among Na⁺, K⁺ and Ca²⁺

Written by Svetlana Kyas (ETH Zurich) on Mar 10th, 2022.
Last revised on Mar 16, 2023 by Allan Leal.

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!

In this tutorial, we show how Reaktoro can be used to model ion exchange processes. We use the PHREEQC database phreeqc.dat for this purpose. We assume that there is an ion exchanger X^- on which the aqueous cations Na⁺, K⁺ and Ca²⁺ can bind and form the exchange phases NaX, KX and CaX₂ respectively. We are interested in calculating the amounts of these phases as the concentrations of K⁺ and Ca²⁺ vary in the aqueous solution.

First, we define a chemical system representative of the process we want to model.

import numpy as np
from reaktoro import *
import pandas as pd

# Initialize the PHREEQC database
db = PhreeqcDatabase("phreeqc.dat")

# Define an aqueous phase
solution = AqueousPhase("H2O Na+ Cl- H+ OH- K+ Ca+2 Mg+2")
solution.set(ActivityModelHKF())

# Define an ion exchange phase
exchange = IonExchangePhase("NaX KX CaX2")
exchange.set(ActivityModelIonExchangeGainesThomas())

# Create the chemical system
system = ChemicalSystem(db, solution, exchange)

As we are going to perform a sequence of chemical equilibrium calculations, we define an object of EquilibriumSolver as well as IonExchangeProps and AqueousProps using the chemical system. These objects will be updated after every equilibrium calculation.

# Define the equilibrium solver
solver = EquilibriumSolver(system)

# Define exchange, aqueous, and chemical properties
exprops = IonExchangeProps(system)
aqprops = AqueousProps(system)

Note

Note that at the end of each equilibrium calculation, state.props() can be used to access a ChemicalProps object with the chemical properties of the system evaluated in the last iteration of the algorithm.

Next, we define the range of Ca²⁺ and K⁺ initial amounts as well as auxiliary pandas.DataFrame to collect the results of the ion exchange simulations.

# Output dataframe
df = pd.DataFrame(columns=[
    "amount_K", 
    "amount_Ca",
    "amount_K_ion", 
    "amount_Na_ion", 
    "amount_Ca_ion",
    "amount_KX", 
    "amount_NaX", 
    "amount_CaX2",
    "pH"
])

# Sampling arrays storing the amounts of ions
steps = 21

mols_K  = np.linspace(0.1, 0.0, num=steps)
mols_Ca = np.linspace(0.0, 0.5, num=steps)

We’ll now perform a sequence of equilibrium calculations by varying the initial amounts of ions Ca²⁺ and K⁺ ions (taken from the arrays mols_Ca and mols_K, respectively). At the end of each calculation, we’ll extract the following properties from the computed chemical state:

• amount of K,

• amount of Ca,

• amount of K⁺,

• amount of Na⁺,

• amount of Ca²⁺,

• amount of KX,

• amount of NaX,

• amount of CaX2,

• pH.

for mol_K, mol_Ca in zip(mols_K, mols_Ca):

    # Define initial chemical state
    state = ChemicalState(system)
    state.setTemperature(25, "celsius")
    state.setPressure(1, "atm")
    state.set("H2O" , 1.0   , "kg")
    state.set("NaX" , 0.4   , "mol")
    state.set("K+"  , mol_K , "mol")
    state.set("Ca+2", mol_Ca, "mol")

    # Equilibrate the chemical state
    res = solver.solve(state)
    
    # Stop if the equilibration did not converge or failed
    assert res.succeeded()

    # Update exchange and aqueous properties
    chemprops = state.props()

    exprops.update(chemprops)
    aqprops.update(chemprops)

    # Update output arrays:
    df.loc[len(df)] = [
        chemprops.elementAmount('K'),   # amount_K
        chemprops.elementAmount('Ca'),  # amount_Ca
        state.speciesAmount('K+'),      # amount_K+
        state.speciesAmount('Na+'),     # amount_Na+
        state.speciesAmount('Ca+2'),    # amount_Ca+2
        state.speciesAmount('KX'),      # amount_KX
        state.speciesAmount('NaX'),     # amount_NaX
        state.speciesAmount('CaX2'),    # amount_CaX2
        aqprops.pH()                    # pH
    ]

The data collected in the previous calculations can be inspected by printing the pandas.DataFrame object df:

df

	amount_K	amount_Ca	amount_K_ion	amount_Na_ion	amount_Ca_ion	amount_KX	amount_NaX	amount_CaX2	pH
0	0.1	2e-16	0.00571081	0.0942892	1e-16	0.0942892	0.305711	1e-16	6.99674
1	0.095	0.025	0.00857519	0.134804	0.000810585	0.0864248	0.265196	0.0241894	6.99694
2	0.09	0.05	0.0114588	0.17168	0.0034304	0.0785412	0.22832	0.0465696	6.99732
3	0.085	0.075	0.014057	0.20273	0.00910655	0.070943	0.19727	0.0658934	6.99787
4	0.08	0.1	0.0160688	0.227017	0.0184573	0.0639312	0.172983	0.0815427	6.9986
5	0.075	0.125	0.0173796	0.245317	0.0311517	0.0576204	0.154683	0.0938483	6.99952
6	0.07	0.15	0.0180442	0.259116	0.04642	0.0519558	0.140884	0.10358	7.00062
7	0.065	0.175	0.0181777	0.269738	0.063542	0.0468223	0.130262	0.111458	7.00187
8	0.06	0.2	0.0178905	0.278134	0.0819878	0.0421095	0.121866	0.118012	7.00326
9	0.055	0.225	0.0172712	0.28494	0.101394	0.0377288	0.11506	0.123606	7.00477
10	0.05	0.25	0.016387	0.290582	0.121516	0.033613	0.109418	0.128484	7.00637
11	0.045	0.275	0.0152885	0.295349	0.142181	0.0297115	0.104651	0.132819	7.00806
12	0.04	0.3	0.0140141	0.299441	0.163273	0.0259859	0.100559	0.136727	7.00982
13	0.035	0.325	0.0125931	0.303002	0.184702	0.0224069	0.0969978	0.140298	7.01164
14	0.03	0.35	0.0110485	0.306137	0.206407	0.0189515	0.0938626	0.143593	7.01353
15	0.025	0.375	0.00939862	0.308925	0.228338	0.0156014	0.0910751	0.146662	7.01546
16	0.02	0.4	0.00765805	0.311425	0.250459	0.0123419	0.0885754	0.149541	7.01744
17	0.015	0.425	0.00583879	0.313683	0.272739	0.00916121	0.0863171	0.152261	7.01945
18	0.01	0.45	0.00395073	0.315736	0.295156	0.00604927	0.0842636	0.154844	7.02151
19	0.005	0.475	0.00200214	0.317615	0.317692	0.00299786	0.0823854	0.157308	7.02359
20	2e-16	0.5	1e-16	0.319341	0.340329	1e-16	0.0806589	0.159671	7.02571

We can now construct an interactive plot from this data:

from reaktplot import *

fig = Figure()
fig.title("AMOUNTS OF EXCHANGE SPECIES KX, CaX<sub>2</sub>, NaX")
fig.xaxisTitle('Ca<sup>2+</sup> [molal]')
fig.yaxisTitle('Amount [mol]')

fig.drawLineWithMarkers(df["amount_Ca"], df["amount_KX"], name="KX")
fig.drawLineWithMarkers(df["amount_Ca"], df["amount_CaX2"], name="CaX<sub>2</sub>")
fig.drawLineWithMarkers(df["amount_Ca"], df["amount_NaX"], name="NaX")

fig.show()