Creating chemical states

Creating chemical states

Written by Allan Leal (ETH Zurich) on Jan 4, 2022. Last revised on Feb 3, 2022.

After learning how to define chemical systems and create a ChemicalSystem object, it is time to create some chemical states. In your applications, you will most often deal with a single ChemicalSystem object and multiple ChemicalState objects in which the system’s temperature, pressure, and species amounts are specified.

Let’s check the following simple example, in which a chemical system is constructed with a gaseous phase only.

from reaktoro import *

db = NasaDatabase("nasa-cea")

gases = GaseousPhase("CO2 O2 H2 H2O CH4 CO")

system = ChemicalSystem(db, gases)

We can now create a ChemicalState object to represent the state of the system:

state = ChemicalState(system)
state.temperature(1000, "celsius")
state.pressure(10, "MPa")
state.set("CO2", 0.1, "mol")
state.set("O2" , 0.2, "mol")
state.set("H2" , 0.3, "mol")
state.set("H2O", 0.4, "mol")
state.set("CH4", 0.5, "mol")
state.set("CO" , 0.6, "mol")

We can print this chemical state to have an idea of what information it contains:

print(state)
+-----------------+------------+------+
| Property        |      Value | Unit |
+-----------------+------------+------+
| Temperature     |  1273.1500 |    K |
| Pressure        |   100.0000 |  bar |
| Charge:         | 0.0000e+00 |  mol |
| Element Amount: |            |      |
| :: H            | 3.4000e+00 |  mol |
| :: C            | 1.2000e+00 |  mol |
| :: O            | 1.6000e+00 |  mol |
| Species Amount: |            |      |
| :: CO2          | 1.0000e-01 |  mol |
| :: O2           | 2.0000e-01 |  mol |
| :: H2           | 3.0000e-01 |  mol |
| :: H2O          | 4.0000e-01 |  mol |
| :: CH4          | 5.0000e-01 |  mol |
| :: CO           | 6.0000e-01 |  mol |
+-----------------+------------+------+

Tip

Access the link ChemicalState to find out all methods available in class ChemicalState!

From the table above we see that ChemicalState objects do not have thermodynamic properties of the system; only temperature, pressure, species amounts, element amounts, and the electric charge. We will learn in the next guide how to evaluate thermodynamic and chemical properties of the system such as Gibbs energy, enthalpy, volume, heat capacities, species activities, etc.

Note

You will often use class ChemicalState for:

  • specifying initial conditions for chemical equilibrium and kinetics calculations

  • storing computed states from the above computations (e.g., chemical equilibrium states)

We will learn all these use cases in subsequent guides.

Let’s consider now a more complicated chemical system and create a chemical state for it:

db = PhreeqcDatabase("pitzer.dat")

solution = AqueousPhase()
minerals = MineralPhases("Halite Calcite Dolomite Quartz")

system = ChemicalSystem(db, solution, minerals)

state = ChemicalState(system)
state.temperature(50.0, "celsius")
state.pressure(5.0, "atm")
state.set("H2O"     , 1.00, "kg")
state.set("Na+"     , 0.10, "mol")
state.set("Cl-"     , 0.10, "mol")
state.set("Ca+2"    , 0.50, "mmol")
state.set("Mg+2"    , 0.50, "mmol")
state.set("Calcite" , 1.00, "mg")
state.set("Dolomite", 1.00, "ug")
state.set("Quartz"  , 1.00, "kmol")

print(state)
+---------------------------+------------+------+
| Property                  |      Value | Unit |
+---------------------------+------------+------+
| Temperature               |   323.1500 |    K |
| Pressure                  |     5.0663 |  bar |
| Charge:                   | 2.0000e-03 |  mol |
| Element Amount:           |            |      |
| :: H                      | 1.1101e+02 |  mol |
| :: C                      | 1.0002e-05 |  mol |
| :: O                      | 2.0555e+03 |  mol |
| :: Na                     | 1.0000e-01 |  mol |
| :: Mg                     | 5.0001e-04 |  mol |
| :: Si                     | 1.0000e+03 |  mol |
| :: Cl                     | 1.0000e-01 |  mol |
| :: Ca                     | 5.1000e-04 |  mol |
| Species Amount:           |            |      |
| :: H+                     | 1.0000e-16 |  mol |
| :: H2O                    | 5.5506e+01 |  mol |
| :: CO3-2                  | 1.0000e-16 |  mol |
| :: CO2                    | 1.0000e-16 |  mol |
| :: Ca+2                   | 5.0000e-04 |  mol |
| :: Cl-                    | 1.0000e-01 |  mol |
| :: H4SiO4                 | 1.0000e-16 |  mol |
| :: H2SiO4-2               | 1.0000e-16 |  mol |
| :: H3SiO4-                | 1.0000e-16 |  mol |
| :: HCO3-                  | 1.0000e-16 |  mol |
| :: Mg+2                   | 5.0000e-04 |  mol |
| :: MgCO3                  | 1.0000e-16 |  mol |
| :: MgOH+                  | 1.0000e-16 |  mol |
| :: Na+                    | 1.0000e-01 |  mol |
| :: OH-                    | 1.0000e-16 |  mol |
| :: Halite :: NaCl         | 1.0000e-16 |  mol |
| :: Calcite :: CaCO3       | 9.9909e-06 |  mol |
| :: Dolomite :: CaMg(CO3)2 | 5.4228e-09 |  mol |
| :: Quartz :: SiO2         | 1.0000e+03 |  mol |
+---------------------------+------------+------+

