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(ActivityModelPhreeqc(db))

# 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.00579097	0.094209	1e-16	0.094209	0.305791	1e-16	6.98671
1	0.095	0.025	0.00873692	0.134683	0.000790228	0.0862631	0.265317	0.0242098	6.98361
2	0.09	0.05	0.011732	0.171596	0.00333599	0.078268	0.228404	0.046664	6.98079
3	0.085	0.075	0.0144722	0.202807	0.00886037	0.0705278	0.197193	0.0661396	6.97806
4	0.08	0.1	0.0166446	0.227346	0.0180048	0.0633554	0.172654	0.0819952	6.97538
5	0.075	0.125	0.0181137	0.245914	0.0304859	0.0568863	0.154086	0.0945141	6.97277
6	0.07	0.15	0.0189176	0.259953	0.0455647	0.0510824	0.140047	0.104435	6.97029
7	0.065	0.175	0.0191628	0.270776	0.0625305	0.0458372	0.129224	0.112469	6.96797
8	0.06	0.2	0.0189568	0.279337	0.0808531	0.0410432	0.120663	0.119147	6.96583
9	0.055	0.225	0.0183882	0.286279	0.100166	0.0366118	0.113721	0.124834	6.96384
10	0.05	0.25	0.0175251	0.292034	0.120221	0.0324749	0.107966	0.129779	6.96199
11	0.045	0.275	0.0164195	0.296895	0.140843	0.0285805	0.103105	0.134157	6.96028
12	0.04	0.3	0.0151111	0.301068	0.161911	0.0248889	0.0989323	0.138089	6.95869
13	0.035	0.325	0.0136307	0.304698	0.183336	0.0213693	0.0953023	0.141664	6.95721
14	0.03	0.35	0.0120026	0.307892	0.205053	0.0179974	0.0921077	0.144947	6.95582
15	0.025	0.375	0.010246	0.310732	0.227011	0.014754	0.0892683	0.147989	6.95452
16	0.02	0.4	0.00837679	0.313277	0.249173	0.0116232	0.0867231	0.150827	6.9533
17	0.015	0.425	0.00640769	0.315575	0.271508	0.00859231	0.0844246	0.153492	6.95214
18	0.01	0.45	0.00434943	0.317665	0.293993	0.00565057	0.0823352	0.156007	6.95105
19	0.005	0.475	0.00221099	0.319575	0.316607	0.00278901	0.0804249	0.158393	6.95002
20	2e-16	0.5	1e-16	0.32133	0.339335	1e-16	0.0786695	0.160665	6.94904

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()