Effect of calcite dissolution on exchangeable cations#

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


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

# Load Phreeqc database
db = PhreeqcDatabase("phreeqc.dat")

# Define an aqueous phase
solution = AqueousPhase(speciate("H O C Na N Cl Ca Mg"))

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

# 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: $$

(5)#\[\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.succeeded(): continue

        # Update aqueous properties and evaluate pH
        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.succeeded(): continue

        # Update aqueous properties to evaluate ionic strength
        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"))]

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.

from reaktplot import * 

fig = Figure()
fig.yaxisTitle('Dissolved Calcite [mol]')

for mol_NaX in mols_NaX:
    df_NaX = data[data['amount_NaX'] == mol_NaX]
    fig.drawLineWithMarkers(df_NaX["pH"], df_NaX["delta_Calcite"], f'{mol_NaX} mol of NaX')