Ion exchange in dune sand and groundwater#

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

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 calculate the exchanger composition in contact with groundwater. To represent the exchanger, we consider 6 g of dune sand with 10 meq/kg of CEC (cation exchange capacity), which results in 0.06 mol of exchanger X (60 meq/L, default 1L).

Our goal is to reproduce the results of the PHREEQC script below.

DATABASE ../database/phreeqc.dat
TITLE Exchange in Dune sand

SOLUTION 1
pressure 1.0
temperature 25.0
Na 1.10006
Mg 0.48
Ca 1.9

EXCHANGE 1
X 0.06	#moles
-equilibrate 1

We start with setting up the chemical system and other means of modeling ion exchange in Reaktoro. For that, we use the PHREEQC database phreeqc.dat that contains ion exchange species, their corresponding reactions, and necessary parameters.

from reaktoro import *

db = PhreeqcDatabase("phreeqc.dat")

solution = AqueousPhase(speciate("H O C Ca Na Mg Cl"))
solution.set(ActivityModelPhreeqc(db))

The ion exchange phase is modelled by the IonExchangePhase class. We select only three ion exchange species: NaX, CaX2, and MgX2, which can be formed from major ions Na+, Ca2+, and Mg2+ in the dunewater. We use the Gaines-Thomas activity model.

exchange = IonExchangePhase("NaX CaX2 MgX2")
exchange.set(ActivityModelIonExchangeGainesThomas())

system = ChemicalSystem(db, solution, exchange)

# Define the equilibrium solver
solver = EquilibriumSolver(system)

Below, we define the initial chemical state for the ion exchange calculation involving dune sand and dunewater at 25 °C and 1 bar.

Note

We deliberately consider a very small initial amount of NaX, 0.06 μmol (and therefore there is 0.06 μmol of X- ion exchanger available for further formation of CaX2 and MgX2 during the equilibrium calculation). The rationale for this is to approximate the condition where there is an abundance of brine solution compared to the available ion exchanger, so that changes in fluid composition are negligible. This is necessary to compare our results with those produced by the PHREEQC script above, which ignores compositional changes in the aqueous solution when determining the distribution of NaX, CaX2 and MgX2 . Reaktoro, on the other hand, calculates not only this distribution, but also its effects on the brine composition, which we want to limit in this example.

Note

Unlike PHREEQC, to model ion exchange in Reaktoro we always assume that the exchanger (X-, in this case) is occupied by ions.

# Define initial equilibrium state
state = ChemicalState(system)
state.setTemperature(25, "celsius")
state.setPressure(1.0, "bar")
# Scale solution recipe to match the values of the PHREEQC examples
state.set("H2O"   , 1.00, "kg")
state.set("Na+" , 1.10, "mmol")
state.set("Mg+2", 0.48, "mmol")
state.set("Ca+2", 1.90, "mmol")
# Set the number of exchange assuming that it is completely occupied by sodium
state.set("NaX" , 0.06, "umol")

We expect the following set of ion exchange reactions happening between dunewater and exchanger site X- (initially fully occupied with sodium ion):

\[\begin{split} \begin{alignat}{4} \tfrac{1}{2} {\rm Ca^{+2}} + {\rm NaX} & \rightleftharpoons {\rm Na}^+ + \tfrac{1}{2} {\rm CaX}_2\\ \tfrac{1}{2} {\rm Mg^{+2}} + {\rm NaX} & \rightleftharpoons {\rm Na}^+ + \tfrac{1}{2} {\rm MgX}_2\\ \end{alignat} \end{split}\]

After the initial chemical state is defined, we can create the equilibrium solver that will perform a Gibbs energy minimization calculation to find the equilibrium state:

solver = EquilibriumSolver(system)
solver.solve(state)

By printing the instance state, we can inspect the content of calculated equilibrium state:

