Solubility of CO₂ in NaCl brines

In this tutorial, we show how Reaktoro can be used to compute the solubility of CO₂ in a 1 molal NaCl brine at temperature 60 °C and pressure 100 bar. We show no magical function to perform such calculation in a single line of code, but instead a sequence of steps using Reaktoro’s components (classes, methods) to enrich your understanding of how Reaktoro can be used for solving this and many other different chemical reaction modeling problems.

To calculate the solubility of CO₂ in the NaCl brine, we need two phases in chemical system: an aqueous phase to represent our NaCl brine, and a gaseous phase for CO₂ gas. Next, we formulate and solve a chemical equilibrium problem, in which 10 mol of CO₂ is mixed with 1 kg of H₂O and 1 mol of NaCl, at 60°C and 100 bar. The solution to this problem is a chemical equilibrium state, from which we can inspect how much CO₂ exists in the gaseous phase and compare this with the initial amount of CO₂ we used to mix with H₂O and NaCl.

Note

Our 1 molal NaCl brine is here represented by the mixture of 1 kg of H₂O and 1 mol of NaCl.

Note

Given that the mole mass of CO₂ is roughly 44 g/mol, 10 mol of CO₂ is approximately 440 g! Mixing this amount of CO₂ with 1 kg of H₂O will most likely result in a chemical equilibrium state in which the aqueous phase is saturated with CO₂ and the remaining CO₂ exists as a gas (or super-critical fluid depending on the temperature and pressure). If we choose a small amount for CO₂, we can end up with an equilibrium state in which all CO₂ has been dissolved in the aqueous phase, and we will then not be able to determine its solubility.

We present below the Python script that performs a multi-phase, multi-species chemical equilibrium calculation using Reaktoro to determine the solubility of CO₂ in a 1 molal NaCl brine at temperature 60 °C and pressure 100 bar.

# Step 1: Import the reaktoro Python package
from reaktoro import *

# Step 2: Initialize a thermodynamic database
db = Database('supcrt98.xml')

# Step 3: Define the chemical system
editor = ChemicalEditor(db)
editor.addAqueousPhaseWithElements('H O Na Cl C')
editor.addGaseousPhase(['CO2(g)'])

# Step 4: Construct the chemical system
system = ChemicalSystem(editor)

# Step 5: Define the chemical equilibrium problem
problem = EquilibriumProblem(system)
problem.setTemperature(60, 'celsius')
problem.setPressure(100, 'bar')
problem.add('H2O', 1.0, 'kg')
problem.add('NaCl', 1.0, 'mol')
problem.add('CO2', 10.0, 'mol')

# Step 6: Calculate the chemical equilibrium state
state = equilibrate(problem)

# Step 7: Output the calculated chemical state to a file
state.output('result.txt')

# Step 8: Print the amounts of some aqueous species
print('Amount of CO2(aq):', state.speciesAmount('CO2(aq)'))
print('Amount of HCO3-:', state.speciesAmount('HCO3-'))
print('Amount of CO3--:', state.speciesAmount('CO3--'))

# Step 9: Print the amounts of element C in both aqueous and gaseous phases
print('Amount of C in aqueous phase:', state.elementAmountInPhase('C', 'Aqueous'))
print('Amount of C in gaseous phase:', state.elementAmountInPhase('C', 'Gaseous'))

You find next a step-by-step explanation of the above script.

Importing the reaktoro Python package

# Step 1: Import the reaktoro Python package
from reaktoro import *

Using Reaktoro in Python requires first an import of the python package reaktoro. From this point on, we are able to use the library components of Reaktoro (classes, methods, constants), which are needed to define our chemical system and chemical reaction modeling problems.

Note

To simplify the tutorials, we use from reaktoro import *, which imports all components of the reaktoro package into the default Python namespace, which can potentially create name conflicts. For your applications, consider instead using import reaktoro as rkt, and then refer to Reaktoro’s classes and methods as rkt.Database, rkt.ChemicalSystem, rkt.equilibrate, and so forth.

