Notebook

Effect of calcite dissolution on exchangeable cations¶

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 effect of calcite dissolution on determination of CaX₂ using 1-molal NH₄-acetate. We will vary the acidity of the solution from 9.3 to 7.0 by adding NH₄⁺ 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.

In [1]:

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"))
solution.set(ActivityModelPitzer())

Out[1]:

<reaktoro.reaktoro4py.AqueousPhase at 0x7f45105ecfb0>

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

In [2]:

# Define an ion exchange phase
exchange = IonExchangePhase("NaX CaX2")
exchange.set(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.

In [3]:

# 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 NH₄⁺ 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.

In [4]:

# 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 NH₄-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⁺²,
amount of CaX₂,
amount of Calcite.

During each equilibrium calculation, the following two major equilbrium calculations will be defining the solubility of calcite and to formation of CaX₂: $$ \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} $$

In [5]:

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

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.

In [6]:

from reaktplot import * 

fig = Figure()
fig.title("DEPENDENCE OF CALCITE SOLUBILITY ON PH")
fig.xaxisTitle('pH')
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')

fig.show()