print(state)
+-----------------+------------+------+
| Property        |      Value | Unit |
+-----------------+------------+------+
| Temperature     |   298.1500 |    K |
| Pressure        |     1.0000 |  bar |
| Charge:         | 5.8600e-03 |  mol |
| Element Amount: |            |      |
| :: X            | 6.0000e-08 |  mol |
| :: H            | 1.1101e+02 |  mol |
| :: C            | 1.2000e-15 |  mol |
| :: O            | 5.5506e+01 |  mol |
| :: Na           | 1.1001e-03 |  mol |
| :: Mg           | 4.8000e-04 |  mol |
| :: Cl           | 1.0000e-16 |  mol |
| :: Ca           | 1.9000e-03 |  mol |
| Species Amount: |            |      |
| :: CO3-2        | 1.0000e-16 |  mol |
| :: H+           | 1.1624e-07 |  mol |
| :: H2O          | 5.5506e+01 |  mol |
| :: CO2          | 1.0000e-16 |  mol |
| :: (CO2)2       | 1.0000e-16 |  mol |
| :: HCO3-        | 1.0000e-16 |  mol |
| :: CH4          | 1.0000e-16 |  mol |
| :: Ca+2         | 1.9000e-03 |  mol |
| :: CaCO3        | 1.0000e-16 |  mol |
| :: CaHCO3+      | 1.0000e-16 |  mol |
| :: CaOH+        | 2.3220e-09 |  mol |
| :: Cl-          | 1.0000e-16 |  mol |
| :: H2           | 1.0000e-16 |  mol |
| :: Mg+2         | 4.7998e-04 |  mol |
| :: MgCO3        | 1.0000e-16 |  mol |
| :: MgHCO3+      | 1.0000e-16 |  mol |
| :: MgOH+        | 1.2831e-08 |  mol |
| :: Na+          | 1.1001e-03 |  mol |
| :: NaCO3-       | 1.0000e-16 |  mol |
| :: NaHCO3       | 1.0000e-16 |  mol |
| :: OH-          | 1.0109e-07 |  mol |
| :: NaOH         | 1.0000e-16 |  mol |
| :: O2           | 1.0000e-16 |  mol |
| :: NaX          | 5.5714e-10 |  mol |
| :: CaX2         | 2.5635e-08 |  mol |
| :: MgX2         | 4.0861e-09 |  mol |
+-----------------+------------+------+

We can now compute the properties of the brine solution after the ion exchange has been carried out. For this purpose, the AqueousProps class is used.

aqprops = AqueousProps(state)
print(f"I  = {aqprops.ionicStrength()} mol/kgw")
print(f"pH = {aqprops.pH()}")
print(f"pE = {aqprops.pE()}")
I  = 0.0053100562282740816 mol/kgw
pH = 6.965204999912109
pE = -3.9972562061620898

To inspect the ion exchange properties of species NaX, CaX2, and MgX2, the class IonExchangeProps can be used, which is done next:

exprops = IonExchangeProps(state)
print(exprops)
+--------------------------------------+------------+------+
| Property                             |      Value | Unit |
+--------------------------------------+------------+------+
| Element Amounts:                     |            |      |
| :: X                                 | 6.0000e-08 |  mol |
| :: Na                                | 5.5714e-10 |  mol |
| :: Mg                                | 4.0861e-09 |  mol |
| :: Ca                                | 2.5635e-08 |  mol |
| Species Amounts:                     |            |      |
| :: NaX                               | 5.5714e-10 |  mol |
| :: CaX2                              | 2.5635e-08 |  mol |
| :: MgX2                              | 4.0861e-09 |  mol |
| Equivalents:                         |            |      |
| :: NaX                               | 5.5714e-10 |   eq |
| :: CaX2                              | 5.1271e-08 |   eq |
| :: MgX2                              | 8.1723e-09 |   eq |
| Equivalent Fractions:                |            |      |
| :: NaX                               |     0.0093 |      |
| :: CaX2                              |     0.8545 |      |
| :: MgX2                              |     0.1362 |      |
| Activity Coefficients (log base 10): |            |      |
| :: NaX                               |    -0.0335 |      |
| :: CaX2                              |    -0.1322 |      |
| :: MgX2                              |    -0.1306 |      |
+--------------------------------------+------------+------+

We see that CaX2 has highest equivalent fraction of all ion exchange species, which causes its highest mole amounts. The table also includes the information on activity coefficient in logarithmic format with teh base 10. Such a detailed summary on the ion exchange properties was designed to compare PHREEQC outputs with those produced by Reaktoro.

To finalise, we have learned in this tutorial how to calculate the cation exchange complex in equilibrium with groundwater using the Gaines-Thomas framework. Similar studies can be conducted for different types of soil and alternative brines (such as seawater). The next tutorial, considers or particular a NaCl extractant to calculate the exchangeable cations in dune sand.