# Chemical equilibrium with fixed pH and charge balance¶

Written by Allan Leal (ETH Zurich) on Jan 7th, 2022

As we saw in the previous tutorial, imposing the pH of a solution in a chemical equilibrium calculation leaves the system open to H+. As a result, the final solution may not have zero electrical charge.

So, in this tutorial we demonstrate how to perform chemical equilibrium calculations with fixed pH values while simultaneously satisfying charge balance.

Let’s create a chemical system (using ChemicalSystem) with just a single aqueous phase containing a few selected species. For this phase, we use a Debye–Hückel activity model. The required thermodynamic data for the species is extracted from the SupcrtDatabase object:

from reaktoro import *

db = SupcrtDatabase("supcrtbl")

solution = AqueousPhase("H2O(aq) H+ OH- Na+ Cl- HCO3- CO2(aq) CO3-2")
solution.setActivityModel(ActivityModelDebyeHuckel())

system = ChemicalSystem(db, solution)


We want to perform a chemical equilibrium calculation in which the following properties are constrained:

• temperature,

• pressure,

• pH, and

• electric charge.

To attain a desired charge value in the solution, we will need to make the system open to a anion (negative ion) to counterbalance the amount of H+ entering/leaving the system. We will use Cl- for this.

Let’s now create a specialized and optimized EquilibriumSolver object to handle equilibrium problems with these specifications (note the use of EquilibriumSpecs below to achieve this):

specs = EquilibriumSpecs(system)
specs.temperature()
specs.pressure()
specs.pH()
specs.charge()
specs.openTo("Cl-")

solver = EquilibriumSolver(specs)


We should now create an initial chemical state for our chemical system:

state = ChemicalState(system)
state.temperature(30.0, "celsius")
state.pressure(1.0, "bar")
state.set("H2O(aq)", 1.00, "kg")
state.set("Na+"    , 0.01, "mol")
state.set("Cl-"    , 0.01, "mol")
state.set("CO2(aq)", 0.10, "mol")

print("INITIAL STATE")
print(state)

INITIAL STATE
+-----------------+-------------+------+
| Property        |       Value | Unit |
+-----------------+-------------+------+
| Temperature     |    303.1500 |    K |
| Pressure        |      1.0000 |  bar |
| Charge:         | -3.0000e-16 |  mol |
| Element Amount: |             |      |
| :: H            |  1.1102e+02 |  mol |
| :: C            |  1.0000e-01 |  mol |
| :: O            |  5.5708e+01 |  mol |
| :: Na           |  1.0000e-02 |  mol |
| :: Cl           |  1.0000e-02 |  mol |
| Species Amount: |             |      |
| :: H2O(aq)      |  5.5508e+01 |  mol |
| :: H+           |  1.0000e-16 |  mol |
| :: OH-          |  1.0000e-16 |  mol |
| :: Na+          |  1.0000e-02 |  mol |
| :: Cl-          |  1.0000e-02 |  mol |
| :: HCO3-        |  1.0000e-16 |  mol |
| :: CO2(aq)      |  1.0000e-01 |  mol |
| :: CO3-2        |  1.0000e-16 |  mol |
+-----------------+-------------+------+


Let’s now define the temperature, pressure, pH, and charge conditions we want to impose at equilibrium for this aqueous solution currently in disequilibrium. We do this by using an object of class EquilibriumConditions:

conditions = EquilibriumConditions(specs)
conditions.temperature(50.0, "celsius")
conditions.pressure(10.0, "bar")
conditions.pH(2.0)
conditions.charge(0.0)


We have everything we need now to perform the equilibrium calculation:

The code below will perform the equilibrium calculation and modify the state object to reflect the computed chemical equilibrium state:

result = solver.solve(state, conditions)

print("FINAL STATE")
print(state)

