Defining materials

Written by Allan Leal (ETH Zurich) on January 31, 2022

When performing chemical equilibrium or chemical kinetics calculations, an initial chemical state is required. From this state, often in disequilibrium, you are interested in:

  1. immediately find a chemical equilibrium state that satisfies the desired constraints; or

  2. compute a sequence of states describing the evolution of the system as it reacts over a period of time.

In most cases, you should be able to construct an initial state in terms of given amounts of chemical species in the system using the ChemicalState class. However, this may not always be convenient or possible for you. In this tutorial, we show an alternative way to define an initial chemical state using the Material class. First, however, we give you an example that shows some disadvantages of working only with the ChemicalState class as a motivation to use Material.

An initial chemical state not so conveniently defined

Consider an aqueous solution with the following molal composition: 1 molal NaCl, 0.1 molal CaCl2, 0.05 molal MgCl2 and 0.1 molal CO2. Consider a rock with the following mass composition: 80%kg quartz (SiO2) and 20%kg calcite (CaCO3). Mix 1 kg of this aqueous solution with 10 g of rock and find the chemical equilibrium state of the system (i.e., solution and minerals) at 60°C and 15 bar.

As usual, we will start with the definition of a suitable chemical system to model this problem, which is done next.

from reaktoro import *

db = PhreeqcDatabase("phreeqc.dat")

solution = AqueousPhase()
solution.setActivityModel(ActivityModelPitzerHMW())

minerals = MineralPhases("Halite Calcite Dolomite Quartz")

system = ChemicalSystem(db, solution, minerals)

Note

In this system, the aqueous species were selected automatically based on the chemical elements composing the explicitly listed minerals (i.e., all aqueous species with elements H, O, Na, Cl, Ca, Mg, and Si). If you are curious, copy the code above and print the species in the chemical system!

for species in system.species():
    print(species.name()) 

Let’s now create our initial chemical containing 1 kg of aqueous solution and 10 g of rock:

state = ChemicalState(system)

# Set initial brine composition
state.set("H2O" , 1.00, "kg")
state.set("Na+" , 1.00, "mol")
state.set("Ca+2", 0.10, "mol")
state.set("Mg+2", 0.05, "mol")
state.set("Cl-" , 1.30, "mol")
state.set("CO2" , 0.10, "mol")

# Set initial rock composition
state.set("Quartz",  80.0, "g")
state.set("Calcite", 20.0, "g")

# Scale fluid and solid masses to 1 kg and 10 g respectively
state.scaleFluidMass(1.0, "kg")
state.scaleSolidMass(10.0, "g")

While defining this initial chemical state, note that:

  • we needed to specify amounts for ions Na+, Ca+2, and Mg+2 that reflect the desired molality of NaCl, CaCl2, and MgCl2.

  • we needed to specify an amount for Cl- that produces zero electric charge in the solution.

  • we needed to provide species names exactly how they exist in the database instead of chemical formulas (e.g., names Quartz and Calcite instead of formulas SiO2 and CaCO3).

  • we needed to scale the mass of fluid and solids in the chemical state.

Let’s now equilibrate this state at 60 °C and 15 bar:

state.temperature(60.0, "celsius")
state.pressure(15.0, "bar")

equilibrate(state)

print(state)
+---------------------------+-------------+------+
| Property                  |       Value | Unit |
+---------------------------+-------------+------+
| Temperature               |    333.1500 |    K |
| Pressure                  |     15.0000 |  bar |
| Charge:                   | -4.0936e-16 |  mol |
| Element Amount:           |             |      |
| :: H                      |  1.0291e+02 |  mol |
| :: C                      |  1.1269e-01 |  mol |
| :: O                      |  5.1968e+01 |  mol |
| :: Na                     |  9.2704e-01 |  mol |
| :: Mg                     |  4.6352e-02 |  mol |
| :: Si                     |  1.3315e-01 |  mol |
| :: Cl                     |  1.2052e+00 |  mol |
| :: Ca                     |  1.1269e-01 |  mol |
| Species Amount:           |             |      |
| :: CO3-2                  |  2.3074e-06 |  mol |
| :: H+                     |  6.6772e-06 |  mol |
| :: H2O                    |  5.1444e+01 |  mol |
| :: CO2                    |  8.1172e-02 |  mol |
| :: (CO2)2                 |  1.0000e-16 |  mol |
| :: HCO3-                  |  1.2608e-02 |  mol |
| :: CH4                    |  1.0000e-16 |  mol |
| :: Ca+2                   |  9.8700e-02 |  mol |
| :: CaCO3                  |  7.4077e-06 |  mol |
| :: CaHCO3+                |  5.5478e-03 |  mol |
| :: CaOH+                  |  1.0217e-09 |  mol |
| :: Cl-                    |  1.2052e+00 |  mol |
| :: H2                     |  1.0000e-16 |  mol |
| :: H4SiO4                 |  2.5791e-04 |  mol |
| :: H2SiO4-2               |  7.4502e-14 |  mol |
| :: H3SiO4-                |  4.5289e-08 |  mol |
| :: Mg+2                   |  4.3873e-02 |  mol |
| :: MgCO3                  |  1.3439e-06 |  mol |
| :: MgHCO3+                |  2.4778e-03 |  mol |
| :: MgOH+                  |  8.7426e-08 |  mol |
| :: Na+                    |  9.2460e-01 |  mol |
| :: NaCO3-                 |  1.4745e-05 |  mol |
| :: NaHCO3                 |  2.4242e-03 |  mol |
| :: OH-                    |  3.2735e-08 |  mol |
| :: NaOH                   |  1.0000e-16 |  mol |
| :: O2                     |  1.0000e-16 |  mol |
| :: Halite :: NaCl         |  1.0000e-16 |  mol |
| :: Calcite :: CaCO3       |  8.4303e-03 |  mol |
| :: Dolomite :: CaMg(CO3)2 |  1.1814e-14 |  mol |
| :: Quartz :: SiO2         |  1.3289e-01 |  mol |
+---------------------------+-------------+------+

