Notebook

Ion exchange competition among NaX, KX, and CaX{{_2}}¶

Written by Svetlana Kyas (ETH Zurich) on Mar 10th, 2022

{attention}
Always make sure you are using the [latest version of Reaktoro](https://anaconda.org/conda-forge/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](updating_reaktoro_via_conda) to get the latest version of Reaktoro!

In this tutorial, we study the dependence of ion exchange species distribution on the initial amounts of K⁺ and Ca²⁺.

First, we set up the chemical system, chemical state, and other means of modeling ion exchange using Reaktoro.

In [2]:

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.

In [3]:

# 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.

In [4]:

# 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.flip(np.linspace(0, 0.1, num=steps))
mols_Ca = np.linspace(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.

In [5]:

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")
    # Exchanger site
    state.set("NaX" , 0.4   , "mol")
    # Changing Ca+2 and K+
    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.optima.succeeded

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

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

We'll use the bokeh plotting library next. First, we need to import it and initialize it to work with Jupyter Notebooks:

In [6]:

from reaktplot import *

fig = Figure()
fig.title("DEPENDENCE OF ION EXCHANGE SPECIES<br>CONCENTRATIONS ON ADDED CA+2")
fig.xaxisTitle('Ca+2 [mol]')
fig.yaxisTitle('Amount [mol]')

fig.drawLineWithMarkers(df["amount_Ca"], df["amount_KX"], name="KX")
fig.drawLineWithMarkers(df["amount_Ca"], df["amount_CaX2"], name="CaX2")
fig.drawLineWithMarkers(df["amount_Ca"], df["amount_NaX"], name="NaX")

fig.height(800)
fig.show()