Initializing a thermodynamic database

Thermodynamic databases are essential for modeling chemically reactive systems using Reaktoro. We need a database from where we collect data of substances that will compose our phases of interest in a multi-phase chemical system. A thermodynamic database also contains model parameters for the evaluation of standard thermodynamic properties of species and/or reactions (e.g., standard Gibbs energies, equilibrium constants).

In this step:

# Step 2: Initialize a thermodynamic database
db = Database('supcrt98.xml')

we initialize a Database object with the supcrt98.xml database file. This database was generated from the original SUPCRT92 database file slop98.dat. You are welcome to inspect these files and learn more about the chemical species available in them. You can also read more about the available thermodynamic databases supported in Reaktoro at Thermodynamic Databases.

Note

All SUPCRT92 thermodynamic databases have been embedded into Reaktoro. Thus, you don’t actually need to have a database file named supcrt98.xml in a local directory when initializing the Database object. If you want to use a customized database file, however, also named supcrt98.xml, then your local file will be used instead.

Tip

If you are using a customized version of a thermodynamic database, consider changing its name (e.g., custom-supcrt98.xml) to avoid accidental use of an embedded database. This can happen if you do not give a correct path to your custom database file.

Defining the chemical system

Reaktoro is a general-purpose chemical solver that avoids as much as possible presuming specific assumptions about your problems. Thus, you need to specify how your chemical system should be defined. This encompasses the specification of all phases in the system as well as the chemical species that compose each phase. By using the ChemicalEditor class, you can conveniently achieve this as shown below:

# Step 3: Define the chemical system
editor = ChemicalEditor(db)
editor.addAqueousPhaseWithElements('H O Na Cl C')
editor.addGaseousPhase(['CO2(g)'])

In this step, we create an object of class ChemicalEditor and specify that an aqueous phase and a gaseous phase should be considered in the chemical system. The aqueous phase is defined such that its species are all those aqueous species in the selected thermodynamic database that can be created by combining the chemical elements H, O, Na, Cl, and C. The gaseous phase is defined with only one gaseous species: CO₂(g).

Note

The selected elements H, O, Na, Cl, and C in the definition of the aqueous phase and the choice of CO₂(g) as the single gaseous species composing the gaseous phase consistently represent our intentions of calculating the solubility of CO₂ in a NaCl brine. If we decide to do a similar computation, but with a NaCl-MgCl₂-CaCl₂ brine, then that original list of elements would need to be incremented with Mg and Ca.

Note

An automatic search for chemical species can result in a large number of species in the phase, potentially causing the chemical reaction calculations to be more computationally expensive. If you are using Reaktoro in demanding applications (e.g., as a chemical solver in a reactive transport simulator), you might want to manually specify the chemical species of each phase in your chemical system. This can be achieved by providing a list of species names as shown below:

editor.addAqueousPhaseWithElements([
    'H2O(l)',
    'H+',
    'OH-',
    'Na+',
    'Cl-',
    'HCO3-',
    'CO3--',
    'CO2(aq)'
    ])

This is exactly what we did for the definition of the gaseous phase. If we had done instead:

editor.addGaseousPhase('C O')

then other gases would be considered, such as CO(g) and O2(g), which are not of interest in our modeling problem.

Caution

If you manually specify the chemical species in a phase, you need to make sure that they exist in the thermodynamic database with the exact same name! Replacing 'CO2(g)' above with 'CO2' will cause an error if the database has no gaseous species with such a name.

Note

By default, activities of the aqueous species are calculated using the HKF extended Debye-Hückel model for solvent water and ionic species, except for the aqueous species CO₂ (aq), for which the Drummond model is used. For gases, Peng-Robinson equation of state is used. For different activity models (e.g., Pitzer), consider inspecting the documentation of classes AqueousPhase and GaseousPhase.

Constructing the chemical system

# Step 4: Construct the chemical system
system = ChemicalSystem(editor)