Tip

By default, a ChemicalState object is initialized with the following conditions:

  • 273.15 K (or 25 °C)

  • 105 Pa (or 1 bar)

  • 10-16 mol as the amounts of every species

Zero is a numerically problematic value for species amounts; that’s why we set it to a very small positive value instead. This should not be an issue for you in most cases, but if for some reason 10-16 is not small enough, you can use the method ChemicalState.setSpeciesAmount(1e-40) to set a common amount value for all species in the system.

Try: Create a ChemicalState object and print it before setting any property.

Let’s now set some surface areas between some phases. You’ll need to do this if you plan to model chemical kinetics in which heterogeneous reactions (i.e., reactions containing species from different phases, such as mineral and gas dissolution reactions) are considered.

state.surfaceArea("AqueousPhase", "Calcite", 1.2, "m2")
state.surfaceArea("AqueousPhase", "Quartz" , 1.4, "m2")

print(state)
+---------------------------+------------+------+
| Property                  |      Value | Unit |
+---------------------------+------------+------+
| Temperature               |   323.1500 |    K |
| Pressure                  |     5.0663 |  bar |
| Charge:                   | 2.0000e-03 |  mol |
| Element Amount:           |            |      |
| :: H                      | 1.1101e+02 |  mol |
| :: C                      | 1.0002e-05 |  mol |
| :: O                      | 2.0555e+03 |  mol |
| :: Na                     | 1.0000e-01 |  mol |
| :: Mg                     | 5.0001e-04 |  mol |
| :: Si                     | 1.0000e+03 |  mol |
| :: Cl                     | 1.0000e-01 |  mol |
| :: Ca                     | 5.1000e-04 |  mol |
| Species Amount:           |            |      |
| :: H+                     | 1.0000e-16 |  mol |
| :: H2O                    | 5.5506e+01 |  mol |
| :: CO3-2                  | 1.0000e-16 |  mol |
| :: CO2                    | 1.0000e-16 |  mol |
| :: Ca+2                   | 5.0000e-04 |  mol |
| :: Cl-                    | 1.0000e-01 |  mol |
| :: H4SiO4                 | 1.0000e-16 |  mol |
| :: H2SiO4-2               | 1.0000e-16 |  mol |
| :: H3SiO4-                | 1.0000e-16 |  mol |
| :: HCO3-                  | 1.0000e-16 |  mol |
| :: Mg+2                   | 5.0000e-04 |  mol |
| :: MgCO3                  | 1.0000e-16 |  mol |
| :: MgOH+                  | 1.0000e-16 |  mol |
| :: Na+                    | 1.0000e-01 |  mol |
| :: OH-                    | 1.0000e-16 |  mol |
| :: Halite :: NaCl         | 1.0000e-16 |  mol |
| :: Calcite :: CaCO3       | 9.9909e-06 |  mol |
| :: Dolomite :: CaMg(CO3)2 | 5.4228e-09 |  mol |
| :: Quartz :: SiO2         | 1.0000e+03 |  mol |
| Surface Area:             |            |      |
| :: AqueousPhase : Calcite |     1.2000 |   m2 |
| :: AqueousPhase : Quartz  |     1.4000 |   m2 |
+---------------------------+------------+------+

ChemicalState objects start with no surface areas. As you set surface areas for some phase pairs, interphase surface information is saved in the ChemicalState object. The code below demonstrates how to access the phase pairs for which surfaces areas have been set and retrieve their values:

print("EXISTING SURFACES IN THE CHEMICAL STATE OBJECT")
for [i, j] in state.surfaces():
    name1 = system.phase(i).name()      # the name of the first phase in the phase pair
    name2 = system.phase(j).name()      # the name of the other phase in the phase pair
    surfarea = state.surfaceArea(i, j)  # the surface area of this phase pair
    print(f"Surface area between {name1} and {name2}: {surfarea} m2")
EXISTING SURFACES IN THE CHEMICAL STATE OBJECT
Surface area between AqueousPhase and Calcite: 1.2 m2
Surface area between AqueousPhase and Quartz: 1.4 m2

That’s it. Now you know how to create ChemicalState objects and how to set their temperature, pressure, species amounts/mass, and surface areas between phases. Continue reading to learn more interesting use cases for ChemicalState, such as computing thermodynamic and chemical properties for the system using class ChemicalProps!