Chemical equilibrium with fixed pH
Chemical equilibrium with fixed pH#
Written by Allan Leal (ETH Zurich) on Jan 7th, 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!
Fixing the pH of an aqueous solution is a common procedure in biogeochemical modeling. This requires, however, the introduction of a new degree of freedom in the problem. This is because we cannot, in general, obtain a desired pH without titrating a substance in the solution. Thus, in the example given below, the system will be open to H+. The unknown amount of H+ necessary to enter or leave the system to attain the requested pH is the new degree of freedom, and is calculated along with the amounts of species in equilibrium.
Note
The choice of H+ as the titrant when fixing the pH of an aqueous solution may seem arbitrary and unintuitive to some coming from a different background or experienced with another convention used by other chemical modeling codes. For example, it can be questioned that a charged species not found isolated in nature (compared to HCl, for example) is chosen as the titrant. However, there are mathematical reasons behind this choice that make the problem to be solved more simply and efficiently. We will give examples later on how to combine this pH constraint with other conditions to demonstrate the full potential of this convention.
For our example, let’s create a ChemicalSystem object containing a simple aqueous phase (with thermodynamic data for the species fetched from a PhreeqcDatabase object):
from reaktoro import *
db = PhreeqcDatabase("pitzer.dat")
solution = AqueousPhase("H2O H+ OH- Na+ Cl- HCO3- CO2 CO3-2")
solution.set(ActivityModelPitzer())
system = ChemicalSystem(db, solution)
Let’s create a chemical equilibrium solver that expects the following properties as input conditions:
temperature,
pressure, and
pH.
The code below will first create an EquilibriumSpecs object in which the above properties are marked as inputs and then use that to create our EquilibriumSolver object:
specs = EquilibriumSpecs(system)
specs.temperature()
specs.pressure()
specs.pH()
solver = EquilibriumSolver(specs)
We create now an initial chemical state for our system (using ChemicalState), which will represent a 1 molal NaCl solution with 0.4 molal dissolved CO2:
state = ChemicalState(system)
state.temperature(30.0, "celsius")
state.pressure(1.0, "atm")
state.set("H2O", 1.0, "kg")
state.set("Na+", 1.0, "mol")
state.set("Cl-", 1.0, "mol")
state.set("CO2", 0.4, "mol")
print("INITIAL STATE")
print(state)
INITIAL STATE
+-----------------+-------------+------+
| Property | Value | Unit |
+-----------------+-------------+------+
| Temperature | 303.1500 | K |
| Pressure | 1.0132 | bar |
| Charge: | -3.0000e-16 | mol |
| Element Amount: | | |
| :: H | 1.1101e+02 | mol |
| :: C | 4.0000e-01 | mol |
| :: O | 5.6306e+01 | mol |
| :: Na | 1.0000e+00 | mol |
| :: Cl | 1.0000e+00 | mol |
| Species Amount: | | |
| :: H2O | 5.5506e+01 | mol |
| :: H+ | 1.0000e-16 | mol |
| :: OH- | 1.0000e-16 | mol |
| :: Na+ | 1.0000e+00 | mol |
| :: Cl- | 1.0000e+00 | mol |
| :: HCO3- | 1.0000e-16 | mol |
| :: CO2 | 4.0000e-01 | mol |
| :: CO3-2 | 1.0000e-16 | mol |
+-----------------+-------------+------+
We need now to specify the desired values for temperature, pressure, and pH in the equilibrium state we want to compute. This is done next using EquilibriumConditions:
conditions = EquilibriumConditions(specs)
conditions.temperature(50.0, "celsius")
conditions.pressure(10.0, "atm")
conditions.pH(1.0)
We are now ready for our chemical equilibrium calculation:
solver.solve(state, conditions)
print("FINAL STATE")
print(state)
FINAL STATE
+-----------------+------------+------+
| Property | Value | Unit |
+-----------------+------------+------+
| Temperature | 323.1500 | K |
| Pressure | 10.1325 | bar |
| Charge: | 1.2314e-01 | mol |
| Element Amount: | | |
| :: H | 1.1114e+02 | mol |
| :: C | 4.0000e-01 | mol |
| :: O | 5.6306e+01 | mol |
| :: Na | 1.0000e+00 | mol |
| :: Cl | 1.0000e+00 | mol |
| Species Amount: | | |
| :: H2O | 5.5506e+01 | mol |
| :: H+ | 1.2315e-01 | mol |
| :: OH- | 9.0556e-13 | mol |
| :: Na+ | 1.0000e+00 | mol |
| :: Cl- | 1.0000e+00 | mol |
| :: HCO3- | 4.9563e-06 | mol |
| :: CO2 | 4.0000e-01 | mol |
| :: CO3-2 | 1.8402e-14 | mol |
+-----------------+------------+------+
The printed state above does not show pH. We can obtain aqueous properties from a chemical state using AqueousProps:
aprops = AqueousProps(state)
print("AQUEOUS PROPERTIES AT EQUILIBRIUM")
print(aprops)
AQUEOUS PROPERTIES AT EQUILIBRIUM
+---------------------------------+------------+-------+
| Property | Value | Unit |
+---------------------------------+------------+-------+
| Temperature | 323.1500 | K |
| Pressure | 10.1325 | bar |
| Ionic Strength (Effective) | 1.0616 | molal |
| Ionic Strength (Stoichiometric) | 1.0616 | molal |
| pH | 1.0000 | |
| pE | -5.1887 | |
| Eh | -0.3327 | V |
| Charge Molality | 1.2314e-01 | molal |
| Element Molality: | | |
| :: C | 4.0000e-01 | molal |
| :: Na | 1.0000e+00 | molal |
| :: Cl | 1.0000e+00 | molal |
| Species Molality: | | |
| :: H+ | 1.2315e-01 | molal |
| :: OH- | 9.0556e-13 | molal |
| :: Na+ | 1.0000e+00 | molal |
| :: Cl- | 1.0000e+00 | molal |
| :: HCO3- | 4.9563e-06 | molal |
| :: CO2 | 4.0000e-01 | molal |
| :: CO3-2 | 1.8402e-14 | molal |
| Saturation Indices: | | |
| :: CO2(g) | 0.3683 | - |
| :: H2O(g) | -1.9384 | - |
| :: Halite :: NaCl | -1.9418 | - |
| :: Nahcolite :: NaHCO3 | -5.2359 | - |
| :: Natron :: Na2CO3:10H2O | -14.4740 | - |
| :: Trona :: Na3H(CO3)2:2H2O | -19.7685 | - |
| Saturation Ratios: | | |
| :: CO2(g) | 2.3353e+00 | - |
| :: H2O(g) | 1.1524e-02 | - |
| :: Halite :: NaCl | 1.1433e-02 | - |
| :: Nahcolite :: NaHCO3 | 5.8085e-06 | - |
| :: Natron :: Na2CO3:10H2O | 3.3574e-15 | - |
| :: Trona :: Na3H(CO3)2:2H2O | 1.7042e-20 | - |
+---------------------------------+------------+-------+
In the table above, we see that temperature, pressure, and pH have values exactly like those given earlier in the EquilibriumConditions object (although with different temperature and pressure units).
As a result of enforcing pH via specs.pH(), this executed equilibrium calculation considered the system to be open to H+. Let’s find out below how much H+ entered/leaved the system so that the desired pH could be reached:
print("Amount of H+ that entered the system:", state.equilibrium().titrantAmount("H+"), "mol")
Amount of H+ that entered the system: 0.123145 mol
Since H+ is a charged species, our final chemical state does not have zero charge as the initial state:
print("Charge at equilibrium state:", state.charge(), "mol")
Charge at equilibrium state: 0.123145 mol
We’ll learn in the next tutorial how a pH constraint can be combined with a charge constraint so that the system has zero charge at equilibrium (or whatever other value you may want to impose).