This step is where we create an object of class ChemicalSystem using the chemical system definition details stored in the object editor.

Note

ChemicalSystem is perhaps the main class in Reaktoro. An object of this class stores the phases, species and elements in our defined chemical system, as well as provides the means to compute many types of thermodynamic properties, such as standard thermodynamic properties (e.g., standard Gibbs energies, standard enthalpies, and standard volumes of species), and thermo-chemical properties (e.g., activity and activity coefficients of species; density, enthalpy and internal energy of phases). As you learn more about other Reaktoro’s classes, you will note that an object of class ChemicalSystem is almost always needed for their initialization!

Defining the chemical equilibrium problem

We have now defined and constructed our chemical system of interest, enabling us to move on to the next step in Reaktoro’s modeling workflow: defining our chemical reaction problems. Below we create an equilibrium problem with our prescribed equilibrium conditions for temperature, pressure, and amounts of elements that are consistent with our intention of calculating the solubility of CO₂ at 60 °C and 100 bar in a 1 molal NaCl brine.

# Step 5: Define the chemical equilibrium problem
problem = EquilibriumProblem(system)
problem.setTemperature(60, 'celsius')
problem.setPressure(100, 'bar')
problem.add('H2O', 1.0, 'kg')
problem.add('NaCl', 1.0, 'mol')
problem.add('CO2', 10.0, 'mol')

Note

Did you pay attention we said prescribed equilibrium conditions for the amounts of elements? Since we actually provided the amounts of substances H₂O, NaCl, and CO₂ in the above code, this statement seems a little bit confusing at least. Here is what happens behind the scenes: Reaktoro parses these chemical formulas and determines the elements and their coefficients. Once this is done, the amount of each element stored inside the object problem is incremented according to the given amount of substance and its coefficient in the formula. The amounts of elements you provide are then used as constraints for the Gibbs energy minimization calculation when computing the state of chemical equilibrium (i.e., when we try to find the amounts of all species in the system that corresponds to a state of minimum Gibbs energy and at the same time satisfying the element amounts constraints).

Danger

Now that you know that an equivalent chemical equilibrium problem could be defined with:

problem.add('H' 111.0, 'mol')
problem.add('O' 75.5, 'mol')
problem.add('Na' 1.0, 'mol')
problem.add('Cl' 1.0, 'mol')
problem.add('C' 10.0, 'mol')

assuming that 1 kg of H₂O is roughly 55.5 mol, you might want to adventure in manually specifying different values for the amounts of elements. Just be extra careful with the values you provide as this could accidentally result in an infeasible chemical equilibrium state.

Here is a simple example of element amount conditions that result in an infeasible equilibrium state. Consider a chemical system containing only a gaseous phase with gases H₂O(g) and CO₂(g). Find non-negative amounts for H₂O(g) and CO₂(g) when the given amounts of elements are: 2 mol of H, 1 mol of O, and 1 mol of C.

Note

Please note that we are not condemning the input form shown above in terms of element amounts, but only telling you to be attentive with the values you input. If you are using Reaktoro as a chemical reaction solver in a reactive transport simulator, for example, you’ll most likely need to work directly with given amounts of elements, which shows that this input form is required in certain cases. For such time-dependent modeling problems, you often only need to ensure that the initial conditions for elements amounts result in feasible initial species’ amounts.

Tip

The substance formulas given in the method add of class EquilibriumProblem can, but do not need to, correspond to names of chemical species in the thermodynamic database. Even unusual, if not strange, substance formulas, such as HCl3(NaO)4C13, would be understood by that method. We do not promise, however, that you will obtain a feasible chemical equilibrium state with unrealistic conditions!

Note

In Reaktoro, the word element is used as a synonym of components of chemical species, and not necessarily chemical elements. Electric charge, for example, is considered as an element, even though it is technically not. Thus, we say the ionic species CO₃^2- is composed of elements C, O, and Z, with coefficients 1, 3, and -2 respectively, where Z is the symbol we use to denote the electric charge element.

Warning