We managed to get the problem defined and solved, but the coding process could be simpler as shown in the next section.

Simplifying with class Material

Let’s now use class Material to obtain the same result with more convenience:

brine = Material(system)
brine.add("H2O"  , 1.00, "kg")
brine.add("NaCl" , 1.00, "mol")
brine.add("CaCl2", 0.10, "mol")
brine.add("MgCl2", 0.05, "mol")
brine.add("CO2"  , 0.10, "mol")

rock = Material(system)
rock.add("SiO2",  80.0, "g")
rock.add("CaCO3", 20.0, "g")

mix = brine(1.0, "kg") + rock(10.0, "g")

state = mix.equilibrate(60, "celsius", 15, "bar")

print(state)
+---------------------------+-------------+------+
| Property                  |       Value | Unit |
+---------------------------+-------------+------+
| Temperature               |    333.1500 |    K |
| Pressure                  |     15.0000 |  bar |
| Charge:                   | -8.9543e-18 |  mol |
| Element Amount:           |             |      |
| :: H                      |  1.0291e+02 |  mol |
| :: C                      |  1.1269e-01 |  mol |
| :: O                      |  5.1968e+01 |  mol |
| :: Na                     |  9.2704e-01 |  mol |
| :: Mg                     |  4.6352e-02 |  mol |
| :: Si                     |  1.3315e-01 |  mol |
| :: Cl                     |  1.2052e+00 |  mol |
| :: Ca                     |  1.1269e-01 |  mol |
| Species Amount:           |             |      |
| :: CO3-2                  |  2.3074e-06 |  mol |
| :: H+                     |  6.6772e-06 |  mol |
| :: H2O                    |  5.1444e+01 |  mol |
| :: CO2                    |  8.1172e-02 |  mol |
| :: (CO2)2                 |  1.0000e-16 |  mol |
| :: HCO3-                  |  1.2608e-02 |  mol |
| :: CH4                    |  1.0000e-16 |  mol |
| :: Ca+2                   |  9.8700e-02 |  mol |
| :: CaCO3                  |  7.4077e-06 |  mol |
| :: CaHCO3+                |  5.5478e-03 |  mol |
| :: CaOH+                  |  1.0217e-09 |  mol |
| :: Cl-                    |  1.2052e+00 |  mol |
| :: H2                     |  1.0000e-16 |  mol |
| :: H4SiO4                 |  2.5791e-04 |  mol |
| :: H2SiO4-2               |  7.4502e-14 |  mol |
| :: H3SiO4-                |  4.5289e-08 |  mol |
| :: Mg+2                   |  4.3873e-02 |  mol |
| :: MgCO3                  |  1.3439e-06 |  mol |
| :: MgHCO3+                |  2.4778e-03 |  mol |
| :: MgOH+                  |  8.7426e-08 |  mol |
| :: Na+                    |  9.2460e-01 |  mol |
| :: NaCO3-                 |  1.4745e-05 |  mol |
| :: NaHCO3                 |  2.4242e-03 |  mol |
| :: OH-                    |  3.2736e-08 |  mol |
| :: NaOH                   |  1.0000e-16 |  mol |
| :: O2                     |  1.0000e-16 |  mol |
| :: Halite :: NaCl         |  1.0000e-16 |  mol |
| :: Calcite :: CaCO3       |  8.4303e-03 |  mol |
| :: Dolomite :: CaMg(CO3)2 |  1.1814e-14 |  mol |
| :: Quartz :: SiO2         |  1.3289e-01 |  mol |
+---------------------------+-------------+------+

The code above should (hopefully!) be intuitive enough and self-explanatory. The table printed above should also be identical (or extremely close!) to the previous one.

That’s it; you now have an extra Reaktoro component in your arsenal when performing your chemical modeling. The Material class should provide that extra level of convenience if things get a little more confusing if only the ChemicalState class is used.

Attention

You may have to use the Material class more carefully if chemical kinetics is considered. In this case, it is critical that species undergoing chemical reactions controlled by kinetic rates must have their initial amounts explicitly specified in the ChemicalState object, using their names instead of their formulas to avoid ambiguity. For example, suppose the minerals calcite and quartz are to react kinetically with the aqueous solution (and the aqueous species react instantaneously at equilibrium). In that case, the initial amounts of these minerals must be given in the ChemicalState object. The code below should give you an idea of how to proceed in this scenario:

brine = Material(system)
brine.add("H2O"  , 1.00, "kg")
brine.add("NaCl" , 1.00, "mol")
brine.add("CaCl2", 0.10, "mol")
brine.add("MgCl2", 0.05, "mol")
brine.add("CO2"  , 0.10, "mol")

state = brine(1.0, "kg").equilibrate(60, "celsius", 15, "bar")

state.set("Quartz",  80.0, "g")
state.set("Calcite", 20.0, "g")
state.scaleSolidMass(10.0, "g")

states = react(state, 10.0, "minutes")

Note: chemical kinetics (available in Reaktoro v1) is still in development in Reaktoro v2 (a completely new kinetics solver with modeling capabilities as extensive as the new equilibrium solver in Reaktoro v2). Thus, the react function above is still not available! If kinetics is already available in v2 and we forgot to remove this note, please report to us in Gitter.