# Effect of calcite dissolution on exchangeable cations¶

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

In this tutorial, we study the effect of calcite dissolution on determination of CaX2 using 1-molal NH4-acetate. We will vary the acidity of the solution from 9.3 to 7.0 by adding NH4+ to the porewater.

As all the tutorials above, the numerical experiment start from importing necessary python-packages, initializing database, and defining aqueous and exchange phases to form chemical system.

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

db = PhreeqcDatabase("phreeqc.dat")

# Define an aqueous phase
solution = AqueousPhase(speciate("H O C Na N Cl Ca Mg"))
solution.setActivityModel(chain(
ActivityModelHKF(),
ActivityModelDrummond("CO2")
))


Since we are mainly concerned with the extraction of Ca+2, we only include two ion exchange species into ion exchange phase. Below, we also define a single mineral phase represented by calcite.

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

# Define an ion exchange phase
mineral = MineralPhase("Calcite")


After all the phases are prepared, the chemical system as well all other building blocks of equilibrium simulation, i.e., equilibrium solver and aqueous properties, can be initialized.

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

# Define equilibrium solver and equilibrate given initial state with input conditions
solver = EquilibriumSolver(system)

# Define ion exchange properties and aqueous properties
aqprops = AqueousProps(system)


Bellow, we define the range of NH4+ amount as well as the array with NaX values in the initial chemical state. We also collect the parameters and arrays of the data we are trying to study.

# Sampling arrays of NH4 ions' amounts and NaX exchange species
num_ph = 31
num_exchangers = 4
mols_NH4 = np.linspace(0, 60.0, num=num_ph)
mols_NaX = 1e-3 * np.array([0.0, 2.5, 12.5, 25])

# Output dataframe
data = pd.DataFrame(columns=["amount_NaX", "amount_NH4", "pH", "I",
"amount_Ca", "amount_CaX2", "delta_Calcite"])


Next, we perform a sequence of equilibrium calculations by varying amounts of NH4-acetate and NaX. At the end of each calculation, we’ll extract the following properties from the computed chemical state:

• pH,

• ionic strength I,

• amount of Ca+2,

• amount of CaX2,

• amount of Calcite.

During each equilibrium calculation, the following two major equilbrium calculations will be defining the solubility of calcite and to formation of CaX2: \begin{alignat}{4} {\rm CaCO{_3}} + {\rm H_2O} & \rightleftharpoons {\rm Ca}^{+2} + {\rm HCO}_3^- + {\rm OH}^- \\ \tfrac{1}{2} {\rm Ca^{+2}} + {\rm NaX} & \rightleftharpoons {\rm Na}^+ + \tfrac{1}{2} {\rm CaX}_2\\ \end{alignat}

for mol_NaX in mols_NaX:

for mol_NH4 in mols_NH4:

# Initial amount of Calcite
m0Calcite = 10.0

# Define initial chemical state
state = ChemicalState(system)
state.setTemperature(25.0, "celsius")
state.setPressure(1.0, "atm")
# Seawater
state.set("H2O" , 1.0    , "kg")
state.set("Na+" , 1.10   , "mmol")
state.set("Mg+2", 0.48   , "mmol")
state.set("Ca+2", 1.90   , "mmol")
state.set("NH4+", mol_NH4, "mmol")
# Ammonia
state.set("Calcite", m0Calcite, "mol")

# Equilibrate chemical state corresponding to the seawater
res = solver.solve(state)
# Stop if the equilibration did not converge or failed
if not res.optima.succeeded: continue

# Update aqueous properties and evaluate pH
aqprops.update(state)
pH = float(aqprops.pH())

# Exchanger
state.set("NaX", mol_NaX, "mol")

# Equilibrate the seawater with carbonates
res = solver.solve(state)
# Stop if the equilibration did not converge or failed
if not res.optima.succeeded: continue

# Update aqueous properties to evaluate ionic strength
aqprops.update(state)
chemprops = state.props()

# Collect the value to be added to the dataframe in the following order
# "amount_NaX", "amount_NH4", "amount_pH", "I", "amount_Ca", "amount_CaX2", "delta_Calcite"
data.loc[len(data)] = [mol_NaX, mol_NH4, pH, float(aqprops.ionicStrength()),
float(chemprops.elementAmountInPhase("Ca", "AqueousPhase")), float(state.speciesAmount("CaX2")),
m0Calcite - float(state.speciesAmount("Calcite"))]


To visualize the obtained results, we use bokeh plotting library.

from bokeh.plotting import figure, show
from bokeh.models import HoverTool
from bokeh.io import output_notebook
from bokeh.models import ColumnDataSource
output_notebook()


First, we plot the dependence of the dissolved calcite on the levels of pH and for different initial values of NaX providing the exchanger X-. We see that the solubility of calcite grows with an increase in acidity.

hovertool = HoverTool()
hovertool.tooltips = [("pH", "@pH"),
("delta(Calcite)", "@delta_Calcite mol"),
("amount(NaX)", "@amount_NaX mol")]

p = figure(
title="DEPENDENCE OF CALCITE SOLUBILITY ON PH",
x_axis_label=r'PH [-]',
y_axis_label='AMOUNT OF DISSOLVED CALCITE [MOL]',
sizing_mode="scale_width",
plot_height=300)

colors = ['teal', 'darkred', 'indigo', 'coral']
for mol_NaX, color in zip(mols_NaX, colors):
df_NaX = ColumnDataSource(data[data['amount_NaX'] == mol_NaX])
p.line("pH", "delta_Calcite", legend_label=f'{mol_NaX} mol of NaX', line_width=3, line_cap="round", line_color=color, source=df_NaX)

p.legend.location = 'top_right'
show(p)


Finally, we represent the amount of CaX2 resulting from the reaction NaX + Ca+2 = CaX2 + Na+, where the level of Ca+2 ions grows from the dissolved calcite. We see that at lower pH values, we obtain higher amounts of precipitated CaX2.

hovertool = HoverTool()
hovertool.tooltips = [("pH", "@pH"),
("delta(Calcite)", "@delta_Calcite mol"),
("amount(NaX)", "@amount_NaX mol")]

p = figure(
title="DEPENDENCE OF CAX2 DISTRIBUTION ON PH",
x_axis_label=r'PH [-]',
y_axis_label='CAX2 AMOUNT [MOL]',
sizing_mode="scale_width",
plot_height=300)