FINAL STATE
+-----------------+-------------+------+
| Property        |       Value | Unit |
+-----------------+-------------+------+
| Temperature     |    323.1500 |    K |
| Pressure        |     10.0000 |  bar |
| Charge:         | -2.2055e-18 |  mol |
| Element Amount: |             |      |
| :: H            |  1.1103e+02 |  mol |
| :: C            |  1.0000e-01 |  mol |
| :: O            |  5.5708e+01 |  mol |
| :: Na           |  1.0000e-02 |  mol |
| :: Cl           |  2.1127e-02 |  mol |
| Species Amount: |             |      |
| :: H2O(aq)      |  5.5508e+01 |  mol |
| :: H+           |  1.1133e-02 |  mol |
| :: OH-          |  6.1358e-12 |  mol |
| :: Na+          |  1.0000e-02 |  mol |
| :: Cl-          |  2.1127e-02 |  mol |
| :: HCO3-        |  6.1507e-06 |  mol |
| :: CO2(aq)      |  9.9994e-02 |  mol |
| :: CO3-2        |  6.0702e-14 |  mol |
+-----------------+-------------+------+


It’s always advisable to verify if the calculation succeeded:

print("Successful computation!" if result.optima.succeeded else "Computation has failed!")

Successful computation!


We can see from the table above that electric charge in the computed equilibrium state is numerically zero (i.e., very close to machine precision, 2.220446e-16).

Let’s check the aqueous properties of this equilibrium state to verify if we obtained the desired pH:

aprops = AqueousProps(state)

print("AQUEOUS PROPERTIES AT EQUILIBRIUM")
print(aprops)

AQUEOUS PROPERTIES AT EQUILIBRIUM
+-----------------------------------+------------+-------+
| Property                          |      Value |  Unit |
+-----------------------------------+------------+-------+
| Temperature                       |   323.1500 |     K |
| Pressure                          |    10.0000 |   bar |
| Ionic Strength (Effective)        |     0.0211 | molal |
| Ionic Strength (Stoichiometric)   |     0.0211 | molal |
| pH                                |     2.0000 |       |
| pE                                |     1.3421 |       |
| Eh                                |     0.0861 |     V |
| Element Molality:                 |            |       |
| :: C                              | 1.0000e-01 | molal |
| :: Na                             | 1.0000e-02 | molal |
| :: Cl                             | 2.1127e-02 | molal |
| Species Molality:                 |            |       |
| :: H+                             | 1.1133e-02 | molal |
| :: OH-                            | 6.1358e-12 | molal |
| :: Na+                            | 1.0000e-02 | molal |
| :: Cl-                            | 2.1127e-02 | molal |
| :: HCO3-                          | 6.1508e-06 | molal |
| :: CO2(aq)                        | 9.9994e-02 | molal |
| :: CO3-2                          | 6.0703e-14 | molal |
| Saturation Indices (log base 10): |            |       |
| :: O2(g)                          |   -64.0586 |     - |
| :: p-Cresol(g) :: CH3C6H4(OH)     |   -53.2650 |     - |
| :: m-Cresol(g) :: CH3C6H4(OH)     |   -51.6682 |     - |
| :: CH4(g)                         |    -5.4132 |     - |
| :: CO(g)                          |    -9.9715 |     - |
| :: CO2(g)                         |    -0.2822 |     - |
| :: H2(g)                          |    -7.1497 |     - |
| :: Ethylene(g) :: C2H4            |   -24.4885 |     - |
| :: H2O(g)                         |    -1.9187 |     - |
| :: Phenol(g) :: C6H5OH            |   -46.3912 |     - |
| :: o-Cresol(g) :: CH3C6H4(OH)     |   -52.2246 |     - |
| :: Graphite :: C                  |     0.0238 |     - |
| :: Diamond :: C                   |    -0.4748 |     - |
| :: Halite :: NaCl                 |    -5.3824 |     - |
+-----------------------------------+------------+-------+


Great! We can now see that pH is exactly 2, as we enforced, and temperature and pressure are 50 °C and 10 bar (but displayed in the table as 323.15 K and 106 Pa respectively).

Let’s now check how much H+ and Cl- entered/leaved the system (compared to the initial state) so that both pH and charge conditions could be reached:

print("Amount of H+ that entered the system:", float(state.equilibrium().implicitTitrantAmounts()), "mol")
print("Amount of Cl- that entered the system:", float(state.equilibrium().explicitTitrantAmounts()), "mol")

Amount of H+ that entered the system: 0.011127035386385114 mol
Amount of Cl- that entered the system: 0.011127035386384815 mol


That’s it! The results shown above demonstrate that our constraints were satisfied at equilibrium and that both H+ and Cl- entered the system with a numerically equal amount.