Prefer the use of neutral substances when using the method add of class EquilibriumProblem, unless you definitely need to add a charged, ionic species in the recipe. The following code:

problem.add('H+' 0.1, 'mmol')

will increment not only the amount of element H by 0.1 mmol, but also the electric charge element Z. As a result, the composition of the aqueous phase at equilibrium will not be electrically neutral, which might not be your intention. The following:

problem.add('H+' 0.1, 'mmol')
problem.add('Cl-' 0.1, 'mmol')

would result in the amount of element Z equal to zero.

Calculating the chemical equilibrium state

# Step 6: Calculate the chemical equilibrium state
state = equilibrate(problem)

In this step, we use the equilibrate function to calculate the chemical equilibrium state of the system with the given equilibrium conditions stored in the object problem. For this calculation, Reaktoro uses an efficient Gibbs energy minimization computation to determine the species amounts that correspond to a state of minimum Gibbs energy in the system, while satisfying the prescribed amount conditions for the temperature, pressure, and element amounts. The result is stored in the object state, of class ChemicalState.

Attention

In the vast majority of cases, you’ll only have one object of class ChemicalSystem in your code and one or more objects of class ChemicalState describing different states of your chemical system! Reaktoro differentiates these two independent concepts: chemical system definition and chemical system state.

Tip

The method equilibrate is a convenient function in Reaktoro. Consider using the class EquilibriumSolver for more advanced requirements. For example, if you have to perform many equilibrium calculations in sequence. The equilibrate method has a computational overhead because every call creates a new object of class EquilibriumSolver. Preferably, this object should be created only once, and then used subsequently for all other equilibrium calculations. Here is a demonstration:

solver = EquilibriumSolver(system) # Our chemical equilibrium solver
state = ChemicalState(system)  # Our chemical state

solver.solve(state, problem)  # Initial equilibrium calculation
state.output('state-initial.txt')  # Output the initial equilibrium state

problem.add('NaCl', 0.1, 'mol')  # Increment the amount of NaCl

solver.solve(state, problem)  # Subsequent equilibrium calculation
state.output('state-modified.txt')  # Output the modified equilibrium state

Outputting the calculated chemical state to a file

We have performed our chemical equilibrium calculation and now we want to inspect the computed compositional state and its thermodynamic properties, which can be done by outputting the chemical state to a file, here named result.txt.

# Step 7: Output the calculated chemical state to a file
state.output('result.txt')

Printing the amounts of some aqueous species

Here is just a small demonstration of getting species amount information from a ChemicalState object using the method speciesAmount of class ChemicalState to extract the amount of a few chemical species. Please inspect the API of ChemicalState class to learn more about its methods.

# Step 8: Print the amounts of some aqueous species
print('Amount of CO2(aq):', state.speciesAmount('CO2(aq)'))
print('Amount of HCO3-:', state.speciesAmount('HCO3-'))
print('Amount of CO3--:', state.speciesAmount('CO3--'))

Printing the amounts of element C in both aqueous and gaseous phases

Finally, we print the amounts of element C in both aqueous and gaseous phases:

# Step 9: Print the amounts of element C in both aqueous and gaseous phases
print('Amount of C in aqueous phase:', state.elementAmountInPhase('C', 'Aqueous'))
print('Amount of C in gaseous phase:', state.elementAmountInPhase('C', 'Gaseous'))

In this specific case, in which there was no initial element C in the aqueous phase, the value corresponding to the amount of element C in the aqueous phase is our solubility of CO₂ in the NaCl brine with the previously prescribed conditions.

Tip

If we had used 2 kg of H₂O, we would have needed to divide the calculated amount of element C in the aqueous phase by 2 to obtain the solubility in molal (mol per kg of H₂O).

Have you got an issue?

Have you found any issue or error in this tutorial? Do you have any recommendations or you think something is not clear enough? Please, let us know by filling a new issue here:

Reaktoro’s GitHub Issues

You’ll need a GitHub account - but this is easy to sort out if you don’t have one yet!